1 <!DOCTYPE Book PUBLIC "-//Davenport//DTD DocBook V3.0//EN">
4 <title>GiNaC 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 GiNaC
48 but already has some background in C++ programming. However, since a
49 hand made documentation like this one is difficult to keep in sync
50 with the development the actual documentation is inside the sources in
51 the form of comments. That documentation may be parsed by one of the
52 many Javadoc-like documentation systems. The generated HTML
53 documenatation is included in the distributed sources (subdir
54 <literal>doc/reference/</literal>) or can be accessed directly at URL
56 url="http://wwwthep.physik.uni-mainz.de/GiNaC/reference/"><literal>http://wwwthep.physik.uni-mainz.de/GiNaC/reference/</literal></ulink>.
57 It is an invaluable resource not only for the advanced user who wishes
58 to extend the system (or chase bugs) but for everybody who wants to
59 comprehend the inner workings of GiNaC. This little tutorial on the
60 other hand only covers the basic things that are unlikely to change in
64 <sect1><title>License</title>
66 <para>The GiNaC framework for symbolic computation within the C++
67 programming language is Copyright (C) 1999 Johannes Gutenberg
68 Universität Mainz, Germany.</para>
70 <para>This program is free software; you can redistribute it and/or
71 modify it under the terms of the GNU General Public License as
72 published by the Free Software Foundation; either version 2 of the
73 License, or (at your option) any later version.</para>
75 <para>This program is distributed in the hope that it will be useful, but
76 WITHOUT ANY WARRANTY; without even the implied warranty of
77 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
78 General Public License for more details.</para>
80 <para>You should have received a copy of the GNU General Public License
81 along with this program; see the file COPYING. If not, write to the
82 Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
83 MA 02111-1307, USA.</para>
88 <title>A Tour of GiNaC</title>
90 <para>This quick tour of GiNaC wants to rise your interest in in the
91 subsequent chapters by showing off a bit. Please excuse us if it
92 leaves many open questions.</para>
94 <sect1><title>How to use it from within C++</title> <para>The GiNaC
95 open framework for symbolic computation within the C++ programming
96 language does not try to define a language of it's own as conventional
97 CAS do. Instead, it extends the capabilities of C++ by symbolic
98 manipulations. Here is how to generate and print a simple (and
99 pointless) bivariate polynomial with some large coefficients:
101 <title>My first GiNaC program (a bivariate polynomial)</title>
103 #include <ginac/ginac.h>
104 using namespace GiNaC;
108 symbol x("x"), y("y");
111 for (int i=0; i<3; ++i)
112 poly += factorial(i+16)*pow(x,i)*pow(y,2-i);
114 cout << poly << endl;
118 <para>Assuming the file is called <literal>hello.cc</literal>, on
119 our system we can compile and run it like this:</para>
121 <prompt>$</prompt> c++ hello.cc -o hello -lcln -lginac
122 <prompt>$</prompt> ./hello
123 355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
128 <para>Next, there is a more meaningful C++ program that calls a
129 function which generates Hermite polynomials in a specified free
132 <title>My second GiNaC program (Hermite polynomials)</title>
134 #include <ginac/ginac.h>
135 using namespace GiNaC;
137 ex HermitePoly(symbol x, int deg)
139 ex HKer=exp(-pow(x,2));
140 // uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2)
141 return normal(pow(-1,deg) * diff(HKer, x, deg) / HKer);
148 for (int i=0; i<6; ++i)
149 cout << "H_" << i << "(z) == " << HermitePoly(z,i) << endl;
154 <para>When run, this will type out</para>
159 H_3(z) == -12*z+8*z^3
160 H_4(z) == -48*z^2+16*z^4+12
161 H_5(z) == 120*z-160*z^3+32*z^5
164 This method of generating the coefficients is of course far from
165 optimal for production purposes.</para>
167 <para>In order to show some more examples of what GiNaC can do we
168 will now use <literal>ginsh</literal>, a simple GiNaC interactive
169 shell that provides a convenient window into GiNaC's capabilities.
172 <sect1><title>What it can do for you</title>
174 <para>After invoking <literal>ginsh</literal> one can test and
175 experiment with GiNaC's features much like in other Computer Algebra
176 Systems except that it does not provide programming constructs like
177 loops or conditionals. For a concise description of the
178 <literal>ginsh</literal> syntax we refer to its accompanied man-page.
179 Suffice to say that assignments and comparisons in
180 <literal>ginsh</literal> are written as they are in C,
181 i.e. <literal>=</literal> assigns and <literal>==</literal>
184 <para>It can manipulate arbitrary precision integers in a very fast
185 way. Rational numbers are automatically converted to fractions of
189 369988485035126972924700782451696644186473100389722973815184405301748249
191 123329495011708990974900260817232214728824366796574324605061468433916083
199 <para>All numbers occuring in GiNaC's expressions can be converted
200 into floating point numbers with the <literal>evalf</literal> method,
201 to arbitrary accuracy:
204 0.14285714285714285714
208 0.1428571428571428571428571428571428571428571428571428571428571428571428
209 5714285714285714285714285714285714285
213 <para>Exact numbers other than rationals that can be manipulated in
214 GiNaC include predefined constants like Archimedes' Pi. They can both
215 be used in symbolic manipulations (as an exact number) as well as in
216 numeric expressions (as an inexact number):
221 x+9.869604401089358619L0
225 11.869604401089358619L0
229 <para>Built-in functions evaluate immediately to exact numbers if
230 this is possible. Conversions that can be safely performed are done
231 immediately; conversions that are not generally valid are not done:
240 (Note that converting the last input to <literal>x</literal> would
241 allow one to conclude that <literal>42*Pi</literal> is equal to
242 <literal>0</literal>.)</para>
244 <para>Linear equation systems can be solved along with basic linear
245 algebra manipulations over symbolic expressions. In C++ there is a
246 matrix class for this purpose but we can see what it can do using
247 <literal>ginsh</literal>'s notation of double brackets to type them in:
249 > lsolve(a+x*y==z,x);
251 lsolve([3*x+5*y == 7, -2*x+10*y == -5], [x, y]);
253 > M = [[ [[1, 3]], [[-3, 2]] ]];
254 [[ [[1,3]], [[-3,2]] ]]
257 > charpoly(M,lambda);
262 <para>Multivariate polynomials and rational functions may be expanded,
263 collected and normalized (i.e. converted to a ratio of two coprime
266 > a = x^4 + 2*x^2*y^2 + 4*x^3*y + 12*x*y^3 - 3*y^4;
267 -3*y^4+x^4+12*x*y^3+2*x^2*y^2+4*x^3*y
268 > b = x^2 + 4*x*y - y^2;
271 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
273 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)
280 You can differentiate functions and expand them as Taylor or Laurent
281 series (the third argument of series is the evaluation point, the
282 fourth defines the order):
286 > series(sin(x),x,0,4);
288 > series(1/tan(x),x,0,4);
289 x^(-1)-1/3*x+Order(x^2)
299 <title>Installation</title>
301 <para>GiNaC's installation follows the spirit of most GNU software. It is
302 easily installed on your system by three steps: configuration, build,
305 <sect1><title id="CLN-main">Prerequistes</title>
307 <para>In order to install GiNaC on your system, some prerequistes need
308 to be met. First of all, you need to have a C++-compiler adhering to
309 the ANSI-standard <citation>ISO/IEC 14882:1998(E)</citation>. We used
310 <literal>GCC</literal> for development so if you have a different
311 compiler you are on your own. For the configuration to succeed you
312 need a Posix compliant shell installed in <literal>/bin/sh</literal>,
313 GNU <literal>bash</literal> is fine. Perl is needed by the built
314 process as well, since some of the source files are automatically
315 generated by Perl scripts. Last but not least, Bruno Haible's library
316 <literal>CLN</literal> is extensively used and needs to be installed
317 on your system. Please get it from <ulink
318 url="ftp://ftp.santafe.edu/pub/gnu/"><literal>ftp://ftp.santafe.edu/pub/gnu/</literal></ulink>
320 url="ftp://ftp.ilog.fr/pub/Users/haible/gnu/"><literal>ftp://ftp.ilog.fr/pub/Users/haible/gnu/</literal></ulink>
321 (it is covered by GPL) and install it prior to trying to install
322 GiNaC. The configure script checks if it can find it and if it cannot
323 it will refuse to continue.</para></sect1>
325 <sect1><title>Configuration</title>
327 <para>To configure GiNaC means to prepare the source distribution for
328 building. It is done via a shell script called
329 <literal>configure</literal> that is shipped with the sources.
330 (Actually, this script is by itself created with GNU Autoconf from the
331 files <literal>configure.in</literal> and
332 <literal>aclocal.m4</literal>.) Since a configure script generated by
333 GNU Autoconf never prompts, all customization must be done either via
334 command line parameters or environment variables. It accepts a list
335 of parameters, the complete set of which can be listed by calling it
336 with the <literal>--help</literal> option. The most important ones
337 will be shortly described in what follows:
340 <para><literal>--disable-shared</literal>: When given, this option
341 switches off the build of a shared library, i.e. a
342 <literal>.so</literal>-file. This may be convenient when developing
343 because it considerably speeds up compilation.</para>
346 <para><literal>--prefix=</literal><emphasis>PREFIX</emphasis>: The
347 directory where the compiled library and headers are installed. It
348 defaults to <literal>/usr/local</literal> which means that the
349 library is installed in the directory
350 <literal>/usr/local/lib</literal> and the header files in
351 <literal>/usr/local/include/GiNaC</literal> and the documentation
352 (like this one) into <literal>/usr/local/share/doc/GiNaC</literal>.</para>
355 <para><literal>--libdir=</literal><emphasis>LIBDIR</emphasis>: Use
356 this option in case you want to have the library installed in some
358 <emphasis>PREFIX</emphasis><literal>/lib/</literal>.</para>
361 <para><literal>--includedir=</literal><emphasis>INCLUDEDIR</emphasis>:
362 Use this option in case you want to have the header files
363 installed in some other directory than
364 <emphasis>PREFIX</emphasis><literal>/include/ginac/</literal>. For
365 instance, if you specify
366 <literal>--includedir=/usr/include</literal> you will end up with
367 the header files sitting in the directory
368 <literal>/usr/include/ginac/</literal>. Note that the subdirectory
369 <literal>GiNaC</literal> is enforced by this process in order to
370 keep the header files separated from others. This avoids some
371 clashes and allows for an easier deinstallation of GiNaC. This ought
372 to be considered A Good Thing (tm).</para>
375 <para><literal>--datadir=</literal><emphasis>DATADIR</emphasis>:
376 This option may be given in case you want to have the documentation
377 installed in some other directory than
378 <emphasis>PREFIX</emphasis><literal>/share/doc/GiNaC/</literal>.
383 <para>In addition, you may specify some environment variables.
384 <literal>CXX</literal> holds the path and the name of the C++ compiler
385 in case you want to override the default in your path. (The
386 <literal>configure</literal> script searches your path for
387 <literal>c++</literal>, <literal>g++</literal>,
388 <literal>gcc</literal>, <literal>CC</literal>, <literal>cxx</literal>
389 and <literal>cc++</literal> in that order.) It may be very useful to
390 define some compiler flags with the <literal>CXXFLAGS</literal>
391 environment variable, like optimization, debugging information and
392 warning levels. If ommitted, it defaults to <literal>-g
393 -O2</literal>.</para>
395 <para>The whole process is illustrated in the following two
396 examples. (Substitute <literal>setenv VARIABLE value</literal> for
397 <literal>export VARIABLE=value</literal> if the Berkeley C shell is
400 <example><title>Sample sessions of how to call the
401 configure-script</title> <para>Simple configuration for a site-wide
402 GiNaC library assuming everything is in default paths:</para>
404 <prompt>$</prompt> export CXXFLAGS="-Wall -O2"
405 <prompt>$</prompt> ./configure
407 <para>Configuration for a private static GiNaC library with several
408 components sitting in custom places (site-wide <literal>GCC</literal>
409 and private <literal>CLN</literal>), the compiler pursueded to be
410 picky and full assertions switched on:</para>
412 <prompt>$</prompt> export CXX=/usr/local/gnu/bin/c++
413 <prompt>$</prompt> export CPPFLAGS="${CPPFLAGS} -I${HOME}/include"
414 <prompt>$</prompt> export CXXFLAGS="${CXXFLAGS} -ggdb -Wall -ansi -pedantic -O2 -DDO_GINAC_ASSERT"
415 <prompt>$</prompt> export LDFLAGS="${LDFLAGS} -L${HOME}/lib"
416 <prompt>$</prompt> ./configure --disable-shared --prefix=${HOME}
423 <sect1><title>Building GiNaC</title>
425 <para>After proper configuration you should just build the whole
426 library by typing <literal>make</literal> at the command
427 prompt and go for a cup of coffee.</para>
429 <para>Just to make sure GiNaC works properly you may run a simple test
430 suite by typing <literal>make check</literal>. This will compile some
431 sample programs, run them and compare the output to reference output.
432 Each of the checks should return a message <literal>passed</literal>
433 together with the CPU time used for that particular test. If it does
434 not, something went wrong. This is mostly intended to be a QA-check
435 if something was broken during the development, but not a sanity check
436 of your system. Another intent is to allow people to fiddle around
437 with optimization. If <literal>CLN</literal> was installed all right
438 this step is unlikely to return any errors.</para>
442 <sect1><title>Installation</title>
444 <para>To install GiNaC on your system, simply type <literal>make
445 install</literal>. As described in the section about configuration
446 the files will be installed in the following directories (the
447 directories will be created if they don't already exist):
450 <para><literal>libginac.a</literal> will go into
451 <emphasis>PREFIX</emphasis><literal>/lib/</literal> (or
452 <emphasis>LIBDIR</emphasis>) which defaults to
453 <literal>/usr/local/lib/</literal>. So will
454 <literal>libginac.so</literal> if the the configure script was
455 given the option <literal>--enable-shared</literal>. In that
456 case, the proper symlinks will be established as well (by running
457 <literal>ldconfig</literal>).</para>
460 <para>All the header files will be installed into
461 <emphasis>PREFIX</emphasis><literal>/include/ginac/</literal> (or
462 <emphasis>INCLUDEDIR</emphasis><literal>/ginac/</literal>, if
466 <para>All documentation (HTML, Postscript and DVI) will be stuffed
468 <emphasis>PREFIX</emphasis><literal>/share/doc/GiNaC/</literal>
469 (or <emphasis>DATADIR</emphasis><literal>/doc/GiNaC/</literal>, if
475 <para>Just for the record we will list some other useful make targets:
476 <literal>make clean</literal> deletes all files generated by
477 <literal>make</literal>, i.e. all the object files. In addition
478 <literal>make distclean</literal> removes all files generated by
479 configuration. And finally <literal>make uninstall</literal> removes
480 the installed library and header files<footnoteref
481 linkend="warnuninstall-1">.
483 <footnote id="warnuninstall-1"><para>Uninstallation does not work
484 after you have called <literal>make distclean</literal> since the
485 <literal>Makefile</literal> is itself generated by the configuration
486 from <literal>Makefile.in</literal> and hence deleted by
487 <literal>make distclean</literal>. There are two obvious ways out
488 of this dilemma. First, you can run the configuration again with
489 the same <emphasis>PREFIX</emphasis> thus creating a
490 <literal>Makefile</literal> with a working
491 <literal>uninstall</literal> target. Second, you can do it by hand
492 since you now know where all the files went during
493 installation.</para></footnote>
501 <title>Basic Concepts</title>
503 <para>This chapter will describe the different fundamental objects
504 that can be handled with GiNaC. But before doing so, it is worthwhile
505 introducing you to the more commonly used class of expressions,
506 representing a flexible meta-class for storing all mathematical
509 <sect1><title>Expressions</title>
511 <para>The most common class of objects a user deals with is the
512 expression <literal>ex</literal>, representing a mathematical object
513 like a variable, number, function, sum, product, etc... Expressions
514 may be put together to form new expressions, passed as arguments to
515 functions, and so on. Here is a little collection of valid
517 <example><title>Examples of expressions</title>
519 ex MyEx1 = 5; // simple number
520 ex MyEx2 = x + 2*y; // polynomial in x and y
521 ex MyEx3 = (x + 1)/(x - 1); // rational expression
522 ex MyEx4 = sin(x + 2*y) + 3*z + 41; // containing a function
523 ex MyEx5 = MyEx4 + 1; // similar to above
526 Before describing the more fundamental objects that form the building
527 blocks of expressions we'll have a quick look under the hood by
528 describing how expressions are internally managed.</para>
530 <sect2><title>Digression: Expressions are reference counted</title>
532 <para>An expression is extremely light-weight since internally it
533 works like a handle to the actual representation and really holds
534 nothing more than a pointer to some other object. What this means in
535 practice is that whenever you create two <literal>ex</literal> and set
536 the second equal to the first no copying process is involved. Instead,
537 the copying takes place as soon as you try to change the second.
538 Consider the simple sequence of code:
539 <example><title>Simple copy-on-write semantics</title>
541 #include <ginac/ginac.h>
542 using namespace GiNaC;
546 symbol x("x"), y("y"), z("z");
549 e1 = sin(x + 2*y) + 3*z + 41;
550 e2 = e1; // e2 points to same object as e1
551 cout << e2 << endl; // prints sin(x+2*y)+3*z+41
552 e2 += 1; // e2 is copied into a new object
553 cout << e2 << endl; // prints sin(x+2*y)+3*z+42
558 The line <literal>e2 = e1;</literal> creates a second expression
559 pointing to the object held already by <literal>e1</literal>. The
560 time involved for this operation is therefore constant, no matter how
561 large <literal>e1</literal> was. Actual copying, however, must take
562 place in the line <literal>e2 += 1</literal> because
563 <literal>e1</literal> and <literal>e2</literal> are not handles for
564 the same object any more. This concept is called
565 <emphasis>copy-on-write semantics</emphasis>. It increases
566 performance considerably whenever one object occurs multiple times and
567 represents a simple garbage collection scheme because when an
568 <literal>ex</literal> runs out of scope its destructor checks whether
569 other expressions handle the object it points to too and deletes the
570 object from memory if that turns out not to be the case. A slightly
571 less trivial example of differentiation using the chain-rule should
572 make clear how powerful this can be. <example><title>Advanced
573 copy-on-write semantics</title>
575 #include <ginac/ginac.h>
576 using namespace GiNaC;
580 symbol x("x"), y("y");
584 ex e3 = diff(sin(e2), x); // first derivative of sin(e2) by x
585 cout << e1 << endl // prints x+3*y
586 << e2 << endl // prints (x+3*y)^3
587 << e3 << endl; // prints 3*(x+3*y)^2*cos((x+3*y)^3)
592 Here, <literal>e1</literal> will actually be referenced three times
593 while <literal>e2</literal> will be referenced two times. When the
594 power of an expression is built, that expression needs not be
595 copied. Likewise, since the derivative of a power of an expression can
596 be easily expressed in terms of that expression, no copying of
597 <literal>e1</literal> is involved when <literal>e3</literal> is
598 constructed. So, when <literal>e3</literal> is constructed it will
599 print as <literal>3*(x+3*y)^2*cos((x+3*y)^3)</literal> but the
600 argument of <literal>cos()</literal> only holds a reference to
601 <literal>e2</literal> and the factor in front is just
602 <literal>3*e1^2</literal>.
605 <para>As a user of GiNaC, you cannot see this mechanism of
606 copy-on-write semantics. When you insert an expression into a second
607 expression, the result behaves exactly as if the contents of the first
608 expression were inserted. But it may be useful to remember that this
609 is not what happens. Knowing this will enable you to write much more
610 efficient code. If you still have an uncertain feeling with
611 copy-on-write semantics, we recommend you have a look at the
612 <emphasis>C++-FAQ lite</emphasis> by Marshall Cline. Chapter 16
613 covers this issue and presents an implementation which is pretty close
614 to the one in GiNaC. You can find it on the Web at <ulink
615 url="http://www.cerfnet.com/~mpcline/c++-faq-lite/"><literal>http://www.cerfnet.com/~mpcline/c++-faq-lite/</literal></ulink>.</para>
617 <para>So much for expressions. But what exactly are these expressions
618 handles of? This will be answered in the following sections.</para>
622 <sect1><title>The Class Hierarchy</title>
624 <para>GiNaC's class hierarchy consists of several classes representing
625 mathematical objects, all of which (except for <literal>ex</literal>
626 and some helpers) are internally derived from one abstract base class
627 called <literal>basic</literal>. You do not have to deal with objects
628 of class <literal>basic</literal>, instead you'll be dealing with
629 symbols and functions of symbols. You'll soon learn in this chapter
630 how many of the functions on symbols are really classes. This is
631 because simple symbolic arithmetic is not supported by languages like
632 C++ so in a certain way GiNaC has to implement its own arithmetic.</para>
634 <para>To give an idea about what kinds of symbolic composits may be
635 built we have a look at the most important classes in the class
636 hierarchy. The dashed line symbolizes a "points to" or "handles"
637 relationship while the solid lines stand for "inherits from"
639 <figure id="classhier-id" float="1">
640 <title>The GiNaC class hierarchy</title>
641 <graphic align="center" fileref="classhierarchy.graext" format="GRAEXT"></graphic>
643 Some of the classes shown here (the ones sitting in white boxes) are
644 abstract base classes that are of no interest at all for the user.
645 They are used internally in order to avoid code duplication if
646 two or more classes derived from them share certain features. An
647 example would be <literal>expairseq</literal>, which is a container
648 for a sequence of pairs each consisting of one expression and a number
649 (<literal>numeric</literal>). What <emphasis>is</emphasis> visible to
650 the user are the derived classes <literal>add</literal> and
651 <literal>mul</literal>, representing sums of terms and products,
652 respectively. We'll come back later to some more details about these
653 two classes and motivate the use of pairs in sums and products here.</para>
655 <sect2><title>Digression: Internal representation of products and sums</title>
657 <para>Although it should be completely transparent for the user of
658 GiNaC a short discussion of this topic helps understanding the sources
659 and also explains performance to a large degree. Consider the
660 symbolic expression <literal>(a+2*b-c)*d</literal>, which could
661 naively be represented by a tree of linear containers for addition and
662 multiplication together with atomic leaves of symbols and integer
663 numbers in this fashion:
664 <figure id="repres-naive-id" float="1">
665 <title>Naive internal representation of <emphasis>d(a+2*b-c)</emphasis></title>
666 <graphic align="center" fileref="rep_naive.graext" format="GRAEXT"></graphic>
668 However, doing so results in a rather deeply nested tree which will
669 quickly become rather slow to manipulate. If we represent the sum
670 instead as a sequence of terms, each having a purely numeric
671 coefficient, the tree becomes much more flat.
672 <figure id="repres-pair-id" float="1">
673 <title>Better internal representation of <emphasis>d(a+2*b-c)</emphasis></title>
674 <graphic align="center" fileref="rep_pair.graext" format="GRAEXT"></graphic>
676 The number <literal>1</literal> under the symbol <literal>d</literal>
677 is a hint that multiplication objects can be treated similarly where
678 the coeffiecients are interpreted as <emphasis>exponents</emphasis>
679 now. Addition of sums of terms or multiplication of products with
680 numerical exponents can be made very efficient with a
681 pair-representation. Internally, this handling is done by most CAS in
682 this way. It typically speeds up manipulations by an order of
683 magnitude. Now it should be clear, why both classes
684 <literal>add</literal> and <literal>mul</literal> are derived from the
685 same abstract class: the representation is the same, only the
686 interpretation differs. </para>
690 <sect1><title>Symbols</title>
692 <para>Symbols are for symbolic manipulation what atoms are for
693 chemistry. You can declare objects of class <literal>symbol</literal>
694 as any other object simply by saying <literal>symbol x,y;</literal>.
695 There is, however, a catch in here having to do with the fact that C++
696 is a compiled language. The information about the symbol's name is
697 thrown away by the compiler but at a later stage you may want to print
698 expressions holding your symbols. In order to avoid confusion GiNaC's
699 symbols are able to know their own name. This is accomplished by
700 declaring its name for output at construction time in the fashion
701 <literal>symbol x("x");</literal>. If you declare a symbol using the
702 default constructor (i.e. without string-argument) the system will
703 deal out a unique name. That name may not be suitable for printing
704 but for internal routines when no output is desired it is often
705 enough. We'll come across examples of such symbols later in this
708 <para>Although symbols can be assigned expressions for internal
709 reasons, you should not do it (and we are not going to tell you how it
710 is done). If you want to replace a symbol with something else in an
711 expression, you can use the expression's <literal>.subs()</literal>
716 <sect1><title>Numbers</title>
718 <para>For storing numerical things, GiNaC uses Bruno Haible's library
719 <literal>CLN</literal>. The classes therein serve as foundation
720 classes for GiNaC. <literal>CLN</literal> stands for Class Library
721 for Numbers or alternatively for Common Lisp Numbers. In order to
722 find out more about <literal>CLN</literal>'s internals the reader is
723 refered to the documentation of that library. Suffice to say that it
724 is by itself build on top of another library, the GNU Multiple
725 Precision library <literal>GMP</literal>, which is an extremely fast
726 library for arbitrary long integers and rationals as well as arbitrary
727 precision floating point numbers. It is very commonly used by several
728 popular cryptographic applications. <literal>CLN</literal> extends
729 <literal>GMP</literal> by several useful things: First, it introduces
730 the complex number field over either reals (i.e. floating point
731 numbers with arbitrary precision) or rationals. Second, it
732 automatically converts rationals to integers if the denominator is
733 unity and complex numbers to real numbers if the imaginary part
734 vanishes and also correctly treats algebraic functions. Third it
735 provides good implementations of state-of-the-art algorithms for all
736 trigonometric and hyperbolic functions as well as for calculation of
737 some useful constants.</para>
739 <para>The user can construct an object of class
740 <literal>numeric</literal> in several ways. The following example
741 shows the four most important constructors: construction from
742 C-integer, construction of fractions from two integers, construction
743 from C-float and construction from a string.
744 <example><title>Construction of numbers</title>
746 #include <ginac/ginac.h>
747 using namespace GiNaC;
751 numeric two(2); // exact integer 2
752 numeric r(2,3); // exact fraction 2/3
753 numeric e(2.71828); // floating point number
754 numeric p("3.1415926535897932385"); // floating point number
756 cout << two*p << endl; // floating point 6.283...
761 Note that all those constructors are <emphasis>explicit</emphasis>
762 which means you are not allowed to write <literal>numeric
763 two=2;</literal>. This is because the basic objects to be handled by
764 GiNaC are the expressions <literal>ex</literal> and we want to keep
765 things simple and wish objects like <literal>pow(x,2)</literal> to be
766 handled the same way as <literal>pow(x,a)</literal>, which means that
767 we need to allow a general <literal>ex</literal> as base and exponent.
768 Therefore there is an implicit constructor from C-integers directly to
769 expressions handling numerics at work in most of our examples. This
770 design really becomes convenient when one declares own functions
771 having more than one parameter but it forbids using implicit
772 constructors because that would lead to ambiguities. </para>
774 <para>It may be tempting to construct numbers writing <literal>numeric
775 r(3/2)</literal>. This would, however, call C's built-in operator
776 <literal>/</literal> for integers first and result in a numeric
777 holding a plain integer 1. <emphasis>Never use
778 </emphasis><literal>/</literal><emphasis> on integers!</emphasis> Use
779 the constructor from two integers instead, as shown in the example
780 above. Writing <literal>numeric(1)/2</literal> may look funny but
783 <para>We have seen now the distinction between exact numbers and
784 floating point numbers. Clearly, the user should never have to worry
785 about dynamically created exact numbers, since their "exactness"
786 always determines how they ought to be handled. The situation is
787 different for floating point numbers. Their accuracy is handled by
788 one <emphasis>global</emphasis> variable, called
789 <literal>Digits</literal>. (For those readers who know about Maple:
790 it behaves very much like Maple's <literal>Digits</literal>). All
791 objects of class numeric that are constructed from then on will be
792 stored with a precision matching that number of decimal digits:
793 <example><title>Controlling the precision of floating point numbers</title>
795 #include <ginac/ginac.h>
796 using namespace GiNaC;
800 numeric three(3.0), one(1.0);
801 numeric x = one/three;
803 cout << "in " << Digits << " digits:" << endl;
804 cout << x << endl;
805 cout << Pi.evalf() << endl;
817 The above example prints the following output to screen:
823 0.333333333333333333333333333333333333333333333333333333333333333333
824 3.14159265358979323846264338327950288419716939937510582097494459231
828 <sect2><title>Tests on numbers</title>
830 <para>Once you have declared some numbers, assigned them to
831 expressions and done some arithmetic with them it is frequently
832 desired to retrieve some kind of information from them like asking
833 whether that number is integer, rational, real or complex. For those
834 cases GiNaC provides several useful methods. (Internally, they fall
835 back to invocations of certain CLN functions.)</para>
837 <para>As an example, let's construct some rational number, multiply it
838 with some multiple of its denominator and check what comes out:
839 <example><title>Sample test on objects of type numeric</title>
841 #include <ginac/ginac.h>
842 using namespace GiNaC;
846 numeric twentyone(21);
848 numeric answer(21,5);
850 cout << answer.is_integer() << endl; // false, it's 21/5
852 cout << answer.is_integer() << endl; // true, it's 42 now!
858 Note that the variable <literal>answer</literal> is constructed here
859 as an integer but in an intermediate step it holds a rational number
860 represented as integer numerator and integer denominator. When
861 multiplied by 10, the denominator becomes unity and the result is
862 automatically converted to a pure integer again. Internally, the
863 underlying <literal>CLN</literal> is responsible for this behaviour
864 and we refer the reader to <literal>CLN</literal>'s documentation.
865 Suffice to say that the same behaviour applies to complex numbers as
866 well as return values of certain functions. Complex numbers are
867 automatically converted to real numbers if the imaginary part becomes
868 zero. The full set of tests that can be applied is listed in the
871 <informaltable colsep="0" frame="topbot" pgwide="1">
873 <colspec colnum="1" colwidth="1*">
874 <colspec colnum="2" colwidth="2*">
877 <entry>Method</entry>
878 <entry>Returns true if...</entry>
883 <entry><literal>.is_zero()</literal></entry>
884 <entry>object is equal to zero</entry>
887 <entry><literal>.is_positive()</literal></entry>
888 <entry>object is not complex and greater than 0</entry>
891 <entry><literal>.is_integer()</literal></entry>
892 <entry>object is a (non-complex) integer</entry>
895 <entry><literal>.is_pos_integer()</literal></entry>
896 <entry>object is an integer and greater than 0</entry>
899 <entry><literal>.is_nonneg_integer()</literal></entry>
900 <entry>object is an integer and greater equal 0</entry>
903 <entry><literal>.is_even()</literal></entry>
904 <entry>object is an even integer</entry>
907 <entry><literal>.is_odd()</literal></entry>
908 <entry>object is an odd integer</entry>
911 <entry><literal>.is_prime()</literal></entry>
912 <entry>object is a prime integer (probabilistic primality test)</entry>
915 <entry><literal>.is_rational()</literal></entry>
916 <entry>object is an exact rational number (integers are rational, too, as are complex extensions like <literal>2/3+7/2*I</literal>)</entry>
919 <entry><literal>.is_real()</literal></entry>
920 <entry>object is a real integer, rational or float (i.e. is not complex)</entry>
932 <sect1><title>Constants</title>
934 <para>Constants behave pretty much like symbols except that that they return
935 some specific number when the method <literal>.evalf()</literal> is called.
938 <para>The predefined known constants are:
939 <informaltable colsep="0" frame="topbot" pgwide="1">
941 <colspec colnum="1" colwidth="1*">
942 <colspec colnum="2" colwidth="2*">
943 <colspec colnum="3" colwidth="4*">
947 <entry>Common Name</entry>
948 <entry>Numerical Value (35 digits)</entry>
953 <entry><literal>Pi</literal></entry>
954 <entry>Archimedes' constant</entry>
955 <entry>3.14159265358979323846264338327950288</entry>
957 <entry><literal>Catalan</literal></entry>
958 <entry>Catalan's constant</entry>
959 <entry>0.91596559417721901505460351493238411</entry>
961 <entry><literal>EulerGamma</literal></entry>
962 <entry>Euler's (or Euler-Mascheroni) constant</entry>
963 <entry>0.57721566490153286060651209008240243</entry>
972 <sect1><title>Fundamental operations: The <literal>power</literal>, <literal>add</literal> and <literal>mul</literal> classes</title>
974 <para>Simple polynomial expressions are written down in GiNaC pretty
975 much like in other CAS. The necessary operators <literal>+</literal>,
976 <literal>-</literal>, <literal>*</literal> and <literal>/</literal>
977 have been overloaded to achieve this goal. When you run the following
978 program, the constructor for an object of type <literal>mul</literal>
979 is automatically called to hold the product of <literal>a</literal>
980 and <literal>b</literal> and then the constructor for an object of
981 type <literal>add</literal> is called to hold the sum of that
982 <literal>mul</literal> object and the number one:
983 <example><title>Construction of <literal>add</literal> and <literal>mul</literal> objects</title>
985 #include <ginac/ginac.h>
986 using namespace GiNaC;
990 symbol a("a"), b("b");
997 <para>For exponentiation, you have already seen the somewhat clumsy
998 (though C-ish) statement <literal>pow(x,2);</literal> to represent
999 <literal>x</literal> squared. This direct construction is necessary
1000 since we cannot safely overload the constructor <literal>^</literal>
1001 in <literal>C++</literal> to construct a <literal>power</literal>
1002 object. If we did, it would have several counterintuitive effects:
1005 <para>Due to <literal>C</literal>'s operator precedence,
1006 <literal>2*x^2</literal> would be parsed as <literal>(2*x)^2</literal>.
1009 <para>Due to the binding of the operator <literal>^</literal>,
1010 <literal>x^a^b</literal> would result in <literal>(x^a)^b</literal>.
1011 This would be confusing since most (though not all) other CAS
1012 interpret this as <literal>x^(a^b)</literal>.
1015 <para>Also, expressions involving integer exponents are very
1016 frequently used, which makes it even more dangerous to overload
1017 <literal>^</literal> since it is then hard to distinguish between the
1018 semantics as exponentiation and the one for exclusive or. (It would
1019 be embarassing to return <literal>1</literal> where one has requested
1020 <literal>2^3</literal>.)</para>
1023 All effects are contrary to mathematical notation and differ from the
1024 way most other CAS handle exponentiation, therefore overloading
1025 <literal>^</literal> is ruled out for GiNaC's C++ part. The situation
1026 is different in <literal>ginsh</literal>, there the
1027 exponentiation-<literal>^</literal> exists. (Also note, that the
1028 other frequently used exponentiation operator <literal>**</literal>
1029 does not exist at all in <literal>C++</literal>).</para>
1031 <para>To be somewhat more precise, objects of the three classes
1032 described here, are all containers for other expressions. An object
1033 of class <literal>power</literal> is best viewed as a container with
1034 two slots, one for the basis, one for the exponent. All valid GiNaC
1035 expressions can be inserted. However, basic transformations like
1036 simplifying <literal>pow(pow(x,2),3)</literal> to
1037 <literal>x^6</literal> automatically are only performed when
1038 this is mathematically possible. If we replace the outer exponent
1039 three in the example by some symbols <literal>a</literal>, the
1040 simplification is not safe and will not be performed, since
1041 <literal>a</literal> might be <literal>1/2</literal> and
1042 <literal>x</literal> negative.</para>
1044 <para>Objects of type <literal>add</literal> and
1045 <literal>mul</literal> are containers with an arbitrary number of
1046 slots for expressions to be inserted. Again, simple and safe
1047 simplifications are carried out like transforming
1048 <literal>3*x+4-x</literal> to <literal>2*x+4</literal>.</para>
1050 <para>The general rule is that when you construct such objects, GiNaC
1051 automatically creates them in canonical form, which might differ from
1052 the form you typed in your program. This allows for rapid comparison
1053 of expressions, since after all <literal>a-a</literal> is simply zero.
1054 Note, that the canonical form is not necessarily lexicographical
1055 ordering or in any way easily guessable. It is only guaranteed that
1056 constructing the same expression twice, either implicitly or
1057 explicitly, results in the same canonical form.</para>
1061 <sect1><title>Built-in Functions</title>
1063 <para>This section is not here yet</para>
1073 <title>Important Algorithms</title>
1075 <para>In this chapter the most important algorithms provided by GiNaC
1076 will be described. Some of them are implemented as functions on
1077 expressions, others are implemented as methods provided by expression
1078 objects. If they are methods, there exists a wrapper function around
1079 it, so you can alternatively call it in a functional way as shown in
1081 <example><title>Methods vs. wrapper functions</title>
1083 #include <ginac/ginac.h>
1084 using namespace GiNaC;
1088 ex x = numeric(1.0);
1090 cout << "As method: " << sin(x).evalf() << endl;
1091 cout << "As function: " << evalf(sin(x)) << endl;
1096 The general rule is that wherever methods accept one or more
1097 parameters (<emphasis>arg1</emphasis>, <emphasis>arg2</emphasis>, ...)
1098 the order of arguments the function wrapper accepts is the same but
1099 preceded by the object to act on (<emphasis>object</emphasis>,
1100 <emphasis>arg1</emphasis>, <emphasis>arg2</emphasis>, ...). This
1101 approach is the most natural one in an OO model but it may lead to
1102 confusion for MapleV users because where they would type
1103 <literal>A:=x+1; subs(x=2,A);</literal> GiNaC would require
1104 <literal>A=x+1; subs(A,x==2);</literal> (after proper declaration of A
1105 and x). On the other hand, since MapleV returns 3 on
1106 <literal>A:=x^2+3; coeff(A,x,0);</literal> (GiNaC:
1107 <literal>A=pow(x,2)+3; coeff(A,x,0);</literal>) it is clear that
1108 MapleV is not trying to be consistent here. Also, users of MuPAD will
1109 in most cases feel more comfortable with GiNaC's convention. All
1110 function wrappers are always implemented as simple inline functions
1111 which just call the corresponding method and are only provided for
1112 users uncomfortable with OO who are dead set to avoid method
1113 invocations. Generally, a chain of function wrappers is much harder
1114 to read than a chain of methods and should therefore be avoided if
1115 possible. On the other hand, not everything in GiNaC is a method on
1116 class <literal>ex</literal> and sometimes calling a function cannot be
1120 <sect1><title>Polynomial Expansion</title>
1122 <para>A polynomial in one or more variables has many equivalent
1123 representations. Some useful ones serve a specific purpose. Consider
1124 for example the trivariate polynomial <literal>4*x*y + x*z + 20*y^2 +
1125 21*y*z + 4*z^2</literal>. It is equivalent to the factorized
1126 polynomial <literal>(x + 5*y + 4*z)*(4*y + z)</literal>. Other
1127 representations are the recursive ones where one collects for
1128 exponents in one of the three variable. Since the factors are
1129 themselves polynomials in the remaining two variables the procedure
1130 can be repeated. In our expample, two possibilies would be
1131 <literal>(4*y + z)*x + 20*y^2 + 21*y*z + 4*z^2</literal> and
1132 <literal>20*y^2 + (21*z + 4*x)*y + 4*z^2 + x*z</literal>.
1135 <para>To bring an expression into expanded form, its method
1136 <function>.expand()</function> may be called. In our example above,
1137 this corresponds to <literal>4*x*y + x*z + 20*y^2 + 21*y*z +
1138 4*z^2</literal>. Again, since the canonical form in GiNaC is not
1139 easily guessable you should be prepared to see different orderings of
1140 terms in such sums!</para>
1144 <sect1><title>Collecting expressions</title>
1146 <para>Another useful representation of multivariate polynomials is as
1147 a univariate polynomial in one of the variables with the coefficients
1148 being polynomials in the remaining variables. The method
1149 <literal>collect()</literal> accomplishes this task:
1151 <funcsynopsisinfo>#include <ginac/ginac.h></funcsynopsisinfo>
1152 <funcdef>ex <function>ex::collect</function></funcdef>
1153 <paramdef>symbol const & <parameter>s</parameter></paramdef>
1155 Note that the original polynomial needs to be in expanded form in
1156 order to be able to find the coefficients properly. The range of
1157 occuring coefficients can be checked using the two methods
1159 <funcsynopsisinfo>#include <ginac/ginac.h></funcsynopsisinfo>
1160 <funcdef>int <function>ex::degree</function></funcdef>
1161 <paramdef>symbol const & <parameter>s</parameter></paramdef>
1164 <funcdef>int <function>ex::ldegree</function></funcdef>
1165 <paramdef>symbol const & <parameter>s</parameter></paramdef>
1167 where <literal>degree()</literal> returns the highest coefficient and
1168 <literal>ldegree()</literal> the lowest one. These two methods work
1169 also reliably on non-expanded input polynomials. This is illustrated
1170 in the following example:
1172 <example><title>Collecting expressions in multivariate polynomials</title>
1174 #include <ginac/ginac.h>
1175 using namespace GiNaC;
1179 symbol x("x"), y("y");
1180 ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
1181 - pow(x+y,2) + 2*pow(y+2,2) - 8;
1182 ex Poly = PolyInp.expand();
1184 for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) {
1185 cout << "The x^" << i << "-coefficient is "
1186 << Poly.coeff(x,i) << endl;
1188 cout << "As polynomial in y: "
1189 << Poly.collect(y) << endl;
1194 When run, it returns an output in the following fashion:
1196 The x^0-coefficient is y^2+11*y
1197 The x^1-coefficient is 5*y^2-2*y
1198 The x^2-coefficient is -1
1199 The x^3-coefficient is 4*y
1200 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
1202 As always, the exact output may vary between different versions of
1203 GiNaC or even from run to run since the internal canonical ordering is
1204 not within the user's sphere of influence.</para>
1208 <sect1><title>Polynomial Arithmetic</title>
1210 <sect2><title>GCD and LCM</title>
1212 <para>The functions for polynomial greatest common divisor and least common
1213 multiple have the synopsis:
1215 <funcsynopsisinfo>#include <GiNaC/normal.h></funcsynopsisinfo>
1216 <funcdef>ex <function>gcd</function></funcdef>
1217 <paramdef>const ex *<parameter>a</parameter>, const ex *<parameter>b</parameter></paramdef>
1220 <funcdef>ex <function>lcm</function></funcdef>
1221 <paramdef>const ex *<parameter>a</parameter>, const ex *<parameter>b</parameter></paramdef>
1222 </funcsynopsis></para>
1224 <para>The functions <function>gcd()</function> and <function
1225 id="lcm-main">lcm()</function> accepts two expressions
1226 <literal>a</literal> and <literal>b</literal> as arguments and return
1227 a new expression, their greatest common divisor or least common
1228 multiple, respectively. If the polynomials <literal>a</literal> and
1229 <literal>b</literal> are coprime <function>gcd(a,b)</function> returns 1
1230 and <function>lcm(a,b)</function> returns the product of
1231 <literal>a</literal> and <literal>b</literal>.
1232 <example><title>Polynomal GCD/LCM</title>
1234 #include <ginac/ginac.h>
1235 using namespace GiNaC;
1239 symbol x("x"), y("y"), z("z");
1240 ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
1241 ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
1243 ex P_gcd = gcd(P_a, P_b);
1245 ex P_lcm = lcm(P_a, P_b);
1246 // 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
1255 <sect2><title>The <function>normal</function> method</title>
1257 <para>While in common symbolic code <function>gcd()</function> and
1258 <function>lcm()</function> are not too heavily used, simplification
1259 occurs frequently. Therefore <function>.normal()</function>, which
1260 provides some basic form of simplification, has become a method of
1261 class <literal>ex</literal>, just like <literal>.expand()</literal>.
1262 It converts a rational function into an equivalent rational function
1263 where numererator and denominator are coprime. This means, it finds
1264 the GCD of numerator and denominator and cancels it. If it encounters
1265 some object which does not belong to the domain of rationals (a
1266 function for instance), that object is replaced by a temporary symbol.
1267 This means that both expressions <literal>t1</literal> and
1268 <literal>t2</literal> are indeed simplified in this little program:
1269 <example><title>Cancellation of polynomial GCD (with obstacles)</title>
1271 #include <ginac/ginac.h>
1272 using namespace GiNaC;
1277 ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
1278 ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
1279 cout << "t1 is " << t1.normal() << endl;
1280 cout << "t2 is " << t2.normal() << endl;
1286 Of course this works for multivariate polynomials too, so the ratio of
1287 the sample-polynomials from the section about GCD and LCM above would
1288 be normalized to <literal>P_a/P_b</literal> =
1289 <literal>(4*y+z)/(y+3*z)</literal>.</para>
1295 <sect1><title>Symbolic Differentiation</title>
1297 <para>GiNaC's objects know how to differentiate themselves. Thus, a
1298 polynomial (class <literal>add</literal>) knows that its derivative is
1299 the sum of the derivatives of all the monomials:
1300 <example><title>Simple polynomial differentiation</title>
1302 #include <ginac/ginac.h>
1303 using namespace GiNaC;
1307 symbol x("x"), y("y"), z("z");
1308 ex P = pow(x, 5) + pow(x, 2) + y;
1310 cout << P.diff(x,2) << endl; // 20*x^3 + 2
1311 cout << P.diff(y) << endl; // 1
1312 cout << P.diff(z) << endl; // 0
1317 If a second integer parameter <emphasis>n</emphasis> is given the
1318 <function>diff</function> method returns the <emphasis>n</emphasis>th
1321 <para>If <emphasis>every</emphasis> object and every function is told
1322 what its derivative is, all derivatives of composed objects can be
1323 calculated using the chain rule and the product rule. Consider, for
1324 instance the expression <literal>1/cosh(x)</literal>. Since the
1325 derivative of <literal>cosh(x)</literal> is <literal>sinh(x)</literal>
1326 and the derivative of <literal>pow(x,-1)</literal> is
1327 <literal>-pow(x,-2)</literal> GiNaC can readily compute the
1328 composition. It turns out that the composition is the generating
1329 function for Euler Numbers, i.e. the so called
1330 <emphasis>n</emphasis>th Euler number is the coefficient of
1331 <literal>x^n/n!</literal> in the expansion of
1332 <literal>1/cosh(x)</literal>. We may use this identity to code a
1333 function that generates Euler numbers in just three lines:
1334 <example><title>Differentiation with nontrivial functions: Euler numbers</title>
1336 #include <ginac/ginac.h>
1337 using namespace GiNaC;
1339 ex EulerNumber(unsigned n)
1342 ex generator = pow(cosh(x),-1);
1343 return generator.diff(x,n).subs(x==0);
1348 for (unsigned i=0; i<11; i+=2)
1349 cout << EulerNumber(i) << endl;
1354 When you run it, it produces the sequence <literal>1</literal>,
1355 <literal>-1</literal>, <literal>5</literal>, <literal>-61</literal>,
1356 <literal>1385</literal>, <literal>-50521</literal>. We increment the
1357 loop variable <literal>i</literal> by two since all odd Euler numbers
1358 vanish anyways.</para>
1362 <sect1><title>Series Expansion</title>
1364 <para>Expressions know how to expand themselves as a Taylor series or
1365 (more generally) a Laurent series. As in most conventional Computer
1366 Algebra Systems no distinction is made between those two. There is a
1367 class of its own for storing such series as well as a class for
1368 storing the order of the series. A sample program could read:
1369 <example><title>Series expansion</title>
1371 #include <ginac/ginac.h>
1372 using namespace GiNaC;
1378 ex MyExpr1 = sin(x);
1379 ex MyExpr2 = 1/(x - pow(x, 2) - pow(x, 3));
1380 ex MyTailor, MySeries;
1382 MyTailor = MyExpr1.series(x, point, 5);
1383 cout << MyExpr1 << " == " << MyTailor
1384 << " for small " << x << endl;
1385 MySeries = MyExpr2.series(x, point, 7);
1386 cout << MyExpr2 << " == " << MySeries
1387 << " for small " << x << endl;
1394 <para>As an instructive application, let us calculate the numerical
1395 value of Archimedes' constant (for which there already exists the
1396 built-in constant <literal>Pi</literal>) using Méchain's
1397 wonderful formula <literal>Pi==16*atan(1/5)-4*atan(1/239)</literal>.
1398 We may expand the arcus tangent around <literal>0</literal> and insert
1399 the fractions <literal>1/5</literal> and <literal>1/239</literal>.
1400 But, as we have seen, a series in GiNaC carries an order term with it.
1401 The function <literal>series_to_poly</literal> may be used to strip
1403 <example><title>Series expansion using Méchain's formula for
1404 <literal>Pi</literal></title>
1406 #include <ginac/ginac.h>
1407 using namespace GiNaC;
1409 ex mechain_pi(int degr)
1412 ex pi_expansion = series_to_poly(atan(x).series(x,0,degr));
1413 ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
1414 -4*pi_expansion.subs(x==numeric(1,239));
1421 for (int i=2; i<12; i+=2) {
1422 pi_frac = mechain_pi(i);
1423 cout << i << ":\t" << pi_frac << endl
1424 << "\t" << pi_frac.evalf() << endl;
1429 <para>When you run this program, it will type out:</para>
1432 3.1832635983263598326
1433 4: 5359397032/1706489875
1434 3.1405970293260603143
1435 6: 38279241713339684/12184551018734375
1436 3.141621029325034425
1437 8: 76528487109180192540976/24359780855939418203125
1438 3.141591772182177295
1439 10: 327853873402258685803048818236/104359128170408663038552734375
1440 3.1415926824043995174
1451 <title>Extending GiNaC</title>
1453 <para>By reading so far you should have gotten a fairly good
1454 understanding of GiNaC's design-patterns. From here on you should
1455 start reading the sources. All we can do now is issue some
1456 recommendations how to tackle the many open ends the system still
1457 has in order to fulfill everybody's dreams.</para>
1459 <sect1><title>What doesn't belong into GiNaC</title>
1461 <para>First of all, GiNaC's name must be read literally. It is
1462 designed to be a library for use within C++. The tiny
1463 <literal>ginsh</literal> accompanying GiNaC makes this even more
1464 clear: it doesn't even attempt to provide a language. There are no
1465 loops or conditional expressions in <literal>ginsh</literal>, it is
1466 merely a window into the library for the programmer to test stuff (or
1467 to show off). Still, the design of a complete CAS with a language of
1468 its own, graphical capabilites and all this on top of GiNaC is
1469 possible and is without doubt a nice project for the future.</para>
1471 <para>There are many built-in functions in GiNaC that do not know how
1472 to evaluate themselves numerically to a precision declared at runtime
1473 (using <literal>Digits</literal>). Some may be evaluated at certain
1474 points, but not generally. This ought to be fixed. However, doing
1475 numerical computations with GiNaC's quite abstract classes is doomed
1476 to be inefficient. For this purpose, the underlying bignum-package
1477 <literal>CLN</literal> is much better suited.</para>
1481 <sect1><title>Other symbolic functions</title>
1483 <para>The easiest and most instructive way to start with is probably
1484 to implement your own function. Objects of class
1485 <literal>function</literal> are inserted into the system via a kind of
1486 "registry". They get a serial number that is used internally to
1487 identify them but you usually need not worry about this. What you
1488 have to care for are functions that are called when the user invokes
1489 certain methods. These are usual C++-functions accepting a number of
1490 <literal>ex</literal> as arguments and returning one
1491 <literal>ex</literal>. As an example, if we have a look at a
1492 simplified implementation of the cosine trigonometric function, we
1493 first need a function that is called when one wishes to
1494 <literal>eval</literal> it. It could look something like this:
1497 static ex cos_eval_method(ex const & x)
1499 // if x%2*Pi return 1
1500 // if x%Pi return -1
1501 // if x%Pi/2 return 0
1502 // care for other cases...
1503 return cos(x).hold();
1506 The last line returns <literal>cos(x)</literal> if we don't know what
1507 else to do and stops a potential recursive evaluation by saying
1508 <literal>.hold()</literal>. We should also implement a method for
1509 numerical evaluation and since we are lazy we sweep the problem under
1510 the rug by calling someone else's function that does so, in this case
1511 the one in class <literal>numeric</literal>:
1513 static ex cos_evalf_method(ex const & x)
1515 return sin(ex_to_numeric(x));
1518 Differentiation will surely turn up and so we need to tell
1519 <literal>sin</literal> how to differentiate itself:
1521 static ex cos_diff_method(ex const & x, unsigned diff_param)
1527 The second parameter is obligatory but uninteresting at this point.
1528 It is used for correct handling of the product rule only. For Taylor
1529 expansion, it is enough to know how to differentiate. But if the
1530 function you want to implement does have a pole somewhere in the
1531 complex plane, you need to write another method for Laurent expansion
1532 around that point.</para>
1534 <para>Now that everything has been written for <literal>cos</literal>,
1535 we need to tell the system about it. This is done by a macro and we
1536 are not going to descibe how it expands, please consult your
1537 preprocessor if you are curious:
1539 REGISTER_FUNCTION(cos, cos_eval_method, cos_evalf_method, cos_diff, NULL);
1541 The first argument is the function's name, the second, third and
1542 fourth bind the corresponding methods to this objects and the fifth is
1543 a slot for inserting a method for series expansion. Also, the new
1544 function needs to be declared somewhere. This may also be done by a
1545 convenient preprocessor macro:
1547 DECLARE_FUNCTION_1P(cos)
1549 The suffix <literal>_1P</literal> stands for <emphasis>one
1550 parameter</emphasis>. Of course, this implementation of
1551 <literal>cos</literal> is very incomplete and lacks several safety
1552 mechanisms. Please, have a look at the real implementation in GiNaC.
1553 (By the way: in case you are worrying about all the macros above we
1554 can assure you that functions are GiNaC's most macro-intense classes.
1555 We have done our best to avoid them where we can.)</para>
1557 <para>That's it. May the source be with you!</para>
1565 <title>A Comparison with other CAS</title>
1567 <para>This chapter will give you some information on how GiNaC
1568 compares to other, traditional Computer Algebra Systems, like
1569 <literal>Maple</literal>, <literal>Mathematica</literal> or
1570 <literal>Reduce</literal>, where it has advantages and disadvantages
1571 over these systems.</para>
1573 <sect1><title>Advantages</title>
1575 <para>GiNaC has several advantages over traditional Computer
1576 Algebra Systems, like
1580 <para>familiar language: all common CAS implement their own
1581 proprietary grammar which you have to learn first (and maybe learn
1582 again when your vendor chooses to "enhance" it). With GiNaC you
1583 can write your program in common <literal>C++</literal>, which is
1584 standardized.</para>
1587 <para>structured data types: you can build up structured data
1588 types using <literal>struct</literal>s or <literal>class</literal>es
1589 together with STL features instead of using unnamed lists of lists
1593 <para>strongly typed: in CAS, you usually have only one kind of
1594 variables which can hold contents of an arbitrary type. This
1595 4GL like feature is nice for novice programmers, but dangerous.
1599 <para>development tools: powerful development tools exist for
1600 <literal>C++</literal>, like fancy editors (e.g. with automatic
1601 indentation and syntax highlighting), debuggers, visualization
1602 tools, documentation tools...</para>
1605 <para>modularization: <literal>C++</literal> programs can
1606 easily be split into modules by separating interface and
1607 implementation.</para>
1610 <para>price: GiNaC is distributed under the GNU Public License
1611 which means that it is free and available with source code. And
1612 there are excellent <literal>C++</literal>-compilers for free, too.
1616 <para>extendable: you can add your own classes to GiNaC, thus
1617 extending it on a very low level. Compare this to a traditional
1618 CAS that you can usually only extend on a high level by writing in
1619 the language defined by the parser. In particular, it turns out
1620 to be almost impossible to fix bugs in a traditional system.
1623 <para>seemless integration: it is somewhere between difficult
1624 and impossible to call CAS functions from within a program
1625 written in <literal>C++</literal> or any other programming
1626 language and vice versa. With GiNaC, your symbolic routines
1627 are part of your program. You can easily call third party
1628 libraries, e.g. for numerical evaluation or graphical
1629 interaction. All other approaches are much more cumbersome: they
1630 range from simply ignoring the problem
1631 (i.e. <literal>Maple</literal>) to providing a
1632 method for "embedding" the system
1633 (i.e. <literal>Yacas</literal>).</para>
1636 <para>efficiency: often large parts of a program do not need
1637 symbolic calculations at all. Why use large integers for loop
1638 variables or arbitrary precision arithmetics where double
1639 accuracy is sufficient? For pure symbolic applications,
1640 GiNaC is comparable in speed with other CAS.
1645 <sect1><title>Disadvantages</title>
1647 <para>Of course it also has some disadvantages
1651 <para>not interactive: GiNaC programs have to be written in
1652 an editor, compiled and executed. You cannot play with
1653 expressions interactively. However, such an extension is not
1654 inherently forbidden by design. In fact, two interactive
1655 interfaces are possible: First, a simple shell that exposes GiNaC's
1656 types to a command line can readily be written (and has been
1657 written) and second, as a more consistent approach we plan
1658 an integration with the <literal>CINT</literal>
1659 <literal>C++</literal> interpreter.</para>
1662 <para>advanced features: GiNaC cannot compete with a program
1663 like <literal>Reduce</literal> which exists for more than
1664 30 years now or <literal>Maple</literal> which grows since
1665 1981 by the work of dozens of programmers, with respect to
1666 mathematical features. Integration, factorization, non-trivial
1667 simplifications, limits etc. are missing in GiNaC (and are not
1668 planned for the near future).</para>
1671 <para>portability: While the GiNaC library itself is designed
1672 to avoid any platform dependent features (it should compile
1673 on any ANSI compliant <literal>C++</literal> compiler), the
1674 currently used version of the CLN library (fast large integer and
1675 arbitrary precision arithmetics) can be compiled only on systems
1676 with a recently new <literal>C++</literal> compiler from the
1677 GNU Compiler Collection (<literal>GCC</literal>). GiNaC uses
1678 recent language features like explicit constructors, mutable
1679 members, RTTI, dynamic_casts and STL, so ANSI compliance is meant
1680 literally. Recent <literal>GCC</literal> versions starting at
1681 2.95, although itself not yet ANSI compliant, support all needed
1688 <sect1><title>Why <literal>C++</literal>?</title>
1690 <para>Why did we choose to implement GiNaC in <literal>C++</literal>
1691 instead of <literal>Java</literal> or any other language?
1692 <literal>C++</literal> is not perfect: type checking is not strict
1693 (casting is possible), separation between interface and implementation
1694 is not complete, object oriented design is not enforced. The main
1695 reason is the often scolded feature of operator overloading in
1696 <literal>C++</literal>. While it may be true that operating on classes
1697 with a <literal>+</literal> operator is rarely meaningful, it is
1698 perfectly suited for algebraic expressions. Writing 3x+5y as
1699 <literal>3*x+5*y</literal> instead of
1700 <literal>x.times(3).plus(y.times(5))</literal> looks much more
1701 natural. Furthermore, the main developers are more familiar with
1702 <literal>C++</literal> than with any other programming
1713 <title>ISO/IEC 14882:1998</title>
1714 <subtitle>Programming Languages: C++</subtitle>
1719 <title>CLN: A Class Library for Numbers</title>
1722 <firstname>Bruno</firstname><surname>Haible</surname>
1723 <affiliation><address><email>haible@ilog.fr</email></address></affiliation>
1730 <title>The C++ Programming Language</title>
1731 <authorgroup><author><firstname>Bjarne</firstname><surname>Stroustrup</surname></author></authorgroup>
1732 <edition>3</edition>
1733 <isbn>0-201-88954-4</isbn>
1734 <publisher><publishername>Addison Wesley</publishername></publisher>
1740 <title>C++ FAQs</title>
1741 <authorgroup><author><firstname>Marshall</firstname><surname>Cline</surname></author></authorgroup>
1742 <isbn>0-201-58958-3</isbn>
1743 <pubdate>1995</pubdate>
1744 <publisher><publishername>Addison Wesley</publishername></publisher>
1750 <title>Algorithms for Computer Algebra</title>
1752 <author><firstname>Keith</firstname><othername>O.</othername><surname>Geddes</surname></author>
1753 <author><firstname>Stephen</firstname><othername>R.</othername><surname>Czapor</surname></author>
1754 <author><firstname>George</firstname><surname>Labahn</surname></author>
1756 <isbn>0-7923-9259-0</isbn>
1757 <pubdate>1992</pubdate>
1759 <publishername>Kluwer Academic Publishers</publishername>
1760 <address><city>Norwell</city>, <state>Massachusetts</state></address>
1767 <title>Computer Algebra</title>
1768 <subtitle>Systems and Algorithms for Algebraic Computation</subtitle>
1770 <author><firstname>J.</firstname><othername>H.</othername><surname>Davenport</surname></author>
1771 <author><firstname>Y.</firstname><surname>Siret</surname></author>
1772 <author><firstname>E.</firstname><surname>Tournier</surname></author>
1774 <isbn>0-12-204230-1</isbn>
1775 <pubdate>1988</pubdate>
1777 <publishername>Academic Press</publishername>
1778 <address><city>London</city></address>