1 <!DOCTYPE Book PUBLIC "-//Davenport//DTD DocBook V3.0//EN">
4 <title>GiNaC MAJOR_VERSION.MINOR_VERSION Tutorial</title>
6 <subtitle>An open framework for symbolic computation within the C++ programming language</subtitle>
10 <collabname>The GiNaC Group</collabname>
13 <firstname>Christian</firstname><surname>Bauer</surname>
15 <address><email>Christian.Bauer@Uni-Mainz.DE</email></address>
19 <firstname>Alexander</firstname><surname>Frink</surname>
21 <address><email>Alexander.Frink@Uni-Mainz.DE</email></address>
25 <firstname>Richard</firstname><othername>B.</othername><surname>Kreckel</surname>
27 <address><email>Richard.Kreckel@Uni-Mainz.DE</email></address>
35 <title>Introduction</title>
37 <para>The motivation behind GiNaC derives from the observation that
38 most present day computer algebra systems (CAS) are linguistically and
39 semantically impoverished. It is an attempt to overcome the current
40 situation by extending a well established and standardized computer
41 language (C++) by some fundamental symbolic capabilities, thus
42 allowing for integrated systems that embed symbolic manipulations
43 together with more established areas of computer science (like
44 computation-intense numeric applications, graphical interfaces, etc.)
45 under one roof.</para>
47 <para>This tutorial is intended for the novice user who is new to
48 GiNaC but already has some background in C++ programming. However,
49 since a hand made documentation like this one is difficult to keep in
50 sync with the development the actual documentation is inside the
51 sources in the form of comments. That documentation may be parsed by
52 one of the many Javadoc-like documentation systems. If you fail at
53 generating it you may access it directly at URL <ulink
54 url="http://www.ginac.de/reference/"><literal>http://www.ginac.de/reference/</literal></ulink>.
55 It is an invaluable resource not only for the advanced user who wishes
56 to extend the system (or chase bugs) but for everybody who wants to
57 comprehend the inner workings of GiNaC. This little tutorial on the
58 other hand only covers the basic things that are unlikely to change in
59 the near future. </para>
61 <sect1><title>License</title>
63 <para>The GiNaC framework for symbolic computation within the C++
64 programming language is Copyright (C) 1999 Johannes Gutenberg
65 Universität Mainz, Germany.</para>
67 <para>This program is free software; you can redistribute it and/or
68 modify it under the terms of the GNU General Public License as
69 published by the Free Software Foundation; either version 2 of the
70 License, or (at your option) any later version.</para>
72 <para>This program is distributed in the hope that it will be useful, but
73 WITHOUT ANY WARRANTY; without even the implied warranty of
74 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
75 General Public License for more details.</para>
77 <para>You should have received a copy of the GNU General Public License
78 along with this program; see the file COPYING. If not, write to the
79 Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
80 MA 02111-1307, USA.</para>
85 <title>A Tour of GiNaC</title>
87 <para>This quick tour of GiNaC wants to rise your interest in in the
88 subsequent chapters by showing off a bit. Please excuse us if it
89 leaves many open questions.</para>
91 <sect1><title>How to use it from within C++</title> <para>The GiNaC
92 open framework for symbolic computation within the C++ programming
93 language does not try to define a language of it's own as conventional
94 CAS do. Instead, it extends the capabilities of C++ by symbolic
95 manipulations. Here is how to generate and print a simple (and
96 pointless) bivariate polynomial with some large coefficients:
98 <title>My first GiNaC program (a bivariate polynomial)</title>
100 #include <ginac/ginac.h>
101 using namespace GiNaC;
105 symbol x("x"), y("y");
108 for (int i=0; i<3; ++i)
109 poly += factorial(i+16)*pow(x,i)*pow(y,2-i);
111 cout << poly << endl;
115 <para>Assuming the file is called <literal>hello.cc</literal>, on
116 our system we can compile and run it like this:</para>
118 <prompt>$</prompt> c++ hello.cc -o hello -lcln -lginac
119 <prompt>$</prompt> ./hello
120 355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
125 <para>Next, there is a more meaningful C++ program that calls a
126 function which generates Hermite polynomials in a specified free
129 <title>My second GiNaC program (Hermite polynomials)</title>
131 #include <ginac/ginac.h>
132 using namespace GiNaC;
134 ex HermitePoly(symbol x, int deg)
136 ex HKer=exp(-pow(x,2));
137 // uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2)
138 return normal(pow(-1,deg) * diff(HKer, x, deg) / HKer);
145 for (int i=0; i<6; ++i)
146 cout << "H_" << i << "(z) == " << HermitePoly(z,i) << endl;
151 <para>When run, this will type out</para>
156 H_3(z) == -12*z+8*z^3
157 H_4(z) == -48*z^2+16*z^4+12
158 H_5(z) == 120*z-160*z^3+32*z^5
161 This method of generating the coefficients is of course far from
162 optimal for production purposes.</para>
164 <para>In order to show some more examples of what GiNaC can do we
165 will now use <literal>ginsh</literal>, a simple GiNaC interactive
166 shell that provides a convenient window into GiNaC's capabilities.
169 <sect1><title>What it can do for you</title>
171 <para>After invoking <literal>ginsh</literal> one can test and
172 experiment with GiNaC's features much like in other Computer Algebra
173 Systems except that it does not provide programming constructs like
174 loops or conditionals. For a concise description of the
175 <literal>ginsh</literal> syntax we refer to its accompanied man-page.
176 Suffice to say that assignments and comparisons in
177 <literal>ginsh</literal> are written as they are in C,
178 i.e. <literal>=</literal> assigns and <literal>==</literal>
181 <para>It can manipulate arbitrary precision integers in a very fast
182 way. Rational numbers are automatically converted to fractions of
186 369988485035126972924700782451696644186473100389722973815184405301748249
188 123329495011708990974900260817232214728824366796574324605061468433916083
196 <para>All numbers occuring in GiNaC's expressions can be converted
197 into floating point numbers with the <literal>evalf</literal> method,
198 to arbitrary accuracy:
201 0.14285714285714285714
205 0.1428571428571428571428571428571428571428571428571428571428571428571428
206 5714285714285714285714285714285714285
210 <para>Exact numbers other than rationals that can be manipulated in
211 GiNaC include predefined constants like Archimedes' Pi. They can both
212 be used in symbolic manipulations (as an exact number) as well as in
213 numeric expressions (as an inexact number):
218 x+9.869604401089358619L0
222 11.869604401089358619L0
226 <para>Built-in functions evaluate immediately to exact numbers if
227 this is possible. Conversions that can be safely performed are done
228 immediately; conversions that are not generally valid are not done:
237 (Note that converting the last input to <literal>x</literal> would
238 allow one to conclude that <literal>42*Pi</literal> is equal to
239 <literal>0</literal>.)</para>
241 <para>Linear equation systems can be solved along with basic linear
242 algebra manipulations over symbolic expressions. In C++ GiNaC offers
243 a matrix class for this purpose but we can see what it can do using
244 <literal>ginsh</literal>'s notation of double brackets to type them
247 > lsolve(a+x*y==z,x);
249 lsolve([3*x+5*y == 7, -2*x+10*y == -5], [x, y]);
251 > M = [[ [[1, 3]], [[-3, 2]] ]];
252 [[ [[1,3]], [[-3,2]] ]]
255 > charpoly(M,lambda);
260 <para>Multivariate polynomials and rational functions may be expanded,
261 collected and normalized (i.e. converted to a ratio of two coprime
264 > a = x^4 + 2*x^2*y^2 + 4*x^3*y + 12*x*y^3 - 3*y^4;
265 -3*y^4+x^4+12*x*y^3+2*x^2*y^2+4*x^3*y
266 > b = x^2 + 4*x*y - y^2;
269 3*y^6+x^6-24*x*y^5+43*x^2*y^4+16*x^3*y^3+17*x^4*y^2+8*x^5*y
271 3*y^6+48*x*y^4+2*x^2*y^2+x^4*(-y^2+x^2+4*x*y)+4*x^3*y*(-y^2+x^2+4*x*y)
277 <para>You can differentiate functions and expand them as Taylor or
278 Laurent series (the third argument of series is the evaluation point,
279 the fourth defines the order):
283 > series(sin(x),x,0,4);
285 > series(1/tan(x),x,0,4);
286 x^(-1)-1/3*x+Order(x^2)
290 <para>If you ever wanted to convert units in C or C++ and found this
291 is cumbersome, here is the solution. Symbolic types can always be
292 used as tags for different types of objects. Converting from wrong
293 units to the metric system is therefore easy:
300 140613.91592783185568*kg*m^(-2)
310 <title>Installation</title>
312 <para>GiNaC's installation follows the spirit of most GNU software. It is
313 easily installed on your system by three steps: configuration, build,
316 <sect1><title id="CLN-main">Prerequistes</title>
318 <para>In order to install GiNaC on your system, some prerequistes need
319 to be met. First of all, you need to have a C++-compiler adhering to
320 the ANSI-standard <citation>ISO/IEC 14882:1998(E)</citation>. We used
321 <literal>GCC</literal> for development so if you have a different
322 compiler you are on your own. For the configuration to succeed you
323 need a Posix compliant shell installed in <literal>/bin/sh</literal>,
324 GNU <literal>bash</literal> is fine. Perl is needed by the built
325 process as well, since some of the source files are automatically
326 generated by Perl scripts. Last but not least, Bruno Haible's library
327 <literal>CLN</literal> is extensively used and needs to be installed
328 on your system. Please get it from <ulink
329 url="ftp://ftp.santafe.edu/pub/gnu/"><literal>ftp://ftp.santafe.edu/pub/gnu/</literal></ulink>
331 url="ftp://ftp.ilog.fr/pub/Users/haible/gnu/"><literal>ftp://ftp.ilog.fr/pub/Users/haible/gnu/</literal></ulink>
332 (it is covered by GPL) and install it prior to trying to install
333 GiNaC. The configure script checks if it can find it and if it cannot
334 it will refuse to continue.</para></sect1>
336 <sect1><title>Configuration</title>
338 <para>To configure GiNaC means to prepare the source distribution for
339 building. It is done via a shell script called
340 <literal>configure</literal> that is shipped with the sources.
341 (Actually, this script is by itself created with GNU Autoconf from the
342 files <literal>configure.in</literal> and
343 <literal>aclocal.m4</literal>.) Since a configure script generated by
344 GNU Autoconf never prompts, all customization must be done either via
345 command line parameters or environment variables. It accepts a list
346 of parameters, the complete set of which can be listed by calling it
347 with the <literal>--help</literal> option. The most important ones
348 will be shortly described in what follows:
351 <para><literal>--disable-shared</literal>: When given, this option
352 switches off the build of a shared library, i.e. a
353 <literal>.so</literal>-file. This may be convenient when developing
354 because it considerably speeds up compilation.</para>
357 <para><literal>--prefix=</literal><emphasis>PREFIX</emphasis>: The
358 directory where the compiled library and headers are installed. It
359 defaults to <literal>/usr/local</literal> which means that the
360 library is installed in the directory
361 <literal>/usr/local/lib</literal> and the header files in
362 <literal>/usr/local/include/GiNaC</literal> and the documentation
363 (like this one) into <literal>/usr/local/share/doc/GiNaC</literal>.</para>
366 <para><literal>--libdir=</literal><emphasis>LIBDIR</emphasis>: Use
367 this option in case you want to have the library installed in some
369 <emphasis>PREFIX</emphasis><literal>/lib/</literal>.</para>
372 <para><literal>--includedir=</literal><emphasis>INCLUDEDIR</emphasis>:
373 Use this option in case you want to have the header files
374 installed in some other directory than
375 <emphasis>PREFIX</emphasis><literal>/include/ginac/</literal>. For
376 instance, if you specify
377 <literal>--includedir=/usr/include</literal> you will end up with
378 the header files sitting in the directory
379 <literal>/usr/include/ginac/</literal>. Note that the subdirectory
380 <literal>GiNaC</literal> is enforced by this process in order to
381 keep the header files separated from others. This avoids some
382 clashes and allows for an easier deinstallation of GiNaC. This ought
383 to be considered A Good Thing (tm).</para>
386 <para><literal>--datadir=</literal><emphasis>DATADIR</emphasis>:
387 This option may be given in case you want to have the documentation
388 installed in some other directory than
389 <emphasis>PREFIX</emphasis><literal>/share/doc/GiNaC/</literal>.
394 <para>In addition, you may specify some environment variables.
395 <literal>CXX</literal> holds the path and the name of the C++ compiler
396 in case you want to override the default in your path. (The
397 <literal>configure</literal> script searches your path for
398 <literal>c++</literal>, <literal>g++</literal>,
399 <literal>gcc</literal>, <literal>CC</literal>, <literal>cxx</literal>
400 and <literal>cc++</literal> in that order.) It may be very useful to
401 define some compiler flags with the <literal>CXXFLAGS</literal>
402 environment variable, like optimization, debugging information and
403 warning levels. If ommitted, it defaults to <literal>-g
404 -O2</literal>.</para>
406 <para>The whole process is illustrated in the following two
407 examples. (Substitute <literal>setenv VARIABLE value</literal> for
408 <literal>export VARIABLE=value</literal> if the Berkeley C shell is
411 <example><title>Sample sessions of how to call the
412 configure-script</title> <para>Simple configuration for a site-wide
413 GiNaC library assuming everything is in default paths:</para>
415 <prompt>$</prompt> export CXXFLAGS="-Wall -O2"
416 <prompt>$</prompt> ./configure
418 <para>Configuration for a private static GiNaC library with several
419 components sitting in custom places (site-wide <literal>GCC</literal>
420 and private <literal>CLN</literal>), the compiler pursueded to be
421 picky and full assertions switched on:</para>
423 <prompt>$</prompt> export CXX=/usr/local/gnu/bin/c++
424 <prompt>$</prompt> export CPPFLAGS="${CPPFLAGS} -I${HOME}/include"
425 <prompt>$</prompt> export CXXFLAGS="${CXXFLAGS} -DDO_GINAC_ASSERT -ggdb -Wall -ansi -pedantic -O2"
426 <prompt>$</prompt> export LDFLAGS="${LDFLAGS} -L${HOME}/lib"
427 <prompt>$</prompt> ./configure --disable-shared --prefix=${HOME}
434 <sect1><title>Building GiNaC</title>
436 <para>After proper configuration you should just build the whole
437 library by typing <literal>make</literal> at the command
438 prompt and go for a cup of coffee.</para>
440 <para>Just to make sure GiNaC works properly you may run a simple test
441 suite by typing <literal>make check</literal>. This will compile some
442 sample programs, run them and compare the output to reference output.
443 Each of the checks should return a message <literal>passed</literal>
444 together with the CPU time used for that particular test. If it does
445 not, something went wrong. This is mostly intended to be a QA-check
446 if something was broken during the development, but not a sanity check
447 of your system. Another intent is to allow people to fiddle around
448 with optimization. If <literal>CLN</literal> was installed all right
449 this step is unlikely to return any errors.</para>
453 <sect1><title>Installation</title>
455 <para>To install GiNaC on your system, simply type <literal>make
456 install</literal>. As described in the section about configuration
457 the files will be installed in the following directories (the
458 directories will be created if they don't already exist):
461 <para><literal>libginac.a</literal> will go into
462 <emphasis>PREFIX</emphasis><literal>/lib/</literal> (or
463 <emphasis>LIBDIR</emphasis>) which defaults to
464 <literal>/usr/local/lib/</literal>. So will
465 <literal>libginac.so</literal> if the the configure script was
466 given the option <literal>--enable-shared</literal>. In that
467 case, the proper symlinks will be established as well (by running
468 <literal>ldconfig</literal>).</para>
471 <para>All the header files will be installed into
472 <emphasis>PREFIX</emphasis><literal>/include/ginac/</literal> (or
473 <emphasis>INCLUDEDIR</emphasis><literal>/ginac/</literal>, if
477 <para>All documentation (HTML, Postscript and DVI) will be stuffed
479 <emphasis>PREFIX</emphasis><literal>/share/doc/GiNaC/</literal>
480 (or <emphasis>DATADIR</emphasis><literal>/doc/GiNaC/</literal>, if
486 <para>Just for the record we will list some other useful make targets:
487 <literal>make clean</literal> deletes all files generated by
488 <literal>make</literal>, i.e. all the object files. In addition
489 <literal>make distclean</literal> removes all files generated by
490 configuration. And finally <literal>make uninstall</literal> removes
491 the installed library and header files<footnoteref
492 linkend="warnuninstall-1">.
494 <footnote id="warnuninstall-1"><para>Uninstallation does not work
495 after you have called <literal>make distclean</literal> since the
496 <literal>Makefile</literal> is itself generated by the configuration
497 from <literal>Makefile.in</literal> and hence deleted by
498 <literal>make distclean</literal>. There are two obvious ways out
499 of this dilemma. First, you can run the configuration again with
500 the same <emphasis>PREFIX</emphasis> thus creating a
501 <literal>Makefile</literal> with a working
502 <literal>uninstall</literal> target. Second, you can do it by hand
503 since you now know where all the files went during
504 installation.</para></footnote>
512 <title>Basic Concepts</title>
514 <para>This chapter will describe the different fundamental objects
515 that can be handled with GiNaC. But before doing so, it is worthwhile
516 introducing you to the more commonly used class of expressions,
517 representing a flexible meta-class for storing all mathematical
520 <sect1><title>Expressions</title>
522 <para>The most common class of objects a user deals with is the
523 expression <literal>ex</literal>, representing a mathematical object
524 like a variable, number, function, sum, product, etc... Expressions
525 may be put together to form new expressions, passed as arguments to
526 functions, and so on. Here is a little collection of valid
528 <example><title>Examples of expressions</title>
530 ex MyEx1 = 5; // simple number
531 ex MyEx2 = x + 2*y; // polynomial in x and y
532 ex MyEx3 = (x + 1)/(x - 1); // rational expression
533 ex MyEx4 = sin(x + 2*y) + 3*z + 41; // containing a function
534 ex MyEx5 = MyEx4 + 1; // similar to above
537 Before describing the more fundamental objects that form the building
538 blocks of expressions we'll have a quick look under the hood by
539 describing how expressions are internally managed.</para>
541 <sect2><title>Digression: Expressions are reference counted</title>
543 <para>An expression is extremely light-weight since internally it
544 works like a handle to the actual representation and really holds
545 nothing more than a pointer to some other object. What this means in
546 practice is that whenever you create two <literal>ex</literal> and set
547 the second equal to the first no copying process is involved. Instead,
548 the copying takes place as soon as you try to change the second.
549 Consider the simple sequence of code:
550 <example><title>Simple copy-on-write semantics</title>
552 #include <ginac/ginac.h>
553 using namespace GiNaC;
557 symbol x("x"), y("y"), z("z");
560 e1 = sin(x + 2*y) + 3*z + 41;
561 e2 = e1; // e2 points to same object as e1
562 cout << e2 << endl; // prints sin(x+2*y)+3*z+41
563 e2 += 1; // e2 is copied into a new object
564 cout << e2 << endl; // prints sin(x+2*y)+3*z+42
569 The line <literal>e2 = e1;</literal> creates a second expression
570 pointing to the object held already by <literal>e1</literal>. The
571 time involved for this operation is therefore constant, no matter how
572 large <literal>e1</literal> was. Actual copying, however, must take
573 place in the line <literal>e2 += 1</literal> because
574 <literal>e1</literal> and <literal>e2</literal> are not handles for
575 the same object any more. This concept is called
576 <emphasis>copy-on-write semantics</emphasis>. It increases
577 performance considerably whenever one object occurs multiple times and
578 represents a simple garbage collection scheme because when an
579 <literal>ex</literal> runs out of scope its destructor checks whether
580 other expressions handle the object it points to too and deletes the
581 object from memory if that turns out not to be the case. A slightly
582 less trivial example of differentiation using the chain-rule should
583 make clear how powerful this can be. <example><title>Advanced
584 copy-on-write semantics</title>
586 #include <ginac/ginac.h>
587 using namespace GiNaC;
591 symbol x("x"), y("y");
595 ex e3 = diff(sin(e2), x); // first derivative of sin(e2) by x
596 cout << e1 << endl // prints x+3*y
597 << e2 << endl // prints (x+3*y)^3
598 << e3 << endl; // prints 3*(x+3*y)^2*cos((x+3*y)^3)
603 Here, <literal>e1</literal> will actually be referenced three times
604 while <literal>e2</literal> will be referenced two times. When the
605 power of an expression is built, that expression needs not be
606 copied. Likewise, since the derivative of a power of an expression can
607 be easily expressed in terms of that expression, no copying of
608 <literal>e1</literal> is involved when <literal>e3</literal> is
609 constructed. So, when <literal>e3</literal> is constructed it will
610 print as <literal>3*(x+3*y)^2*cos((x+3*y)^3)</literal> but the
611 argument of <literal>cos()</literal> only holds a reference to
612 <literal>e2</literal> and the factor in front is just
613 <literal>3*e1^2</literal>.
616 <para>As a user of GiNaC, you cannot see this mechanism of
617 copy-on-write semantics. When you insert an expression into a second
618 expression, the result behaves exactly as if the contents of the first
619 expression were inserted. But it may be useful to remember that this
620 is not what happens. Knowing this will enable you to write much more
621 efficient code. If you still have an uncertain feeling with
622 copy-on-write semantics, we recommend you have a look at the
623 <emphasis>C++-FAQ lite</emphasis> by Marshall Cline. Chapter 16
624 covers this issue and presents an implementation which is pretty close
625 to the one in GiNaC. You can find it on the Web at <ulink
626 url="http://www.cerfnet.com/~mpcline/c++-faq-lite/"><literal>http://www.cerfnet.com/~mpcline/c++-faq-lite/</literal></ulink>.</para>
628 <para>So much for expressions. But what exactly are these expressions
629 handles of? This will be answered in the following sections.</para>
633 <sect1><title>The Class Hierarchy</title>
635 <para>GiNaC's class hierarchy consists of several classes representing
636 mathematical objects, all of which (except for <literal>ex</literal>
637 and some helpers) are internally derived from one abstract base class
638 called <literal>basic</literal>. You do not have to deal with objects
639 of class <literal>basic</literal>, instead you'll be dealing with
640 symbols and functions of symbols. You'll soon learn in this chapter
641 how many of the functions on symbols are really classes. This is
642 because simple symbolic arithmetic is not supported by languages like
643 C++ so in a certain way GiNaC has to implement its own arithmetic.</para>
645 <para>To give an idea about what kinds of symbolic composits may be
646 built we have a look at the most important classes in the class
647 hierarchy. The dashed line symbolizes a "points to" or "handles"
648 relationship while the solid lines stand for "inherits from"
650 <figure id="classhier-id" float="1">
651 <title>The GiNaC class hierarchy</title>
652 <graphic align="center" fileref="classhierarchy.graext" format="GRAEXT"></graphic>
654 Some of the classes shown here (the ones sitting in white boxes) are
655 abstract base classes that are of no interest at all for the user.
656 They are used internally in order to avoid code duplication if
657 two or more classes derived from them share certain features. An
658 example would be <literal>expairseq</literal>, which is a container
659 for a sequence of pairs each consisting of one expression and a number
660 (<literal>numeric</literal>). What <emphasis>is</emphasis> visible to
661 the user are the derived classes <literal>add</literal> and
662 <literal>mul</literal>, representing sums of terms and products,
663 respectively. We'll come back later to some more details about these
664 two classes and motivate the use of pairs in sums and products here.</para>
666 <sect2><title>Digression: Internal representation of products and sums</title>
668 <para>Although it should be completely transparent for the user of
669 GiNaC a short discussion of this topic helps to understand the sources
670 and also explain performance to a large degree. Consider the symbolic
672 <emphasis>2*d<superscript>3</superscript>*(4*a+5*b-3)</emphasis>,
673 which could naively be represented by a tree of linear containers for
674 addition and multiplication, one container for exponentiation with
675 base and exponent and some atomic leaves of symbols and numbers in
677 <figure id="repres-naive-id" float="1">
678 <title>Naive internal representation-tree for <emphasis>2*d<superscript>3</superscript>*(4*a+5*b-3)</emphasis></title>
679 <graphic align="center" fileref="rep_naive.graext" format="GRAEXT"></graphic>
681 However, doing so results in a rather deeply nested tree which will
682 quickly become inefficient to manipulate. If we represent the sum
683 instead as a sequence of terms, each having a purely numeric
684 multiplicative coefficient and the multiplication as a sequence of
685 terms, each having a numeric exponent, the tree becomes much more
687 <figure id="repres-pair-id" float="1">
688 <title>Pair-wise internal representation-tree for <emphasis>2*d<superscript>3</superscript>*(4*a+5*b-3)</emphasis></title>
689 <graphic align="center" fileref="rep_pair.graext" format="GRAEXT"></graphic>
691 The number <literal>3</literal> above the symbol <literal>d</literal>
692 shows that <literal>mul</literal> objects are treated similarly where
693 the coefficients are interpreted as <emphasis>exponents</emphasis>
694 now. Addition of sums of terms or multiplication of products with
695 numerical exponents can be coded to be very efficient with such a
696 pair-representation. Internally, this handling is done by many CAS in
697 this way. It typically speeds up manipulations by an order of
698 magnitude. The overall multiplicative factor <literal>2</literal> and
699 the additive term <literal>-3</literal> look somewhat cumbersome in
700 this representation, however, since they are still carrying a trivial
701 exponent and multiplicative factor <literal>1</literal> respectively.
702 Within GiNaC, this is avoided by adding a field that carries overall
704 <figure id="repres-real-id" float="1">
705 <title>Realistic picture of GiNaC's representation-tree for <emphasis>2*d<superscript>3</superscript>*(4*a+5*b-3)</emphasis></title>
706 <graphic align="center" fileref="rep_real.graext" format="GRAEXT"></graphic>
708 This also allows for a better handling of numeric radicals, since
709 <literal>sqrt(2)</literal> can now be carried along calculations. Now
710 it should be clear, why both classes <literal>add</literal> and
711 <literal>mul</literal> are derived from the same abstract class: the
712 data representation is the same, only the semantics differs. In the
713 class hierarchy, methods for polynomial expansion and such are
714 reimplemented for <literal>add</literal> and <literal>mul</literal>,
715 but the data structure is inherited from
716 <literal>expairseq</literal>.</para>
720 <sect1><title>Symbols</title>
722 <para>Symbols are for symbolic manipulation what atoms are for
723 chemistry. You can declare objects of class <literal>symbol</literal>
724 as any other object simply by saying <literal>symbol x,y;</literal>.
725 There is, however, a catch in here having to do with the fact that C++
726 is a compiled language. The information about the symbol's name is
727 thrown away by the compiler but at a later stage you may want to print
728 expressions holding your symbols. In order to avoid confusion GiNaC's
729 symbols are able to know their own name. This is accomplished by
730 declaring its name for output at construction time in the fashion
731 <literal>symbol x("x");</literal>. If you declare a symbol using the
732 default constructor (i.e. without string-argument) the system will
733 deal out a unique name. That name may not be suitable for printing
734 but for internal routines when no output is desired it is often
735 enough. We'll come across examples of such symbols later in this
738 <para>This implies that the stings passed to symbols at construction
739 time may not be used for comparing two of them. It is perfectly
740 legitimate to write <literal>symbol x("x"),y("x");</literal> but it is
741 likely to lead into trouble. Here, <literal>x</literal> and
742 <literal>y</literal> are different symbols and statements like
743 <literal>x-y</literal> will not be simplified to zero although the
744 output <literal>x-x</literal> looks funny. Such output may also occur
745 when there are two different symbols in two scopes, for instance when
746 you call a function that declares a symbol with a name already
747 existent in a symbol in the calling function. Again, comparing them
748 (using <literal>operator==</literal> for instance) will always reveal
749 their difference. Watch out, please.</para>
751 <para>Although symbols can be assigned expressions for internal
752 reasons, you should not do it (and we are not going to tell you how it
753 is done). If you want to replace a symbol with something else in an
754 expression, you can use the expression's <literal>.subs()</literal>
759 <sect1><title>Numbers</title>
761 <para>For storing numerical things, GiNaC uses Bruno Haible's library
762 <literal>CLN</literal>. The classes therein serve as foundation
763 classes for GiNaC. <literal>CLN</literal> stands for Class Library
764 for Numbers or alternatively for Common Lisp Numbers. In order to
765 find out more about <literal>CLN</literal>'s internals the reader is
766 refered to the documentation of that library. Suffice to say that it
767 is by itself build on top of another library, the GNU Multiple
768 Precision library <literal>GMP</literal>, which is an extremely fast
769 library for arbitrary long integers and rationals as well as arbitrary
770 precision floating point numbers. It is very commonly used by several
771 popular cryptographic applications. <literal>CLN</literal> extends
772 <literal>GMP</literal> by several useful things: First, it introduces
773 the complex number field over either reals (i.e. floating point
774 numbers with arbitrary precision) or rationals. Second, it
775 automatically converts rationals to integers if the denominator is
776 unity and complex numbers to real numbers if the imaginary part
777 vanishes and also correctly treats algebraic functions. Third it
778 provides good implementations of state-of-the-art algorithms for all
779 trigonometric and hyperbolic functions as well as for calculation of
780 some useful constants.</para>
782 <para>The user can construct an object of class
783 <literal>numeric</literal> in several ways. The following example
784 shows the four most important constructors: construction from
785 C-integer, construction of fractions from two integers, construction
786 from C-float and construction from a string.
787 <example><title>Construction of numbers</title>
789 #include <ginac/ginac.h>
790 using namespace GiNaC;
794 numeric two(2); // exact integer 2
795 numeric r(2,3); // exact fraction 2/3
796 numeric e(2.71828); // floating point number
797 numeric p("3.1415926535897932385"); // floating point number
799 cout << two*p << endl; // floating point 6.283...
804 Note that all those constructors are <emphasis>explicit</emphasis>
805 which means you are not allowed to write <literal>numeric
806 two=2;</literal>. This is because the basic objects to be handled by
807 GiNaC are the expressions <literal>ex</literal> and we want to keep
808 things simple and wish objects like <literal>pow(x,2)</literal> to be
809 handled the same way as <literal>pow(x,a)</literal>, which means that
810 we need to allow a general <literal>ex</literal> as base and exponent.
811 Therefore there is an implicit constructor from C-integers directly to
812 expressions handling numerics at work in most of our examples. This
813 design really becomes convenient when one declares own functions
814 having more than one parameter but it forbids using implicit
815 constructors because that would lead to ambiguities. </para>
817 <para>It may be tempting to construct numbers writing <literal>numeric
818 r(3/2)</literal>. This would, however, call C's built-in operator
819 <literal>/</literal> for integers first and result in a numeric
820 holding a plain integer 1. <emphasis>Never use
821 </emphasis><literal>/</literal><emphasis> on integers!</emphasis> Use
822 the constructor from two integers instead, as shown in the example
823 above. Writing <literal>numeric(1)/2</literal> may look funny but
826 <para>We have seen now the distinction between exact numbers and
827 floating point numbers. Clearly, the user should never have to worry
828 about dynamically created exact numbers, since their "exactness"
829 always determines how they ought to be handled. The situation is
830 different for floating point numbers. Their accuracy is handled by
831 one <emphasis>global</emphasis> variable, called
832 <literal>Digits</literal>. (For those readers who know about Maple:
833 it behaves very much like Maple's <literal>Digits</literal>). All
834 objects of class numeric that are constructed from then on will be
835 stored with a precision matching that number of decimal digits:
836 <example><title>Controlling the precision of floating point numbers</title>
838 #include <ginac/ginac.h>
839 using namespace GiNaC;
843 numeric three(3.0), one(1.0);
844 numeric x = one/three;
846 cout << "in " << Digits << " digits:" << endl;
847 cout << x << endl;
848 cout << Pi.evalf() << endl;
860 The above example prints the following output to screen:
866 0.333333333333333333333333333333333333333333333333333333333333333333
867 3.14159265358979323846264338327950288419716939937510582097494459231
871 <para>It should be clear that objects of class
872 <literal>numeric</literal> should be used for constructing numbers or
873 for doing arithmetic with them. The objects one deals with most of
874 the time are the polymorphic expressions <literal>ex</literal>.</para>
876 <sect2><title>Tests on numbers</title>
878 <para>Once you have declared some numbers, assigned them to
879 expressions and done some arithmetic with them it is frequently
880 desired to retrieve some kind of information from them like asking
881 whether that number is integer, rational, real or complex. For those
882 cases GiNaC provides several useful methods. (Internally, they fall
883 back to invocations of certain CLN functions.)</para>
885 <para>As an example, let's construct some rational number, multiply it
886 with some multiple of its denominator and check what comes out:
887 <example><title>Sample test on objects of type numeric</title>
889 #include <ginac/ginac.h>
890 using namespace GiNaC;
892 // some very important constants:
893 const numeric twentyone(21);
894 const numeric ten(10);
895 const numeric fife(5);
899 numeric answer = twentyone;
902 cout << answer.is_integer() << endl; // false, it's 21/5
904 cout << answer.is_integer() << endl; // true, it's 42 now!
910 Note that the variable <literal>answer</literal> is constructed here
911 as an integer by <literal>numeric</literal>'s copy constructor but in
912 an intermediate step it holds a rational number represented as integer
913 numerator and integer denominator. When multiplied by 10, the
914 denominator becomes unity and the result is automatically converted to
915 a pure integer again. Internally, the underlying
916 <literal>CLN</literal> is responsible for this behaviour and we refer
917 the reader to <literal>CLN</literal>'s documentation. Suffice to say
918 that the same behaviour applies to complex numbers as well as return
919 values of certain functions. Complex numbers are automatically
920 converted to real numbers if the imaginary part becomes zero. The
921 full set of tests that can be applied is listed in the following
924 <informaltable colsep="0" frame="topbot" pgwide="1">
926 <colspec colnum="1" colwidth="1*">
927 <colspec colnum="2" colwidth="2*">
930 <entry>Method</entry>
931 <entry>Returns true if...</entry>
936 <entry><literal>.is_zero()</literal></entry>
937 <entry>object is equal to zero</entry>
940 <entry><literal>.is_positive()</literal></entry>
941 <entry>object is not complex and greater than 0</entry>
944 <entry><literal>.is_integer()</literal></entry>
945 <entry>object is a (non-complex) integer</entry>
948 <entry><literal>.is_pos_integer()</literal></entry>
949 <entry>object is an integer and greater than 0</entry>
952 <entry><literal>.is_nonneg_integer()</literal></entry>
953 <entry>object is an integer and greater equal 0</entry>
956 <entry><literal>.is_even()</literal></entry>
957 <entry>object is an even integer</entry>
960 <entry><literal>.is_odd()</literal></entry>
961 <entry>object is an odd integer</entry>
964 <entry><literal>.is_prime()</literal></entry>
965 <entry>object is a prime integer (probabilistic primality test)</entry>
968 <entry><literal>.is_rational()</literal></entry>
969 <entry>object is an exact rational number (integers are rational, too, as are complex extensions like <literal>2/3+7/2*I</literal>)</entry>
972 <entry><literal>.is_real()</literal></entry>
973 <entry>object is a real integer, rational or float (i.e. is not complex)</entry>
985 <sect1><title>Constants</title>
987 <para>Constants behave pretty much like symbols except that that they return
988 some specific number when the method <literal>.evalf()</literal> is called.
991 <para>The predefined known constants are:
992 <informaltable colsep="0" frame="topbot" pgwide="1">
994 <colspec colnum="1" colwidth="1*">
995 <colspec colnum="2" colwidth="2*">
996 <colspec colnum="3" colwidth="4*">
1000 <entry>Common Name</entry>
1001 <entry>Numerical Value (35 digits)</entry>
1006 <entry><literal>Pi</literal></entry>
1007 <entry>Archimedes' constant</entry>
1008 <entry>3.14159265358979323846264338327950288</entry>
1010 <entry><literal>Catalan</literal></entry>
1011 <entry>Catalan's constant</entry>
1012 <entry>0.91596559417721901505460351493238411</entry>
1014 <entry><literal>EulerGamma</literal></entry>
1015 <entry>Euler's (or Euler-Mascheroni) constant</entry>
1016 <entry>0.57721566490153286060651209008240243</entry>
1025 <sect1><title>Fundamental operations: The <literal>power</literal>, <literal>add</literal> and <literal>mul</literal> classes</title>
1027 <para>Simple polynomial expressions are written down in GiNaC pretty
1028 much like in other CAS. The necessary operators <literal>+</literal>,
1029 <literal>-</literal>, <literal>*</literal> and <literal>/</literal>
1030 have been overloaded to achieve this goal. When you run the following
1031 program, the constructor for an object of type <literal>mul</literal>
1032 is automatically called to hold the product of <literal>a</literal>
1033 and <literal>b</literal> and then the constructor for an object of
1034 type <literal>add</literal> is called to hold the sum of that
1035 <literal>mul</literal> object and the number one:
1036 <example><title>Construction of <literal>add</literal> and <literal>mul</literal> objects</title>
1038 #include <ginac/ginac.h>
1039 using namespace GiNaC;
1043 symbol a("a"), b("b");
1050 <para>For exponentiation, you have already seen the somewhat clumsy
1051 (though C-ish) statement <literal>pow(x,2);</literal> to represent
1052 <literal>x</literal> squared. This direct construction is necessary
1053 since we cannot safely overload the constructor <literal>^</literal>
1054 in <literal>C++</literal> to construct a <literal>power</literal>
1055 object. If we did, it would have several counterintuitive effects:
1058 <para>Due to <literal>C</literal>'s operator precedence,
1059 <literal>2*x^2</literal> would be parsed as <literal>(2*x)^2</literal>.
1062 <para>Due to the binding of the operator <literal>^</literal>,
1063 <literal>x^a^b</literal> would result in <literal>(x^a)^b</literal>.
1064 This would be confusing since most (though not all) other CAS
1065 interpret this as <literal>x^(a^b)</literal>.
1068 <para>Also, expressions involving integer exponents are very
1069 frequently used, which makes it even more dangerous to overload
1070 <literal>^</literal> since it is then hard to distinguish between the
1071 semantics as exponentiation and the one for exclusive or. (It would
1072 be embarassing to return <literal>1</literal> where one has requested
1073 <literal>2^3</literal>.)</para>
1076 All effects are contrary to mathematical notation and differ from the
1077 way most other CAS handle exponentiation, therefore overloading
1078 <literal>^</literal> is ruled out for GiNaC's C++ part. The situation
1079 is different in <literal>ginsh</literal>, there the
1080 exponentiation-<literal>^</literal> exists. (Also note, that the
1081 other frequently used exponentiation operator <literal>**</literal>
1082 does not exist at all in <literal>C++</literal>).</para>
1084 <para>To be somewhat more precise, objects of the three classes
1085 described here, are all containers for other expressions. An object
1086 of class <literal>power</literal> is best viewed as a container with
1087 two slots, one for the basis, one for the exponent. All valid GiNaC
1088 expressions can be inserted. However, basic transformations like
1089 simplifying <literal>pow(pow(x,2),3)</literal> to
1090 <literal>x^6</literal> automatically are only performed when
1091 this is mathematically possible. If we replace the outer exponent
1092 three in the example by some symbols <literal>a</literal>, the
1093 simplification is not safe and will not be performed, since
1094 <literal>a</literal> might be <literal>1/2</literal> and
1095 <literal>x</literal> negative.</para>
1097 <para>Objects of type <literal>add</literal> and
1098 <literal>mul</literal> are containers with an arbitrary number of
1099 slots for expressions to be inserted. Again, simple and safe
1100 simplifications are carried out like transforming
1101 <literal>3*x+4-x</literal> to <literal>2*x+4</literal>.</para>
1103 <para>The general rule is that when you construct such objects, GiNaC
1104 automatically creates them in canonical form, which might differ from
1105 the form you typed in your program. This allows for rapid comparison
1106 of expressions, since after all <literal>a-a</literal> is simply zero.
1107 Note, that the canonical form is not necessarily lexicographical
1108 ordering or in any way easily guessable. It is only guaranteed that
1109 constructing the same expression twice, either implicitly or
1110 explicitly, results in the same canonical form.</para>
1114 <sect1><title>Built-in Functions</title>
1116 <para>There are quite a number of useful functions built into GiNaC.
1117 They are all objects of class <literal>function</literal>. They
1118 accept one or more expressions as arguments and return one expression.
1119 If the arguments are not numerical, the evaluation of the functions
1120 may be halted, as it does in the next example:
1121 <example><title>Evaluation of built-in functions</title>
1123 #include <ginac/ginac.h>
1124 using namespace GiNaC;
1128 symbol x("x"), y("y");
1131 cout << "gamma(" << foo << ") -> " << gamma(foo) << endl;
1132 ex bar = foo.subs(y==1);
1133 cout << "gamma(" << bar << ") -> " << gamma(bar) << endl;
1134 ex foobar= bar.subs(x==7);
1135 cout << "gamma(" << foobar << ") -> " << gamma(foobar) << endl;
1139 <para>This program will type out two times a function and then an
1140 expression that may be really useful:
1142 gamma(x+(1/2)*y) -> gamma(x+(1/2)*y)
1143 gamma(x+1/2) -> gamma(x+1/2)
1144 gamma(15/2) -> (135135/128)*Pi^(1/2)
1147 Most of these functions can be differentiated, series expanded so on.
1148 Read the next chapter in order to learn more about this..</para>
1156 <title>Important Algorithms</title>
1158 <para>In this chapter the most important algorithms provided by GiNaC
1159 will be described. Some of them are implemented as functions on
1160 expressions, others are implemented as methods provided by expression
1161 objects. If they are methods, there exists a wrapper function around
1162 it, so you can alternatively call it in a functional way as shown in
1164 <example><title>Methods vs. wrapper functions</title>
1166 #include <ginac/ginac.h>
1167 using namespace GiNaC;
1171 ex x = numeric(1.0);
1173 cout << "As method: " << sin(x).evalf() << endl;
1174 cout << "As function: " << evalf(sin(x)) << endl;
1179 The general rule is that wherever methods accept one or more
1180 parameters (<emphasis>arg1</emphasis>, <emphasis>arg2</emphasis>, ...)
1181 the order of arguments the function wrapper accepts is the same but
1182 preceded by the object to act on (<emphasis>object</emphasis>,
1183 <emphasis>arg1</emphasis>, <emphasis>arg2</emphasis>, ...). This
1184 approach is the most natural one in an OO model but it may lead to
1185 confusion for MapleV users because where they would type
1186 <literal>A:=x+1; subs(x=2,A);</literal> GiNaC would require
1187 <literal>A=x+1; subs(A,x==2);</literal> (after proper declaration of A
1188 and x). On the other hand, since MapleV returns 3 on
1189 <literal>A:=x^2+3; coeff(A,x,0);</literal> (GiNaC:
1190 <literal>A=pow(x,2)+3; coeff(A,x,0);</literal>) it is clear that
1191 MapleV is not trying to be consistent here. Also, users of MuPAD will
1192 in most cases feel more comfortable with GiNaC's convention. All
1193 function wrappers are always implemented as simple inline functions
1194 which just call the corresponding method and are only provided for
1195 users uncomfortable with OO who are dead set to avoid method
1196 invocations. Generally, a chain of function wrappers is much harder
1197 to read than a chain of methods and should therefore be avoided if
1198 possible. On the other hand, not everything in GiNaC is a method on
1199 class <literal>ex</literal> and sometimes calling a function cannot be
1203 <sect1><title>Polynomial Expansion</title>
1205 <para>A polynomial in one or more variables has many equivalent
1206 representations. Some useful ones serve a specific purpose. Consider
1207 for example the trivariate polynomial <literal>4*x*y + x*z + 20*y^2 +
1208 21*y*z + 4*z^2</literal>. It is equivalent to the factorized
1209 polynomial <literal>(x + 5*y + 4*z)*(4*y + z)</literal>. Other
1210 representations are the recursive ones where one collects for
1211 exponents in one of the three variable. Since the factors are
1212 themselves polynomials in the remaining two variables the procedure
1213 can be repeated. In our expample, two possibilies would be
1214 <literal>(4*y + z)*x + 20*y^2 + 21*y*z + 4*z^2</literal> and
1215 <literal>20*y^2 + (21*z + 4*x)*y + 4*z^2 + x*z</literal>.
1218 <para>To bring an expression into expanded form, its method
1219 <function>.expand()</function> may be called. In our example above,
1220 this corresponds to <literal>4*x*y + x*z + 20*y^2 + 21*y*z +
1221 4*z^2</literal>. Again, since the canonical form in GiNaC is not
1222 easily guessable you should be prepared to see different orderings of
1223 terms in such sums!</para>
1227 <sect1><title>Collecting expressions</title>
1229 <para>Another useful representation of multivariate polynomials is as
1230 a univariate polynomial in one of the variables with the coefficients
1231 being polynomials in the remaining variables. The method
1232 <literal>collect()</literal> accomplishes this task:
1234 <funcsynopsisinfo>#include <ginac/ginac.h></funcsynopsisinfo>
1235 <funcdef>ex <function>ex::collect</function></funcdef>
1236 <paramdef>symbol const & <parameter>s</parameter></paramdef>
1238 Note that the original polynomial needs to be in expanded form in
1239 order to be able to find the coefficients properly. The range of
1240 occuring coefficients can be checked using the two methods
1242 <funcsynopsisinfo>#include <ginac/ginac.h></funcsynopsisinfo>
1243 <funcdef>int <function>ex::degree</function></funcdef>
1244 <paramdef>symbol const & <parameter>s</parameter></paramdef>
1247 <funcdef>int <function>ex::ldegree</function></funcdef>
1248 <paramdef>symbol const & <parameter>s</parameter></paramdef>
1250 where <literal>degree()</literal> returns the highest coefficient and
1251 <literal>ldegree()</literal> the lowest one. These two methods work
1252 also reliably on non-expanded input polynomials. This is illustrated
1253 in the following example:
1255 <example><title>Collecting expressions in multivariate polynomials</title>
1257 #include <ginac/ginac.h>
1258 using namespace GiNaC;
1262 symbol x("x"), y("y");
1263 ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
1264 - pow(x+y,2) + 2*pow(y+2,2) - 8;
1265 ex Poly = PolyInp.expand();
1267 for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) {
1268 cout << "The x^" << i << "-coefficient is "
1269 << Poly.coeff(x,i) << endl;
1271 cout << "As polynomial in y: "
1272 << Poly.collect(y) << endl;
1277 When run, it returns an output in the following fashion:
1279 The x^0-coefficient is y^2+11*y
1280 The x^1-coefficient is 5*y^2-2*y
1281 The x^2-coefficient is -1
1282 The x^3-coefficient is 4*y
1283 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
1285 As always, the exact output may vary between different versions of
1286 GiNaC or even from run to run since the internal canonical ordering is
1287 not within the user's sphere of influence.</para>
1291 <sect1><title>Polynomial Arithmetic</title>
1293 <sect2><title>GCD and LCM</title>
1295 <para>The functions for polynomial greatest common divisor and least common
1296 multiple have the synopsis:
1298 <funcsynopsisinfo>#include <GiNaC/normal.h></funcsynopsisinfo>
1299 <funcdef>ex <function>gcd</function></funcdef>
1300 <paramdef>const ex *<parameter>a</parameter>, const ex *<parameter>b</parameter></paramdef>
1303 <funcdef>ex <function>lcm</function></funcdef>
1304 <paramdef>const ex *<parameter>a</parameter>, const ex *<parameter>b</parameter></paramdef>
1305 </funcsynopsis></para>
1307 <para>The functions <function>gcd()</function> and <function
1308 id="lcm-main">lcm()</function> accepts two expressions
1309 <literal>a</literal> and <literal>b</literal> as arguments and return
1310 a new expression, their greatest common divisor or least common
1311 multiple, respectively. If the polynomials <literal>a</literal> and
1312 <literal>b</literal> are coprime <function>gcd(a,b)</function> returns 1
1313 and <function>lcm(a,b)</function> returns the product of
1314 <literal>a</literal> and <literal>b</literal>.
1315 <example><title>Polynomal GCD/LCM</title>
1317 #include <ginac/ginac.h>
1318 using namespace GiNaC;
1322 symbol x("x"), y("y"), z("z");
1323 ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
1324 ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
1326 ex P_gcd = gcd(P_a, P_b);
1328 ex P_lcm = lcm(P_a, P_b);
1329 // 4*x*y^2 + 13*y*x*z + 20*y^3 + 81*y^2*z + 67*y*z^2 + 3*x*z^2 + 12*z^3
1338 <sect2><title>The <function>normal</function> method</title>
1340 <para>While in common symbolic code <function>gcd()</function> and
1341 <function>lcm()</function> are not too heavily used, simplification
1342 occurs frequently. Therefore <function>.normal()</function>, which
1343 provides some basic form of simplification, has become a method of
1344 class <literal>ex</literal>, just like <literal>.expand()</literal>.
1345 It converts a rational function into an equivalent rational function
1346 where numererator and denominator are coprime. This means, it finds
1347 the GCD of numerator and denominator and cancels it. If it encounters
1348 some object which does not belong to the domain of rationals (a
1349 function for instance), that object is replaced by a temporary symbol.
1350 This means that both expressions <literal>t1</literal> and
1351 <literal>t2</literal> are indeed simplified in this little program:
1352 <example><title>Cancellation of polynomial GCD (with obstacles)</title>
1354 #include <ginac/ginac.h>
1355 using namespace GiNaC;
1360 ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
1361 ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
1362 cout << "t1 is " << t1.normal() << endl;
1363 cout << "t2 is " << t2.normal() << endl;
1369 Of course this works for multivariate polynomials too, so the ratio of
1370 the sample-polynomials from the section about GCD and LCM above would
1371 be normalized to <literal>P_a/P_b</literal> =
1372 <literal>(4*y+z)/(y+3*z)</literal>.</para>
1378 <sect1><title>Symbolic Differentiation</title>
1380 <para>GiNaC's objects know how to differentiate themselves. Thus, a
1381 polynomial (class <literal>add</literal>) knows that its derivative is
1382 the sum of the derivatives of all the monomials:
1383 <example><title>Simple polynomial differentiation</title>
1385 #include <ginac/ginac.h>
1386 using namespace GiNaC;
1390 symbol x("x"), y("y"), z("z");
1391 ex P = pow(x, 5) + pow(x, 2) + y;
1393 cout << P.diff(x,2) << endl; // 20*x^3 + 2
1394 cout << P.diff(y) << endl; // 1
1395 cout << P.diff(z) << endl; // 0
1400 If a second integer parameter <emphasis>n</emphasis> is given the
1401 <function>diff</function> method returns the <emphasis>n</emphasis>th
1404 <para>If <emphasis>every</emphasis> object and every function is told
1405 what its derivative is, all derivatives of composed objects can be
1406 calculated using the chain rule and the product rule. Consider, for
1407 instance the expression <literal>1/cosh(x)</literal>. Since the
1408 derivative of <literal>cosh(x)</literal> is <literal>sinh(x)</literal>
1409 and the derivative of <literal>pow(x,-1)</literal> is
1410 <literal>-pow(x,-2)</literal>, GiNaC can readily compute the
1411 composition. It turns out that the composition is the generating
1412 function for Euler Numbers, i.e. the so called
1413 <emphasis>n</emphasis>th Euler number is the coefficient of
1414 <literal>x^n/n!</literal> in the expansion of
1415 <literal>1/cosh(x)</literal>. We may use this identity to code a
1416 function that generates Euler numbers in just three lines:
1417 <example><title>Differentiation with nontrivial functions: Euler numbers</title>
1419 #include <ginac/ginac.h>
1420 using namespace GiNaC;
1422 ex EulerNumber(unsigned n)
1425 ex generator = pow(cosh(x),-1);
1426 return generator.diff(x,n).subs(x==0);
1431 for (unsigned i=0; i<11; i+=2)
1432 cout << EulerNumber(i) << endl;
1437 When you run it, it produces the sequence <literal>1</literal>,
1438 <literal>-1</literal>, <literal>5</literal>, <literal>-61</literal>,
1439 <literal>1385</literal>, <literal>-50521</literal>. We increment the
1440 loop variable <literal>i</literal> by two since all odd Euler numbers
1441 vanish anyways.</para>
1445 <sect1><title>Series Expansion</title>
1447 <para>Expressions know how to expand themselves as a Taylor series or
1448 (more generally) a Laurent series. As in most conventional Computer
1449 Algebra Systems no distinction is made between those two. There is a
1450 class of its own for storing such series as well as a class for
1451 storing the order of the series. A sample program could read:
1452 <example><title>Series expansion</title>
1454 #include <ginac/ginac.h>
1455 using namespace GiNaC;
1461 ex MyExpr1 = sin(x);
1462 ex MyExpr2 = 1/(x - pow(x, 2) - pow(x, 3));
1463 ex MyTailor, MySeries;
1465 MyTailor = MyExpr1.series(x, point, 5);
1466 cout << MyExpr1 << " == " << MyTailor
1467 << " for small " << x << endl;
1468 MySeries = MyExpr2.series(x, point, 7);
1469 cout << MyExpr2 << " == " << MySeries
1470 << " for small " << x << endl;
1477 <para>As an instructive application, let us calculate the numerical
1478 value of Archimedes' constant (for which there already exists the
1479 built-in constant <literal>Pi</literal>) using Méchain's
1480 mysterious formula <literal>Pi==16*atan(1/5)-4*atan(1/239)</literal>.
1481 We may expand the arcus tangent around <literal>0</literal> and insert
1482 the fractions <literal>1/5</literal> and <literal>1/239</literal>.
1483 But, as we have seen, a series in GiNaC carries an order term with it.
1484 The function <literal>series_to_poly</literal> may be used to strip
1486 <example><title>Series expansion using Méchain's formula for
1487 <literal>Pi</literal></title>
1489 #include <ginac/ginac.h>
1490 using namespace GiNaC;
1492 ex mechain_pi(int degr)
1495 ex pi_expansion = series_to_poly(atan(x).series(x,0,degr));
1496 ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
1497 -4*pi_expansion.subs(x==numeric(1,239));
1504 for (int i=2; i<12; i+=2) {
1505 pi_frac = mechain_pi(i);
1506 cout << i << ":\t" << pi_frac << endl
1507 << "\t" << pi_frac.evalf() << endl;
1512 <para>When you run this program, it will type out:</para>
1515 3.1832635983263598326
1516 4: 5359397032/1706489875
1517 3.1405970293260603143
1518 6: 38279241713339684/12184551018734375
1519 3.141621029325034425
1520 8: 76528487109180192540976/24359780855939418203125
1521 3.141591772182177295
1522 10: 327853873402258685803048818236/104359128170408663038552734375
1523 3.1415926824043995174
1534 <title>Extending GiNaC</title>
1536 <para>By reading so far you should have gotten a fairly good
1537 understanding of GiNaC's design-patterns. From here on you should
1538 start reading the sources. All we can do now is issue some
1539 recommendations how to tackle GiNaC's many loose ends in order to
1540 fulfill everybody's dreams.</para>
1542 <sect1><title>What doesn't belong into GiNaC</title>
1544 <para>First of all, GiNaC's name must be read literally. It is
1545 designed to be a library for use within C++. The tiny
1546 <literal>ginsh</literal> accompanying GiNaC makes this even more
1547 clear: it doesn't even attempt to provide a language. There are no
1548 loops or conditional expressions in <literal>ginsh</literal>, it is
1549 merely a window into the library for the programmer to test stuff (or
1550 to show off). Still, the design of a complete CAS with a language of
1551 its own, graphical capabilites and all this on top of GiNaC is
1552 possible and is without doubt a nice project for the future.</para>
1554 <para>There are many built-in functions in GiNaC that do not know how
1555 to evaluate themselves numerically to a precision declared at runtime
1556 (using <literal>Digits</literal>). Some may be evaluated at certain
1557 points, but not generally. This ought to be fixed. However, doing
1558 numerical computations with GiNaC's quite abstract classes is doomed
1559 to be inefficient. For this purpose, the underlying bignum-package
1560 <literal>CLN</literal> is much better suited.</para>
1564 <sect1><title>Other symbolic functions</title>
1566 <para>The easiest and most instructive way to start with is probably
1567 to implement your own function. Objects of class
1568 <literal>function</literal> are inserted into the system via a kind of
1569 "registry". They get a serial number that is used internally to
1570 identify them but you usually need not worry about this. What you
1571 have to care for are functions that are called when the user invokes
1572 certain methods. These are usual C++-functions accepting a number of
1573 <literal>ex</literal> as arguments and returning one
1574 <literal>ex</literal>. As an example, if we have a look at a
1575 simplified implementation of the cosine trigonometric function, we
1576 first need a function that is called when one wishes to
1577 <literal>eval</literal> it. It could look something like this:
1580 static ex cos_eval_method(ex const & x)
1582 // if x%2*Pi return 1
1583 // if x%Pi return -1
1584 // if x%Pi/2 return 0
1585 // care for other cases...
1586 return cos(x).hold();
1589 The last line returns <literal>cos(x)</literal> if we don't know what
1590 else to do and stops a potential recursive evaluation by saying
1591 <literal>.hold()</literal>. We should also implement a method for
1592 numerical evaluation and since we are lazy we sweep the problem under
1593 the rug by calling someone else's function that does so, in this case
1594 the one in class <literal>numeric</literal>:
1596 static ex cos_evalf_method(ex const & x)
1598 return sin(ex_to_numeric(x));
1601 Differentiation will surely turn up and so we need to tell
1602 <literal>sin</literal> how to differentiate itself:
1604 static ex cos_diff_method(ex const & x, unsigned diff_param)
1610 The second parameter is obligatory but uninteresting at this point.
1611 It is used for correct handling of the product rule only. For Taylor
1612 expansion, it is enough to know how to differentiate. But if the
1613 function you want to implement does have a pole somewhere in the
1614 complex plane, you need to write another method for Laurent expansion
1615 around that point.</para>
1617 <para>Now that everything has been written for <literal>cos</literal>,
1618 we need to tell the system about it. This is done by a macro and we
1619 are not going to descibe how it expands, please consult your
1620 preprocessor if you are curious:
1622 REGISTER_FUNCTION(cos, cos_eval_method, cos_evalf_method, cos_diff, NULL);
1624 The first argument is the function's name, the second, third and
1625 fourth bind the corresponding methods to this objects and the fifth is
1626 a slot for inserting a method for series expansion. Also, the new
1627 function needs to be declared somewhere. This may also be done by a
1628 convenient preprocessor macro:
1630 DECLARE_FUNCTION_1P(cos)
1632 The suffix <literal>_1P</literal> stands for <emphasis>one
1633 parameter</emphasis>. Of course, this implementation of
1634 <literal>cos</literal> is very incomplete and lacks several safety
1635 mechanisms. Please, have a look at the real implementation in GiNaC.
1636 (By the way: in case you are worrying about all the macros above we
1637 can assure you that functions are GiNaC's most macro-intense classes.
1638 We have done our best to avoid them where we can.)</para>
1640 <para>That's it. May the source be with you!</para>
1648 <title>A Comparison with other CAS</title>
1650 <para>This chapter will give you some information on how GiNaC
1651 compares to other, traditional Computer Algebra Systems, like
1652 <literal>Maple</literal>, <literal>Mathematica</literal> or
1653 <literal>Reduce</literal>, where it has advantages and disadvantages
1654 over these systems.</para>
1656 <sect1><title>Advantages</title>
1658 <para>GiNaC has several advantages over traditional Computer
1659 Algebra Systems, like
1663 <para>familiar language: all common CAS implement their own
1664 proprietary grammar which you have to learn first (and maybe learn
1665 again when your vendor chooses to "enhance" it). With GiNaC you
1666 can write your program in common <literal>C++</literal>, which is
1667 standardized.</para>
1670 <para>structured data types: you can build up structured data
1671 types using <literal>struct</literal>s or <literal>class</literal>es
1672 together with STL features instead of using unnamed lists of lists
1676 <para>strongly typed: in CAS, you usually have only one kind of
1677 variables which can hold contents of an arbitrary type. This
1678 4GL like feature is nice for novice programmers, but dangerous.
1682 <para>development tools: powerful development tools exist for
1683 <literal>C++</literal>, like fancy editors (e.g. with automatic
1684 indentation and syntax highlighting), debuggers, visualization
1685 tools, documentation tools...</para>
1688 <para>modularization: <literal>C++</literal> programs can
1689 easily be split into modules by separating interface and
1690 implementation.</para>
1693 <para>price: GiNaC is distributed under the GNU Public License
1694 which means that it is free and available with source code. And
1695 there are excellent <literal>C++</literal>-compilers for free, too.
1699 <para>extendable: you can add your own classes to GiNaC, thus
1700 extending it on a very low level. Compare this to a traditional
1701 CAS that you can usually only extend on a high level by writing in
1702 the language defined by the parser. In particular, it turns out
1703 to be almost impossible to fix bugs in a traditional system.
1706 <para>seemless integration: it is somewhere between difficult
1707 and impossible to call CAS functions from within a program
1708 written in <literal>C++</literal> or any other programming
1709 language and vice versa. With GiNaC, your symbolic routines
1710 are part of your program. You can easily call third party
1711 libraries, e.g. for numerical evaluation or graphical
1712 interaction. All other approaches are much more cumbersome: they
1713 range from simply ignoring the problem
1714 (i.e. <literal>Maple</literal>) to providing a
1715 method for "embedding" the system
1716 (i.e. <literal>Yacas</literal>).</para>
1719 <para>efficiency: often large parts of a program do not need
1720 symbolic calculations at all. Why use large integers for loop
1721 variables or arbitrary precision arithmetics where double
1722 accuracy is sufficient? For pure symbolic applications,
1723 GiNaC is comparable in speed with other CAS.
1728 <sect1><title>Disadvantages</title>
1730 <para>Of course it also has some disadvantages
1734 <para>not interactive: GiNaC programs have to be written in
1735 an editor, compiled and executed. You cannot play with
1736 expressions interactively. However, such an extension is not
1737 inherently forbidden by design. In fact, two interactive
1738 interfaces are possible: First, a simple shell that exposes GiNaC's
1739 types to a command line can readily be written (and has been
1740 written) and second, as a more consistent approach we plan
1741 an integration with the <literal>CINT</literal>
1742 <literal>C++</literal> interpreter.</para>
1745 <para>advanced features: GiNaC cannot compete with a program
1746 like <literal>Reduce</literal> which exists for more than
1747 30 years now or <literal>Maple</literal> which grows since
1748 1981 by the work of dozens of programmers, with respect to
1749 mathematical features. Integration, factorization, non-trivial
1750 simplifications, limits etc. are missing in GiNaC (and are not
1751 planned for the near future).</para>
1754 <para>portability: While the GiNaC library itself is designed
1755 to avoid any platform dependent features (it should compile
1756 on any ANSI compliant <literal>C++</literal> compiler), the
1757 currently used version of the CLN library (fast large integer and
1758 arbitrary precision arithmetics) can be compiled only on systems
1759 with a recently new <literal>C++</literal> compiler from the
1760 GNU Compiler Collection (<literal>GCC</literal>). GiNaC uses
1761 recent language features like explicit constructors, mutable
1762 members, RTTI, dynamic_casts and STL, so ANSI compliance is meant
1763 literally. Recent <literal>GCC</literal> versions starting at
1764 2.95, although itself not yet ANSI compliant, support all needed
1771 <sect1><title>Why <literal>C++</literal>?</title>
1773 <para>Why did we choose to implement GiNaC in <literal>C++</literal>
1774 instead of <literal>Java</literal> or any other language?
1775 <literal>C++</literal> is not perfect: type checking is not strict
1776 (casting is possible), separation between interface and implementation
1777 is not complete, object oriented design is not enforced. The main
1778 reason is the often scolded feature of operator overloading in
1779 <literal>C++</literal>. While it may be true that operating on classes
1780 with a <literal>+</literal> operator is rarely meaningful, it is
1781 perfectly suited for algebraic expressions. Writing 3x+5y as
1782 <literal>3*x+5*y</literal> instead of
1783 <literal>x.times(3).plus(y.times(5))</literal> looks much more
1784 natural. Furthermore, the main developers are more familiar with
1785 <literal>C++</literal> than with any other programming
1796 <title>ISO/IEC 14882:1998</title>
1797 <subtitle>Programming Languages: C++</subtitle>
1802 <title>CLN: A Class Library for Numbers</title>
1805 <firstname>Bruno</firstname><surname>Haible</surname>
1806 <affiliation><address><email>haible@ilog.fr</email></address></affiliation>
1813 <title>The C++ Programming Language</title>
1814 <authorgroup><author><firstname>Bjarne</firstname><surname>Stroustrup</surname></author></authorgroup>
1815 <edition>3</edition>
1816 <isbn>0-201-88954-4</isbn>
1817 <publisher><publishername>Addison Wesley</publishername></publisher>
1823 <title>C++ FAQs</title>
1824 <authorgroup><author><firstname>Marshall</firstname><surname>Cline</surname></author></authorgroup>
1825 <isbn>0-201-58958-3</isbn>
1826 <pubdate>1995</pubdate>
1827 <publisher><publishername>Addison Wesley</publishername></publisher>
1833 <title>Algorithms for Computer Algebra</title>
1835 <author><firstname>Keith</firstname><othername>O.</othername><surname>Geddes</surname></author>
1836 <author><firstname>Stephen</firstname><othername>R.</othername><surname>Czapor</surname></author>
1837 <author><firstname>George</firstname><surname>Labahn</surname></author>
1839 <isbn>0-7923-9259-0</isbn>
1840 <pubdate>1992</pubdate>
1842 <publishername>Kluwer Academic Publishers</publishername>
1843 <address><city>Norwell</city>, <state>Massachusetts</state></address>
1850 <title>Computer Algebra</title>
1851 <subtitle>Systems and Algorithms for Algebraic Computation</subtitle>
1853 <author><firstname>J.</firstname><othername>H.</othername><surname>Davenport</surname></author>
1854 <author><firstname>Y.</firstname><surname>Siret</surname></author>
1855 <author><firstname>E.</firstname><surname>Tournier</surname></author>
1857 <isbn>0-12-204230-1</isbn>
1858 <pubdate>1988</pubdate>
1860 <publishername>Academic Press</publishername>
1861 <address><city>London</city></address>