3869936c2401e200eff1811757e3afdaa91d3f91
[ginac.git] / doc / tutorial / tutorial.sgml.in
1 <!DOCTYPE Book PUBLIC "-//Davenport//DTD DocBook V3.0//EN">
2
3 <book>
4 <title>GiNaC Tutorial</title>
5 <bookinfo>
6 <subtitle>An open framework for symbolic computation within the C++ programming language</subtitle>
7 <bookbiblio>
8 <authorgroup>
9   <collab>
10     <collabname>The GiNaC Group</collabname>
11   </collab>
12   <author>
13     <firstname>Christian</firstname><surname>Bauer</surname>
14     <affiliation>
15       <address><email>Christian.Bauer@Uni-Mainz.DE</email></address>
16     </affiliation>
17   </author>
18   <author>
19     <firstname>Alexander</firstname><surname>Frink</surname>
20     <affiliation>
21       <address><email>Alexander.Frink@Uni-Mainz.DE</email></address>
22     </affiliation>
23   </author>
24   <author>
25     <firstname>Richard</firstname><othername>B.</othername><surname>Kreckel</surname>
26     <affiliation>
27       <address><email>Richard.Kreckel@Uni-Mainz.DE</email></address>
28     </affiliation>
29   </author>
30   <author>
31     <surname>Others</surname>
32     <affiliation>
33       <address><email>whoever@ThEP.Physik.Uni-Mainz.DE</email></address>
34     </affiliation>
35   </author>
36 </authorgroup>
37 </bookbiblio>
38 </bookinfo>
39
40 <preface>
41 <title>Introduction</title>
42
43 <para>The motivation behind GiNaC derives from the observation that
44 most present day computer algebra systems (CAS) are linguistically and
45 semantically impoverished.  It is an attempt to overcome the current
46 situation by extending a well established and standardized computer
47 language (C++) by some fundamental symbolic capabilities, thus
48 allowing for integrated systems that embed symbolic manipulations
49 together with more established areas of computer science (like
50 computation-intense numeric applications, graphical interfaces, etc.)
51 under one roof.</para>
52
53 <para>This tutorial is intended for the novice user who is new to GiNaC
54 but already has some background in C++ programming.  However, since a
55 hand made documentation like this one is difficult to keep in sync
56 with the development the actual documentation is inside the sources in
57 the form of comments.  That documentation may be parsed by one of the
58 many Javadoc-like documentation systems.  The generated HTML
59 documenatation is included in the distributed sources (subdir
60 <literal>doc/reference/</literal>) or can be accessed directly at URL
61 <ulink
62 url="http://wwwthep.physik.uni-mainz.de/GiNaC/reference/"><literal>http://wwwthep.physik.uni-mainz.de/GiNaC/reference/</literal></ulink>.
63 It is an invaluable resource not only for the advanced user who wishes
64 to extend the system (or chase bugs) but for everybody who wants to
65 comprehend the inner workings of GiNaC.  This little tutorial on the
66 other hand only covers the basic things that are unlikely to change in
67 the near future.
68 </para>
69
70 <sect1><title>License</title>
71
72 <para>The GiNaC framework for symbolic computation within the C++
73 programming language is Copyright (C) 1999 Johannes Gutenberg
74 Universit&auml;t Mainz, Germany.</para>
75
76 <para>This program is free software; you can redistribute it and/or
77 modify it under the terms of the GNU General Public License as
78 published by the Free Software Foundation; either version 2 of the
79 License, or (at your option) any later version.</para>
80
81 <para>This program is distributed in the hope that it will be useful, but
82 WITHOUT ANY WARRANTY; without even the implied warranty of
83 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
84 General Public License for more details.</para>
85
86 <para>You should have received a copy of the GNU General Public License
87 along with this program; see the file COPYING.  If not, write to the
88 Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
89 MA 02111-1307, USA.</para>
90
91 </preface>
92
93 <chapter>
94 <title>A Tour of GiNaC</title>
95
96 <para>This quick tour of GiNaC wants to rise your interest in in the
97 subsequent chapters by showing off a bit.  Please excuse us if it
98 leaves many open questions.</para>
99
100 <sect1><title>How to use it from within C++</title> <para>The GiNaC
101 open framework for symbolic computation within the C++ programming
102 language does not try to define a language of it's own as conventional
103 CAS do.  Instead, it extends the capabilities of C++ by symbolic
104 manipulations.  Here is how to generate and print a simple (and
105 pointless) bivariate polynomial with some large coefficients:
106 <example>
107 <title>My first GiNaC program (a bivariate polynomial)</title>
108 <programlisting>
109 #include &lt;ginac/ginac.h&gt;
110 using namespace GiNaC;
111
112 int main()
113 {
114     symbol x("x"), y("y");
115     ex poly;
116
117     for (int i=0; i<3; ++i)
118         poly += factorial(i+16)*pow(x,i)*pow(y,2-i);
119
120     cout &lt;&lt; poly &lt;&lt; endl;
121     return 0;
122 }
123 </programlisting>
124 <para>Assuming the file is called <literal>hello.cc</literal>, on 
125 our system we can compile and run it like this:</para>
126 <screen>
127 <prompt>$</prompt> c++ hello.cc -o hello -lcln -lginac
128 <prompt>$</prompt> ./hello
129 355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
130 </screen>
131 </example>
132 </para>
133
134 <para>Next, there is a more meaningful C++ program that calls a
135 function which generates Hermite polynomials in a specified free
136 variable.
137 <example>
138 <title>My second GiNaC program (Hermite polynomials)</title>
139 <programlisting>
140 #include &lt;ginac/ginac.h&gt;
141 using namespace GiNaC;
142
143 ex HermitePoly(symbol x, int deg)
144 {
145     ex HKer=exp(-pow(x,2));
146     // uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2) 
147     return normal(pow(-1,deg) * diff(HKer, x, deg) / HKer);
148 }
149
150 int main()
151 {
152     symbol z("z");
153
154     for (int i=0; i<6; ++i)
155         cout &lt;&lt; "H_" &lt;&lt; i &lt;&lt; "(z) == " &lt;&lt; HermitePoly(z,i) &lt;&lt; endl;
156
157     return 0;
158 }
159 </programlisting>
160 <para>When run, this will type out</para>
161 <screen>
162 H_0(z) == 1
163 H_1(z) == 2*z
164 H_2(z) == 4*z^2-2
165 H_3(z) == -12*z+8*z^3
166 H_4(z) == -48*z^2+16*z^4+12
167 H_5(z) == 120*z-160*z^3+32*z^5
168 </screen>
169 </example>
170 This method of generating the coefficients is of course far from
171 optimal for production purposes.</para>
172
173 <para>In order to show some more examples of what GiNaC can do we
174 will now use <literal>ginsh</literal>, a simple GiNaC interactive
175 shell that provides a convenient window into GiNaC's capabilities.
176 </para></sect1>
177
178 <sect1><title>What it can do for you</title>
179
180 <para>After invoking <literal>ginsh</literal> one can test and
181 experiment with GiNaC's features much like in other Computer Algebra
182 Systems except that it does not provide programming constructs like
183 loops or conditionals.  For a concise description of the
184 <literal>ginsh</literal> syntax we refer to its accompanied 
185 man-page.</para>
186
187 <para>It can manipulate arbitrary precision integers in a very fast
188 way.  Rational numbers are automatically converted to fractions of
189 coprime integers:
190 <screen>
191 > x=3^150;
192 369988485035126972924700782451696644186473100389722973815184405301748249
193 > y=3^149;
194 123329495011708990974900260817232214728824366796574324605061468433916083
195 > x/y;
196 3
197 > y/x;
198 1/3
199 </screen>
200 </para>
201
202 <para>All numbers occuring in GiNaC's expressions can be converted
203 into floating point numbers with the <literal>evalf</literal> method,
204 to arbitrary accuracy:
205 <screen>
206 > evalf(1/7);
207 0.14285714285714285714
208 > Digits=150;
209 150
210 > evalf(1/7);
211 0.1428571428571428571428571428571428571428571428571428571428571428571428
212 5714285714285714285714285714285714285
213 </screen>
214 </para>
215
216 <para>Exact numbers other than rationals that can be manipulated in
217 GiNaC include predefined constants like Archimedes' Pi.  They can both
218 be used in symbolic manipulations (as an exact number) as well as in
219 numeric expressions (as an inexact number):
220 <screen>
221 > a=Pi^2+x;
222 x+Pi^2
223 > evalf(a);
224 x+9.869604401089358619L0
225 > x=2;
226 2
227 > evalf(a);
228 11.869604401089358619L0
229 </screen>
230 </para>
231
232 <para>Built-in functions evaluate immediately to exact numbers if
233 this is possible.  Conversions that can be safely performed are done
234 immediately; conversions that are not generally valid are not done:
235 <screen>
236 > cos(42*Pi);
237 1
238 > cos(acos(x));
239 x
240 > acos(cos(x));
241 acos(cos(x))
242 </screen>
243 (Note that converting the last input to <literal>x</literal> would
244 allow one to conclude that <literal>42*Pi</literal> is equal to
245 <literal>0</literal>.)</para>
246
247 <para>Linear equation systems can be solved along with basic linear
248 algebra manipulations over symbolic expressions.  In C++ there is a
249 matrix class for this purpose but we can see what it can do using 
250 <literal>ginsh</literal>'s notation of double brackets to type them in:
251 <screen>
252 > lsolve(a+x*y==z,x);
253 y^(-1)*(z-a);
254 lsolve([3*x+5*y == 7, -2*x+10*y == -5], [x, y]);
255 [x==19/8,y==-1/40]
256 > M = [[ [[1, 3]], [[-3, 2]] ]];
257 [[ [[1,3]], [[-3,2]] ]]
258 > determinant(M);
259 11
260 > charpoly(M,lambda);
261 lambda^2-3*lambda+11
262 </screen>
263 </para>
264
265 <para>Multivariate polynomials and rational functions may be expanded,
266 collected and normalized (i.e. converted to a ratio of two coprime 
267 polynomials):
268 <screen>
269 > a = x^4 + 2*x^2*y^2 + 4*x^3*y + 12*x*y^3 - 3*y^4;
270 -3*y^4+x^4+12*x*y^3+2*x^2*y^2+4*x^3*y
271 > b = x^2 + 4*x*y - y^2;
272 -y^2+x^2+4*x*y
273 > expand(a*b);
274 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
275 > collect(a*b,x);
276 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 > normal(a/b);
278 3*y^2+x^2
279 </screen>
280 </para>
281
282 <para>
283 You can differentiate functions and expand them as Taylor or Laurent 
284 series (the third argument of series is the evaluation point, the 
285 fourth defines the order):
286 <screen>
287 > diff(tan(x),x);
288 tan(x)^2+1
289 > series(sin(x),x,0,4);
290 x-1/6*x^3+Order(x^4)
291 > series(1/tan(x),x,0,4);
292 x^(-1)-1/3*x+Order(x^2)
293 </screen>
294 </para>
295
296 </sect1>
297
298 </chapter>
299
300
301 <chapter>
302 <title>Installation</title>
303
304 <para>GiNaC's installation follows the spirit of most GNU software. It is
305 easily installed on your system by three steps: configuration, build,
306 installation.</para>
307
308 <sect1><title id="CLN-main">Prerequistes</title>
309
310 <para>In order to install GiNaC on your system, some prerequistes need
311 to be met.  First of all, you need to have a C++-compiler adhering to
312 the ANSI-standard <citation>ISO/IEC 14882:1998(E)</citation>.  We used
313 <literal>GCC</literal> for development so if you have a different
314 compiler you are on your own.  For the configuration to succeed you
315 need a Posix compliant shell installed in <literal>/bin/sh</literal>,
316 GNU <literal>bash</literal> is fine.  Perl is needed by the built
317 process as well, since some of the source files are automatically
318 generated by Perl scripts.  Last but not least, Bruno Haible's library
319 <literal>CLN</literal> is extensively used and needs to be installed
320 on your system.  Please get it from <ulink
321 url="ftp://ftp.santafe.edu/pub/gnu/"><literal>ftp://ftp.santafe.edu/pub/gnu/</literal></ulink>
322 or from <ulink
323 url="ftp://ftp.ilog.fr/pub/Users/haible/gnu/"><literal>ftp://ftp.ilog.fr/pub/Users/haible/gnu/</literal></ulink>
324 (it is covered by GPL) and install it prior to trying to install
325 GiNaC.  The configure script checks if it can find it and if it cannot
326 it will refuse to continue.</para></sect1>
327
328 <sect1><title>Configuration</title>
329
330 <para>To configure GiNaC means to prepare the source distribution for
331 building.  It is done via a shell script called
332 <literal>configure</literal> that is shipped with the sources.
333 (Actually, this script is by itself created with GNU Autoconf from the
334 files <literal>configure.in</literal> and
335 <literal>aclocal.m4</literal>.)  Since a configure script generated by
336 GNU Autoconf never prompts, all customization must be done either via
337 command line parameters or environment variables.  It accepts a list
338 of parameters, the complete set of which can be listed by calling it
339 with the <literal>--help</literal> option.  The most important ones
340 will be shortly described in what follows:
341 <itemizedlist>
342   <listitem>
343     <para><literal>--disable-shared</literal>: When given, this option
344     switches off the build of a shared library, i.e. a 
345     <literal>.so</literal>-file.  This may be convenient when developing
346     because it considerably speeds up compilation.</para>
347   </listitem>
348   <listitem>
349     <para><literal>--prefix=</literal><emphasis>PREFIX</emphasis>: The
350     directory where the compiled library and headers are installed. It
351     defaults to <literal>/usr/local</literal> which means that the
352     library is installed in the directory
353     <literal>/usr/local/lib</literal> and the header files in
354     <literal>/usr/local/include/GiNaC</literal> and the documentation
355     (like this one) into <literal>/usr/local/share/doc/GiNaC</literal>.</para>
356   </listitem>
357   <listitem>
358     <para><literal>--libdir=</literal><emphasis>LIBDIR</emphasis>: Use
359     this option in case you want to have the library installed in some
360     other directory than
361     <emphasis>PREFIX</emphasis><literal>/lib/</literal>.</para>
362   </listitem>
363   <listitem>
364     <para><literal>--includedir=</literal><emphasis>INCLUDEDIR</emphasis>:
365     Use this option in case you want to have the header files
366     installed in some other directory than
367     <emphasis>PREFIX</emphasis><literal>/include/ginac/</literal>. For
368     instance, if you specify
369     <literal>--includedir=/usr/include</literal> you will end up with
370     the header files sitting in the directory
371     <literal>/usr/include/ginac/</literal>. Note that the subdirectory
372     <literal>GiNaC</literal> is enforced by this process in order to
373     keep the header files separated from others.  This avoids some
374     clashes and allows for an easier deinstallation of GiNaC. This ought
375     to be considered A Good Thing (tm).</para>
376   </listitem>
377   <listitem>
378     <para><literal>--datadir=</literal><emphasis>DATADIR</emphasis>:
379     This option may be given in case you want to have the documentation 
380     installed in some other directory than
381     <emphasis>PREFIX</emphasis><literal>/share/doc/GiNaC/</literal>.
382   </listitem>
383 </itemizedlist>
384 </para>
385
386 <para>In addition, you may specify some environment variables.
387 <literal>CXX</literal> holds the path and the name of the C++ compiler
388 in case you want to override the default in your path.  (The
389 <literal>configure</literal> script searches your path for
390 <literal>c++</literal>, <literal>g++</literal>,
391 <literal>gcc</literal>, <literal>CC</literal>, <literal>cxx</literal>
392 and <literal>cc++</literal> in that order.)  It may be very useful to
393 define some compiler flags with the <literal>CXXFLAGS</literal>
394 environment variable, like optimization, debugging information and
395 warning levels.  If ommitted, it defaults to <literal>-g
396 -O2</literal>.</para>
397
398 <para>The whole process is illustrated in the following two
399 examples. (Substitute <literal>setenv VARIABLE value</literal> for
400 <literal>export VARIABLE=value</literal> if the Berkeley C shell is
401 your login shell.)
402
403 <example><title>Sample sessions of how to call the
404 configure-script</title> <para>Simple configuration for a site-wide
405 GiNaC library assuming everything is in default paths:</para>
406 <screen>
407 <prompt>$</prompt> export CXXFLAGS="-Wall -O2"
408 <prompt>$</prompt> ./configure
409 </screen>
410 <para>Configuration for a private static GiNaC library with several
411 components sitting in custom places (site-wide <literal>GCC</literal>
412 and private <literal>CLN</literal>):</para>
413 <screen>
414 <prompt>$</prompt> export CXX=/usr/local/gnu/bin/c++
415 <prompt>$</prompt> export CPPFLAGS="${CPPFLAGS} -I${HOME}/include"
416 <prompt>$</prompt> export CXXFLAGS="${CXXFLAGS} -ggdb -Wall -ansi -pedantic -O2"
417 <prompt>$</prompt> export LDFLAGS="${LDFLAGS} -L${HOME}/lib"
418 <prompt>$</prompt> ./configure --disable-shared --prefix=${HOME}
419 </screen>
420 </example>
421 </para>
422
423 </sect1>
424
425 <sect1><title>Building GiNaC</title>
426
427 <para>After proper configuration you should just build the whole
428 library by typing <literal>make</literal> at the command
429 prompt and go for a cup of coffee.</para>
430
431 <para>Just to make sure GiNaC works properly you may run a simple test
432 suite by typing <literal>make check</literal>.  This will compile some
433 sample programs, run them and compare the output to reference output.
434 Each of the checks should return a message <literal>passed</literal>
435 together with the CPU time used for that particular test.  If it does
436 not, something went wrong.  This is mostly intended to be a QA-check
437 if something was broken during the development, but not a sanity check
438 of your system.  Another intent is to allow people to fiddle around
439 with optimization.  If <literal>CLN</literal> was installed all right
440 this step is unlikely to return any errors.</para>
441
442 </sect1>
443
444 <sect1><title>Installation</title>
445
446 <para>To install GiNaC on your system, simply type <literal>make
447 install</literal>.  As described in the section about configuration
448 the files will be installed in the following directories (the
449 directories will be created if they don't already exist):
450 <itemizedlist>
451   <listitem>
452     <para><literal>libginac.a</literal> will go into
453     <emphasis>PREFIX</emphasis><literal>/lib/</literal> (or
454     <emphasis>LIBDIR</emphasis>) which defaults to
455     <literal>/usr/local/lib/</literal>.  So will
456     <literal>libginac.so</literal> if the the configure script was
457     given the option <literal>--enable-shared</literal>.  In that
458     case, the proper symlinks will be established as well (by running
459     <literal>ldconfig</literal>).</para>
460   </listitem>
461   <listitem>
462     <para>All the header files will be installed into
463     <emphasis>PREFIX</emphasis><literal>/include/ginac/</literal> (or
464     <emphasis>INCLUDEDIR</emphasis><literal>/ginac/</literal>, if
465     specified).</para>
466   </listitem>
467   <listitem>
468     <para>All documentation (HTML, Postscript and DVI) will be stuffed
469     into
470     <emphasis>PREFIX</emphasis><literal>/share/doc/GiNaC/</literal>
471     (or <emphasis>DATADIR</emphasis><literal>/doc/GiNaC/</literal>, if
472     specified).</para>
473   </listitem>
474 </itemizedlist>
475 </para>
476
477 <para>Just for the record we will list some other useful make targets:
478 <literal>make clean</literal> deletes all files generated by
479 <literal>make</literal>, i.e. all the object files.  In addition
480 <literal>make distclean</literal> removes all files generated by
481 configuration.  And finally <literal>make uninstall</literal> removes
482 the installed library and header files<footnoteref
483 linkend="warnuninstall-1">.  
484
485   <footnote id="warnuninstall-1"><para>Uninstallation does not work
486   after you have called <literal>make distclean</literal> since the
487   <literal>Makefile</literal> is itself generated by the configuration
488   from <literal>Makefile.in</literal> and hence deleted by
489   <literal>make distclean</literal>.  There are two obvious ways out
490   of this dilemma.  First, you can run the configuration again with
491   the same <emphasis>PREFIX</emphasis> thus creating a
492   <literal>Makefile</literal> with a working
493   <literal>uninstall</literal> target.  Second, you can do it by hand
494   since you now know where all the files went during
495   installation.</para></footnote>
496 </para>
497
498 </sect1>
499 </chapter>
500
501
502 <chapter>
503 <title>Basic Concepts</title>
504
505 <para>This chapter will describe the different fundamental objects
506 that can be handled with GiNaC.  But before doing so, it is worthwhile
507 introducing you to the more commonly used class of expressions,
508 representing a flexible meta-class for storing all mathematical
509 objects.</para>
510
511 <sect1><title>Expressions</title>
512
513 <para>The most common class of objects a user deals with is the
514 expression <literal>ex</literal>, representing a mathematical object
515 like a variable, number, function, sum, product, etc...  Expressions
516 may be put together to form new expressions, passed as arguments to
517 functions, and so on.  Here is a little collection of valid
518 expressions:
519 <example><title>Examples of expressions</title>
520 <programlisting>
521     ex MyEx1 = 5;                       // simple number
522     ex MyEx2 = x + 2*y;                 // polynomial in x and y
523     ex MyEx3 = (x + 1)/(x - 1);         // rational expression
524     ex MyEx4 = sin(x + 2*y) + 3*z + 41; // containing a function
525     ex MyEx5 = MyEx4 + 1;               // similar to above
526 </programlisting>
527 </example>
528 Before describing the more fundamental objects that form the building
529 blocks of expressions we'll have a quick look under the hood by
530 describing how expressions are internally managed.</para>
531
532 <sect2><title>Digression: Expressions are reference counted</title>
533
534 <para>An expression is extremely light-weight since internally it
535 works like a handle to the actual representation and really holds
536 nothing more than a pointer to some other object. What this means in
537 practice is that whenever you create two <literal>ex</literal> and set
538 the second equal to the first no copying process is involved. Instead,
539 the copying takes place as soon as you try to change the second.
540 Consider the simple sequence of code:
541 <example><title>Simple copy-on-write semantics</title>
542 <programlisting>
543 #include &lt;ginac/ginac.h&gt;
544 using namespace GiNaC;
545
546 int main()
547 {
548     symbol x("x"), y("y"), z("z");
549     ex e1, e2;
550
551     e1 = sin(x + 2*y) + 3*z + 41;
552     e2 = e1;                // e2 points to same object as e1
553     cout &lt;&lt; e2 &lt;&lt; endl;     // prints sin(x+2*y)+3*z+41
554     e2 += 1;                // e2 is copied into a new object
555     cout &lt;&lt; e2 &lt;&lt; endl;     // prints sin(x+2*y)+3*z+42
556     // ...
557 }
558 </programlisting>
559 </example>
560 The line <literal>e2 = e1;</literal> creates a second expression
561 pointing to the object held already by <literal>e1</literal>.  The
562 time involved for this operation is therefore constant, no matter how
563 large <literal>e1</literal> was.  Actual copying, however, must take
564 place in the line <literal>e2 += 1</literal> because
565 <literal>e1</literal> and <literal>e2</literal> are not handles for
566 the same object any more.  This concept is called
567 <emphasis>copy-on-write semantics</emphasis>.  It increases
568 performance considerably whenever one object occurs multiple times and
569 represents a simple garbage collection scheme because when an
570 <literal>ex</literal> runs out of scope its destructor checks whether
571 other expressions handle the object it points to too and deletes the
572 object from memory if that turns out not to be the case.  A slightly
573 less trivial example of differentiation using the chain-rule should
574 make clear how powerful this can be.  <example><title>Advanced
575 copy-on-write semantics</title>
576 <programlisting>
577 #include &lt;ginac/ginac.h&gt;
578 using namespace GiNaC;
579
580 int main()
581 {
582     symbol x("x"), y("y");
583
584     ex e1 = x + 3*y;
585     ex e2 = pow(e1, 3);
586     ex e3 = diff(sin(e2), x);   // first derivative of sin(e2) by x
587     cout &lt;&lt; e1 &lt;&lt; endl          // prints x+3*y
588          &lt;&lt; e2 &lt;&lt; endl          // prints (x+3*y)^3
589          &lt;&lt; e3 &lt;&lt; endl;         // prints 3*(x+3*y)^2*cos((x+3*y)^3)
590     // ...
591 }
592 </programlisting>
593 </example>
594 Here, <literal>e1</literal> will actually be referenced three times
595 while <literal>e2</literal> will be referenced two times.  When the
596 power of an expression is built, that expression needs not be
597 copied. Likewise, since the derivative of a power of an expression can
598 be easily expressed in terms of that expression, no copying of
599 <literal>e1</literal> is involved when <literal>e3</literal> is
600 constructed.  So, when <literal>e3</literal> is constructed it will
601 print as <literal>3*(x+3*y)^2*cos((x+3*y)^3)</literal> but the
602 argument of <literal>cos()</literal> only holds a reference to
603 <literal>e2</literal> and the factor in front is just
604 <literal>3*e1^2</literal>.
605 </para>
606
607 <para>As a user of GiNaC, you cannot see this mechanism of
608 copy-on-write semantics.  When you insert an expression into a
609 second expression, the result behaves exactly as if the contents of
610 the first expression were inserted.  But it may be useful to remember
611 that this is not what happens.  Knowing this will enable you to write
612 much more efficient code.</para>
613
614 <para>So much for expressions.  But what exactly are these expressions
615 handles of?  This will be answered in the following sections.</para>
616 </sect2>
617 </sect1>
618
619 <sect1><title>The Class Hierarchy</title>
620
621 <para>GiNaC's class hierarchy consists of several classes representing
622 mathematical objects, all of which (except for <literal>ex</literal>
623 and some helpers) are internally derived from one abstract base class
624 called <literal>basic</literal>.  You do not have to deal with objects
625 of class <literal>basic</literal>, instead you'll be dealing with
626 symbols and functions of symbols.  You'll soon learn in this chapter
627 how many of the functions on symbols are really classes.  This is
628 because simple symbolic arithmetic is not supported by languages like
629 C++ so in a certain way GiNaC has to implement its own arithmetic.</para>
630
631 <para>To give an idea about what kinds of symbolic composits may be
632 built we have a look at the most important classes in the class
633 hierarchy.  The dashed line symbolizes a "points to" or "handles"
634 relationship while the solid lines stand for "inherits from"
635 relationships.
636 <figure id="classhier-id" float="1">
637 <title>The GiNaC class hierarchy</title>
638   <graphic align="center" fileref="classhierarchy.graext" format="GRAEXT"></graphic>
639 </figure>
640 Some of the classes shown here (the ones sitting in white boxes) are
641 abstract base classes that are of no interest at all for the user.
642 They are used internally in order to avoid code duplication if
643 two or more classes derived from them share certain features.  An
644 example would be <literal>expairseq</literal>, which is a container
645 for a sequence of pairs each consisting of one expression and a number
646 (<literal>numeric</literal>).  What <emphasis>is</emphasis> visible to
647 the user are the derived classes <literal>add</literal> and
648 <literal>mul</literal>, representing sums of terms and products,
649 respectively.  We'll come back later to some more details about these
650 two classes and motivate the use of pairs in sums and products here.</para>
651
652 <sect2><title>Digression: Internal representation of products and sums</title>
653
654 <para>Although it should be completely transparent for the user of
655 GiNaC a short discussion of this topic helps understanding the sources
656 and also explains performance to a large degree.  Consider the
657 symbolic expression <literal>(a+2*b-c)*d</literal>, which could
658 naively be represented by a tree of linear containers for addition and
659 multiplication together with atomic leaves of symbols and integer
660 numbers in this fashion:
661 <figure id="repres-naive-id" float="1">
662 <title>Naive internal representation of <emphasis>d(a+2*b-c)</emphasis></title>
663   <graphic align="center" fileref="rep_naive.graext" format="GRAEXT"></graphic>
664 </figure>
665 However, doing so results in a rather deeply nested tree which will
666 quickly become rather slow to manipulate.  If we represent the sum
667 instead as a sequence of terms, each having a purely numeric
668 coefficient, the tree becomes much more flat.
669 <figure id="repres-pair-id" float="1">
670 <title>Better internal representation of <emphasis>d(a+2*b-c)</emphasis></title>
671   <graphic align="center" fileref="rep_pair.graext" format="GRAEXT"></graphic>
672 </figure>
673 The number <literal>1</literal> under the symbol <literal>d</literal>
674 is a hint that multiplication objects can be treated similarly where
675 the coeffiecients are interpreted as <emphasis>exponents</emphasis>
676 now.  Addition of sums of terms or multiplication of products with
677 numerical exponents can be made very efficient with a
678 pair-representation.  Internally, this handling is done by most CAS in
679 this way.  It typically speeds up manipulations by an order of
680 magnitude.  Now it should be clear, why both classes
681 <literal>add</literal> and <literal>mul</literal> are derived from the
682 same abstract class: the representation is the same, only the
683 interpretation differs.  </para>
684
685 </sect1>
686
687 <sect1><title>Symbols</title>
688
689 <para>Symbols are for symbolic manipulation what atoms are for
690 chemistry.  You can declare objects of class <literal>symbol</literal>
691 as any other object simply by saying <literal>symbol x,y;</literal>.
692 There is, however, a catch in here having to do with the fact that C++
693 is a compiled language.  The information about the symbol's name is
694 thrown away by the compiler but at a later stage you may want to print
695 expressions holding your symbols.  In order to avoid confusion GiNaC's
696 symbols are able to know their own name.  This is accomplished by
697 declaring its name for output at construction time in the fashion
698 <literal>symbol x("x");</literal>.  If you declare a symbol using the
699 default constructor (i.e. without string-argument) the system will
700 deal out a unique name.  That name may not be suitable for printing
701 but for internal routines when no output is desired it is often
702 enough.  We'll come across examples of such symbols later in this
703 tutorial.  </para>
704
705 <para>Although symbols can be assigned expressions for internal
706 reasons, you should not do it (and we are not going to tell you how it
707 is done).  If you want to replace a symbol with something else in an
708 expression, you can use the expression's <literal>.subs()</literal>
709 method.</para>
710
711 </sect1>
712
713 <sect1><title>Numbers</title>
714
715 <para>For storing numerical things, GiNaC uses Bruno Haible's library
716 <literal>CLN</literal>.  The classes therein serve as foundation
717 classes for GiNaC.  <literal>CLN</literal> stands for Class Library
718 for Numbers or alternatively for Common Lisp Numbers.  In order to
719 find out more about <literal>CLN</literal>'s internals the reader is
720 refered to the documentation of that library.  Suffice to say that it
721 is by itself build on top of another library, the GNU Multiple
722 Precision library <literal>GMP</literal>, which is a extremely fast
723 library for arbitrary long integers and rationals as well as arbitrary
724 precision floating point numbers.  It is very commonly used by several
725 popular cryptographic applications.  <literal>CLN</literal> extends
726 <literal>GMP</literal> by several useful things: First, it introduces
727 the complex number field over either reals (i.e. floating point
728 numbers with arbitrary precision) or rationals.  Second, it
729 automatically converts rationals to integers if the denominator is
730 unity and complex numbers to real numbers if the imaginary part
731 vanishes and also correctly treats algebraic functions.  Third it
732 provides good implementations of state-of-the-art algorithms for all
733 trigonometric and hyperbolic functions as well as for calculation of
734 some useful constants.</para>
735
736 <para>The user can construct an object of class
737 <literal>numeric</literal> in several ways.  The following example
738 shows the four most important constructors: construction from
739 C-integer, construction of fractions from two integers, construction
740 from C-float and construction from a string.
741 <example><title>Construction of numbers</title>
742 <programlisting>
743 #include &lt;ginac/ginac.h&gt;
744 using namespace GiNaC;
745
746 int main()
747 {
748     numeric two(2);                     // exact integer 2
749     numeric r(2,3);                     // exact fraction 2/3
750     numeric e(2.71828);                 // floating point number
751     numeric p("3.1415926535897932385"); // floating point number
752
753     cout &lt;&lt; two*p &lt;&lt; endl;  // floating point 6.283...
754     // ...
755 }
756 </programlisting>
757 </example>
758 Note that all those constructors are <emphasis>explicit</emphasis>
759 which means you are not allowed to write <literal>numeric
760 two=2;</literal>.  This is because the basic objects to be handled by
761 GiNaC are the expressions and we want to keep things simple and wish
762 objects like <literal>pow(x,2)</literal> to be handled the same way
763 as <literal>pow(x,a)</literal>, which means that we need to allow a
764 general expression as base and exponent.  Therefore there is an
765 implicit construction from a C-integer directly to an expression
766 handling a numeric in the first example.  This design really becomes
767 convenient when one declares own functions having more than one
768 parameter but it forbids using implicit constructors because that
769 would lead to ambiguities.
770 </para>
771
772 <para>We have seen now the distinction between exact numbers and
773 floating point numbers.  Clearly, the user should never have to worry
774 about dynamically created exact numbers, since their "exactness"
775 always determines how they ought to be handled.  The situation is
776 different for floating point numbers.  Their accuracy is handled by
777 one <emphasis>global</emphasis> variable, called
778 <literal>Digits</literal>.  (For those readers who know about Maple:
779 it behaves very much like Maple's <literal>Digits</literal>).  All
780 objects of class numeric that are constructed from then on will be
781 stored with a precision matching that number of decimal digits:
782 <example><title>Controlling the precision of floating point numbers</title>
783 <programlisting> 
784 #include &lt;ginac/ginac.h&gt;
785 using namespace GiNaC;
786
787 void foo()
788 {
789     numeric three(3.0), one(1.0);
790     numeric x = one/three;
791
792     cout &lt;&lt; "in " &lt;&lt; Digits &lt;&lt; " digits:" &lt;&lt; endl;
793     cout &lt;&lt; x &lt;&lt; endl;
794     cout &lt;&lt; Pi.evalf() &lt;&lt; endl;
795 }
796
797 int main()
798 {
799     foo();
800     Digits = 60;
801     foo();
802     return 0;
803 }
804 </programlisting>
805 </example>
806 The above example prints the following output to screen:
807 <screen>
808 in 17 digits:
809 0.333333333333333333
810 3.14159265358979324
811 in 60 digits:
812 0.333333333333333333333333333333333333333333333333333333333333333333
813 3.14159265358979323846264338327950288419716939937510582097494459231
814 </screen>
815 </para>
816
817 <sect2><title>Tests on numbers</title>
818
819 <para>Once you have declared some numbers, assigned them to
820 expressions and done some arithmetic with them it is frequently
821 desired to retrieve some kind of information from them like asking
822 whether that number is integer, rational, real or complex.  For those
823 cases GiNaC provides several useful methods.  (Internally, they fall
824 back to invocations of certain CLN functions.)</para>
825
826 <para>As an example, let's construct some rational number, multiply it
827 with some multiple of its denominator and check what comes out:
828 <example><title>Sample test on objects of type numeric</title>
829 <programlisting>
830 #include &lt;ginac/ginac.h&gt;
831 using namespace GiNaC;
832
833 int main()
834 {
835     numeric twentyone(21);
836     numeric ten(10);
837     numeric answer(21,5);
838
839     cout &lt;&lt; answer.is_integer() &lt;&lt; endl;  // false, it's 21/5
840     answer *= ten;
841     cout &lt;&lt; answer.is_integer() &lt;&lt; endl;  // true, it's 42 now!
842     // ...
843 }
844 </programlisting>
845 </example>
846 Note that the variable <literal>answer</literal> is constructed here
847 as an integer but in an intermediate step it holds a rational number
848 represented as integer numerator and denominator.  When multiplied by
849 10, the denominator becomes unity and the result is automatically
850 converted to a pure integer again.  Internally, the underlying
851 <literal>CLN</literal> is responsible for this behaviour and we refer
852 the reader to <literal>CLN</literal>'s documentation.  Suffice to say
853 that the same behaviour applies to complex numbers as well as return
854 values of certain functions.  Complex numbers are automatically
855 converted to real numbers if the imaginary part becomes zero.  The
856 full set of tests that can be applied is listed in the following
857 table.
858 <informaltable colsep="0" frame="topbot" pgwide="1">
859 <tgroup cols="2">
860 <colspec colnum="1" colwidth="1*">
861 <colspec colnum="2" colwidth="2*">
862 <thead>
863   <row>
864     <entry>Method</entry>
865     <entry>Returns true if...</entry>
866   </row>
867 </thead>
868 <tbody>
869   <row>
870     <entry><literal>.is_zero()</literal></entry>
871     <entry>object is equal to zero</entry>
872   </row>
873   <row>
874     <entry><literal>.is_positive()</literal></entry>
875     <entry>object is not complex and greater than 0</entry>
876   </row>
877   <row>
878     <entry><literal>.is_integer()</literal></entry>
879     <entry>object is a (non-complex) integer</entry>
880   </row>
881   <row>
882     <entry><literal>.is_pos_integer()</literal></entry>
883     <entry>object is an integer and greater than 0</entry>
884   </row>
885   <row>
886     <entry><literal>.is_nonneg_integer()</literal></entry>
887     <entry>object is an integer and greater equal 0</entry>
888   </row>
889   <row>
890     <entry><literal>.is_even()</literal></entry>
891     <entry>object is an even integer</entry>
892   </row>
893   <row>
894     <entry><literal>.is_odd()</literal></entry>
895     <entry>object is an odd integer</entry>
896   </row>
897   <row>
898     <entry><literal>.is_prime()</literal></entry>
899     <entry>object is a prime integer (probabilistic primality test)</entry>
900   </row>
901   <row>
902     <entry><literal>.is_rational()</literal></entry>
903     <entry>object is an exact rational number (integers are rational, too, as are complex extensions like <literal>2/3+7/2*I</literal>)</entry>
904   </row>
905   <row>
906     <entry><literal>.is_real()</literal></entry>
907     <entry>object is a real integer, rational or float (i.e. is not complex)</entry>
908   </row>
909 </tbody>
910 </tgroup>
911 </informaltable>
912 </para>
913
914 </sect2>
915
916 </sect1>
917
918
919 <sect1><title>Constants</title>
920
921 <para>Constants behave pretty much like symbols except that that they return
922 some specific number when the method <literal>.evalf()</literal> is called.
923 </para>
924
925 <para>The predefined known constants are:
926 <informaltable colsep="0" frame="topbot" pgwide="1">
927 <tgroup cols="3">
928 <colspec colnum="1" colwidth="1*">
929 <colspec colnum="2" colwidth="2*">
930 <colspec colnum="3" colwidth="4*">
931 <thead>
932   <row>
933     <entry>Name</entry>
934     <entry>Common Name</entry>
935     <entry>Numerical Value (35 digits)</entry>
936   </row>
937 </thead>
938 <tbody>
939   <row>
940     <entry><literal>Pi</literal></entry>
941     <entry>Archimedes' constant</entry>
942     <entry>3.14159265358979323846264338327950288</entry>
943   </row><row>
944     <entry><literal>Catalan</literal></entry>
945     <entry>Catalan's constant</entry>
946     <entry>0.91596559417721901505460351493238411</entry>
947   </row><row>
948     <entry><literal>EulerGamma</literal></entry>
949     <entry>Euler's (or Euler-Mascheroni) constant</entry>
950     <entry>0.57721566490153286060651209008240243</entry>
951   </row>
952 </tbody>
953 </tgroup>
954 </informaltable>
955 </para>
956
957 </sect1>
958
959 <sect1><title>Fundamental operations: The <literal>power</literal>, <literal>add</literal> and <literal>mul</literal> classes</title>
960
961 <para>Simple polynomial expressions are written down in GiNaC pretty
962 much like in other CAS.  The necessary operators <literal>+</literal>,
963 <literal>-</literal>, <literal>*</literal> and <literal>/</literal>
964 have been overloaded to achieve this goal.  When you run the following
965 program, the constructor for an object of type <literal>mul</literal>
966 is automatically called to hold the product of <literal>a</literal>
967 and <literal>b</literal> and then the constructor for an object of
968 type <literal>add</literal> is called to hold the sum of that
969 <literal>mul</literal> object and the number one:
970 <example><title>Construction of <literal>add</literal> and <literal>mul</literal> objects</title>
971 <programlisting>
972 #include &lt;ginac/ginac.h&gt;
973 using namespace GiNaC;
974
975 int main()
976 {
977     symbol a("a"), b("b");
978     ex MyTerm = 1+a*b;
979     // ...
980 }
981 </programlisting>
982 </example></para>
983
984 <para>For exponentiation, you have already seen the somewhat clumsy
985 (though C-ish) statement <literal>pow(x,2);</literal> to represent
986 <literal>x</literal> squared.  This direct construction is necessary
987 since we cannot safely overload the constructor <literal>^</literal>
988 in <literal>C++</literal> to construct a <literal>power</literal>
989 object.  If we did, it would have several counterintuitive effects:
990 <itemizedlist>
991   <listitem>
992     <para>Due to <literal>C</literal>'s operator precedence,
993     <literal>2*x^2</literal> would be parsed as <literal>(2*x)^2</literal>.
994   </listitem>
995   <listitem>
996     <para>Due to the binding of the operator <literal>^</literal>, 
997     <literal>x^a^b</literal> would result in <literal>(x^a)^b</literal>. 
998     This would be confusing since most (though not all) other CAS 
999     interpret this as <literal>x^(a^b)</literal>.
1000   </listitem>
1001   <listitem>
1002     <para>Also, expressions involving integer exponents are very 
1003     frequently used, which makes it even more dangerous to overload 
1004     <literal>^</literal> since it is then hard to distinguish between the
1005     semantics as exponentiation and the one for exclusive or.  (It would
1006     be embarassing to return <literal>1</literal> where one has requested 
1007     <literal>2^3</literal>.)</para>
1008   </listitem>
1009 </itemizedlist>
1010 All effects are contrary to mathematical notation and differ from the
1011 way most other CAS handle exponentiation, therefore overloading
1012 <literal>^</literal> is ruled out for GiNaC's C++ part.  The situation
1013 is different in <literal>ginsh</literal>, there the
1014 exponentiation-<literal>^</literal> exists.  (Also note, that the
1015 other frequently used exponentiation operator <literal>**</literal>
1016 does not exist at all in <literal>C++</literal>).</para>
1017
1018 <para>To be somewhat more precise, objects of the three classes
1019 described here, are all containers for other expressions.  An object
1020 of class <literal>power</literal> is best viewed as a container with
1021 two slots, one for the basis, one for the exponent.  All valid GiNaC
1022 expressions can be inserted.  However, basic transformations like
1023 simplifying <literal>pow(pow(x,2),3)</literal> to
1024 <literal>x^6</literal> automatically are only performed when
1025 this is mathematically possible.  If we replace the outer exponent
1026 three in the example by some symbols <literal>a</literal>, the
1027 simplification is not safe and will not be performed, since
1028 <literal>a</literal> might be <literal>1/2</literal> and
1029 <literal>x</literal> negative.</para>
1030
1031 <para>Objects of type <literal>add</literal> and
1032 <literal>mul</literal> are containers with an arbitrary number of
1033 slots for expressions to be inserted.  Again, simple and safe
1034 simplifications are carried out like transforming
1035 <literal>3*x+4-x</literal> to <literal>2*x+4</literal>.</para>
1036
1037 <para>The general rule is that when you construct such objects, GiNaC
1038 automatically creates them in canonical form, which might differ from
1039 the form you typed in your program.  This allows for rapid comparison
1040 of expressions, since after all <literal>a-a</literal> is simply zero.
1041 Note, that the canonical form is not necessarily lexicographical
1042 ordering or in any way easily guessable.  It is only guaranteed that
1043 constructing the same expression twice, either implicitly or
1044 explicitly, results in the same canonical form.</para>
1045
1046 </sect1>
1047
1048 <sect1><title>Built-in Functions</title>
1049
1050 <para>This section is not here yet</para>
1051
1052
1053
1054 </sect1>
1055
1056 </chapter>
1057
1058
1059 <chapter>
1060 <title>Important Algorithms</title>
1061
1062 <para>In this chapter the most important algorithms provided by GiNaC
1063 will be described.  Some of them are implemented as functions on
1064 expressions, others are implemented as methods provided by expression
1065 objects.  If they are methods, there exists a wrapper function around
1066 it, so you can alternatively call it in a functional way as shown in
1067 the simple example:
1068 <example><title>Methods vs. wrapper functions</title>
1069 <programlisting>
1070 #include &lt;ginac/ginac.h&gt;
1071 using namespace GiNaC;
1072
1073 int main()
1074 {
1075     ex x = numeric(1.0);
1076     
1077     cout &lt;&lt; "As method:   " &lt;&lt; sin(x).evalf() &lt;&lt; endl;
1078     cout &lt;&lt; "As function: " &lt;&lt; evalf(sin(x)) &lt;&lt; endl;
1079     // ...
1080 }
1081 </programlisting>
1082 </example>
1083 The general rule is that wherever methods accept one or more
1084 parameters (<emphasis>arg1</emphasis>, <emphasis>arg2</emphasis>, ...)
1085 the order of arguments the function wrapper accepts is the same but
1086 preceded by the object to act on (<emphasis>object</emphasis>,
1087 <emphasis>arg1</emphasis>, <emphasis>arg2</emphasis>, ...).  This
1088 approach is the most natural one in an OO model but it may lead to
1089 confusion for MapleV users because where they would type
1090 <literal>A:=x+1; subs(x=2,A);</literal> GiNaC would require
1091 <literal>A=x+1; subs(A,x==2);</literal> (after proper declaration of A
1092 and x).  On the other hand, since MapleV returns 3 on
1093 <literal>A:=x^2+3; coeff(A,x,0);</literal> (GiNaC:
1094 <literal>A=pow(x,2)+3; coeff(A,x,0);</literal>) it is clear that
1095 MapleV is not trying to be consistent here.  Also, users of MuPAD will
1096 in most cases feel more comfortable with GiNaC's convention.  All
1097 function wrappers are always implemented as simple inline functions
1098 which just call the corresponding method and are only provided for
1099 users uncomfortable with OO who are dead set to avoid method
1100 invocations.  Generally, a chain of function wrappers is much harder
1101 to read than a chain of methods and should therefore be avoided if
1102 possible.  On the other hand, not everything in GiNaC is a method on
1103 class <literal>ex</literal> and sometimes calling a function cannot be
1104 avoided.
1105 </para>
1106
1107 <sect1><title>Polynomial Expansion</title>
1108
1109 <para>A polynomial in one or more variables has many equivalent
1110 representations.  Some useful ones serve a specific purpose.  Consider
1111 for example the trivariate polynomial <literal>4*x*y + x*z + 20*y^2 +
1112 21*y*z + 4*z^2</literal>.  It is equivalent to the factorized
1113 polynomial <literal>(x + 5*y + 4*z)*(4*y + z)</literal>.  Other
1114 representations are the recursive ones where one collects for
1115 exponents in one of the three variable.  Since the factors are
1116 themselves polynomials in the remaining two variables the procedure
1117 can be repeated.  In our expample, two possibilies would be
1118 <literal>(4*y + z)*x + 20*y^2 + 21*y*z + 4*z^2</literal> and
1119 <literal>20*y^2 + (21*z + 4*x)*y + 4*z^2 + x*z</literal>.
1120 </para>
1121
1122 <para>To bring an expression into expanded form, its method
1123 <function>.expand()</function> may be called.  In our example above,
1124 this corresponds to <literal>4*x*y + x*z + 20*y^2 + 21*y*z +
1125 4*z^2</literal>.  Again, since the canonical form in GiNaC is not
1126 easily guessable you should be prepared to see different orderings of
1127 terms in such sums!</para>
1128
1129 </sect1>
1130
1131 <sect1><title>Collecting expressions</title>
1132
1133 <para>Another useful representation of multivariate polynomials is as
1134 a univariate polynomial in one of the variables with the coefficients
1135 being polynomials in the remaining variables.  The method
1136 <literal>collect()</literal> accomplishes this task:
1137 <funcsynopsis>
1138   <funcsynopsisinfo>#include &lt;ginac/ginac.h></funcsynopsisinfo>
1139   <funcdef>ex <function>ex::collect</function></funcdef>
1140   <paramdef>symbol const & <parameter>s</parameter></paramdef>
1141 </funcsynopsis>
1142 Note that the original polynomial needs to be in expanded form in
1143 order to be able to find the coefficients properly.  The range of
1144 occuring coefficients can be checked using the two methods
1145 <funcsynopsis>
1146   <funcsynopsisinfo>#include &lt;ginac/ginac.h></funcsynopsisinfo>
1147   <funcdef>int <function>ex::degree</function></funcdef>
1148   <paramdef>symbol const & <parameter>s</parameter></paramdef>
1149 </funcsynopsis>
1150 <funcsynopsis>
1151   <funcdef>int <function>ex::ldegree</function></funcdef>
1152   <paramdef>symbol const & <parameter>s</parameter></paramdef>
1153 </funcsynopsis>
1154 where <literal>degree()</literal> returns the highest coefficient and
1155 <literal>ldegree()</literal> the lowest one.  These two methods work
1156 also reliably on non-expanded input polynomials.  This is illustrated
1157 in the following example: 
1158
1159 <example><title>Collecting expressions in multivariate polynomials</title>
1160 <programlisting>
1161 #include &lt;ginac/ginac.h&gt;
1162 using namespace GiNaC;
1163
1164 int main()
1165 {
1166     symbol x("x"), y("y");
1167     ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
1168                  - pow(x+y,2) + 2*pow(y+2,2) - 8;
1169     ex Poly = PolyInp.expand();
1170     
1171     for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) {
1172         cout &lt;&lt; "The x^" &lt;&lt; i &lt;&lt; "-coefficient is "
1173              &lt;&lt; Poly.coeff(x,i) &lt;&lt; endl;
1174     }
1175     cout &lt;&lt; "As polynomial in y: " 
1176          &lt;&lt; Poly.collect(y) &lt;&lt; endl;
1177     // ...
1178 }
1179 </programlisting>
1180 </example>
1181 When run, it returns an output in the following fashion:
1182 <screen>
1183 The x^0-coefficient is y^2+11*y
1184 The x^1-coefficient is 5*y^2-2*y
1185 The x^2-coefficient is -1
1186 The x^3-coefficient is 4*y
1187 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
1188 </screen>
1189 As always, the exact output may vary between different versions of
1190 GiNaC or even from run to run since the internal canonical ordering is
1191 not within the user's sphere of influence.</para>
1192
1193 </sect1>
1194
1195 <sect1><title>Polynomial Arithmetic</title>
1196
1197 <sect2><title>GCD and LCM</title>
1198
1199 <para>The functions for polynomial greatest common divisor and least common
1200 multiple have the synopsis:
1201 <funcsynopsis>
1202   <funcsynopsisinfo>#include &lt;GiNaC/normal.h></funcsynopsisinfo>
1203   <funcdef>ex <function>gcd</function></funcdef>
1204   <paramdef>const ex *<parameter>a</parameter>, const ex *<parameter>b</parameter></paramdef>
1205 </funcsynopsis>
1206 <funcsynopsis>
1207   <funcdef>ex <function>lcm</function></funcdef>
1208   <paramdef>const ex *<parameter>a</parameter>, const ex *<parameter>b</parameter></paramdef>
1209 </funcsynopsis></para>
1210
1211 <para>The functions <function>gcd()</function> and <function
1212 id="lcm-main">lcm()</function> accepts two expressions
1213 <literal>a</literal> and <literal>b</literal> as arguments and return
1214 a new expression, their greatest common divisor or least common
1215 multiple, respectively.  If the polynomials <literal>a</literal> and
1216 <literal>b</literal> are coprime <function>gcd(a,b)</function> returns 1
1217 and <function>lcm(a,b)</function> returns the product of
1218 <literal>a</literal> and <literal>b</literal>.
1219 <example><title>Polynomal GCD/LCM</title>
1220 <programlisting>
1221 #include &lt;ginac/ginac.h&gt;
1222 using namespace GiNaC;
1223
1224 int main()
1225 {
1226     symbol x("x"), y("y"), z("z");
1227     ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
1228     ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
1229
1230     ex P_gcd = gcd(P_a, P_b);
1231     // x + 5*y + 4*z
1232     ex P_lcm = lcm(P_a, P_b);
1233     // 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
1234     // ...
1235 }
1236 </programlisting>
1237 </example>
1238 </para>
1239
1240 </sect2>
1241
1242 <sect2><title>The <function>normal</function> method</title>
1243
1244 <para>While in common symbolic code <function>gcd()</function> and
1245 <function>lcm()</function> are not too heavily used, simplification
1246 occurs frequently.  Therefore <function>.normal()</function>, which
1247 provides some basic form of simplification, has become a method of
1248 class <literal>ex</literal>, just like <literal>.expand()</literal>.
1249 It converts a rational function into an equivalent rational function
1250 where numererator and denominator are coprime.  This means, it finds
1251 the GCD of numerator and denominator and cancels it.  If it encounters
1252 some object which does not belong to the domain of rationals (a
1253 function for instance), that object is replaced by a temporary symbol.
1254 This means that both expressions <literal>t1</literal> and
1255 <literal>t2</literal> are indeed simplified in this little program:
1256 <example><title>Cancellation of polynomial GCD (with obstacles)</title>
1257 <programlisting>
1258 #include &lt;ginac/ginac.h&gt;
1259 using namespace GiNaC;
1260
1261 int main()
1262 {
1263     symbol x("x");
1264     ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
1265     ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
1266     cout &lt;&lt; "t1 is " &lt;&lt; t1.normal() &lt;&lt; endl;
1267     cout &lt;&lt; "t2 is " &lt;&lt; t2.normal() &lt;&lt; endl;
1268     // ...
1269 }
1270 </programlisting>
1271 </example>
1272
1273 Of course this works for multivariate polynomials too, so the ratio of
1274 the sample-polynomials from the section about GCD and LCM above would
1275 be normalized to <literal>P_a/P_b</literal> =
1276 <literal>(4*y+z)/(y+3*z)</literal>.</para>
1277
1278 </sect2>
1279
1280 </sect1>
1281
1282 <sect1><title>Symbolic Differentiation</title>
1283
1284 <para>GiNaC's objects know how to differentiate themselves.  Thus, a
1285 polynomial (class <literal>add</literal>) knows that its derivative is
1286 the sum of the derivatives of all the monomials:
1287 <example><title>Simple polynomial differentiation</title>
1288 <programlisting>
1289 #include &lt;ginac/ginac.h&gt;
1290 using namespace GiNaC;
1291
1292 int main()
1293 {
1294     symbol x("x"), y("y"), z("z");
1295     ex P = pow(x, 5) + pow(x, 2) + y;
1296
1297     cout &lt;&lt; P.diff(x,2) &lt;&lt; endl;  // 20*x^3 + 2
1298     cout &lt;&lt; P.diff(y) &lt;&lt; endl;    // 1
1299     cout &lt;&lt; P.diff(z) &lt;&lt; endl;    // 0
1300     // ...
1301 }
1302 </programlisting>
1303 </example>
1304 If a second integer parameter <emphasis>n</emphasis> is given the
1305 <function>diff</function> method returns the <emphasis>n</emphasis>th
1306 derivative.</para>
1307
1308 <para>If <emphasis>every</emphasis> object and every function is told
1309 what its derivative is, all derivatives of composed objects can be
1310 calculated using the chain rule and the product rule.  Consider, for
1311 instance the expression <literal>1/cosh(x)</literal>.  Since the
1312 derivative of <literal>cosh(x)</literal> is <literal>sinh(x)</literal>
1313 and the derivative of <literal>pow(x,-1)</literal> is
1314 <literal>-pow(x,-2)</literal> GiNaC can readily compute the
1315 composition.  It turns out that the composition is the generating
1316 function for Euler Numbers, i.e. the so called
1317 <emphasis>n</emphasis>th Euler number is the coefficient of
1318 <literal>x^n/n!</literal> in the expansion of
1319 <literal>1/cosh(x)</literal>.  We may use this identity to code a
1320 function that generates Euler numbers in just three lines:
1321 <example><title>Differentiation with nontrivial functions: Euler numbers</title>
1322 <programlisting>
1323 #include &lt;ginac/ginac.h&gt;
1324 using namespace GiNaC;
1325
1326 ex EulerNumber(unsigned n)
1327 {
1328     symbol x;
1329     ex generator = pow(cosh(x),-1);
1330     return generator.diff(x,n).subs(x==0);
1331 }
1332
1333 int main()
1334 {
1335     for (unsigned i=0; i&lt;11; i+=2)
1336         cout &lt;&lt; EulerNumber(i) &lt;&lt; endl;
1337     return 0;
1338 }
1339 </programlisting>
1340 </example>
1341 When you run it, it produces the sequence <literal>1</literal>,
1342 <literal>-1</literal>, <literal>5</literal>, <literal>-61</literal>,
1343 <literal>1385</literal>, <literal>-50521</literal>.  We increment the
1344 loop variable <literal>i</literal> by two since all odd Euler numbers
1345 vanish anyways.</para>
1346
1347 </sect1>
1348
1349 <sect1><title>Series Expansion</title>
1350
1351 <para>Expressions know how to expand themselves as a Taylor series or
1352 (more generally) a Laurent series.  As in most conventional Computer
1353 Algebra Systems no distinction is made between those two.  There is a
1354 class of its own for storing such series as well as a class for
1355 storing the order of the series.  A sample program could read:
1356 <example><title>Series expansion</title>
1357 <programlisting>
1358 #include &lt;ginac/ginac.h&gt;
1359 using namespace GiNaC;
1360
1361 int main()
1362 {
1363     symbol x("x");
1364     numeric point(0);
1365     ex MyExpr1 = sin(x);
1366     ex MyExpr2 = 1/(x - pow(x, 2) - pow(x, 3));
1367     ex MyTailor, MySeries;
1368     
1369     MyTailor = MyExpr1.series(x, point, 5);
1370     cout &lt;&lt; MyExpr1 &lt;&lt; " == " &lt;&lt; MyTailor
1371          &lt;&lt; " for small " &lt;&lt; x &lt;&lt; endl;
1372     MySeries = MyExpr2.series(x, point, 7);
1373     cout &lt;&lt; MyExpr2 &lt;&lt; " == " &lt;&lt; MySeries
1374          &lt;&lt; " for small " &lt;&lt; x &lt;&lt; endl;
1375     // ...
1376 }
1377 </programlisting>
1378 </example>
1379 </para>
1380
1381 <para>As an instructive application, let us calculate the numerical
1382 value of Archimedes' constant (for which there already exists the
1383 built-in constant <literal>Pi</literal>) using M&eacute;chain's
1384 wonderful formula <literal>Pi==16*atan(1/5)-4*atan(1/239)</literal>.
1385 We may expand the arcus tangent around <literal>0</literal> and insert
1386 the fractions <literal>1/5</literal> and <literal>1/239</literal>.
1387 But, as we have seen, a series in GiNaC carries an order term with it.
1388 The function <literal>series_to_poly</literal> may be used to strip 
1389 this off:
1390 <example><title>Series expansion using M&eacute;chain's formula for 
1391 <literal>Pi</literal></title>
1392 <programlisting>
1393 #include &lt;ginac/ginac.h&gt;
1394 using namespace GiNaC;
1395
1396 ex mechain_pi(int degr)
1397 {
1398     symbol x("x");
1399     ex pi_expansion = series_to_poly(atan(x).series(x,0,degr));
1400     ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
1401                    -4*pi_expansion.subs(x==numeric(1,239));
1402     return pi_approx;
1403 }
1404
1405 int main()
1406 {
1407     ex pi_frac;
1408     for (int i=2; i&lt;12; i+=2) {
1409         pi_frac = mechain_pi(i);
1410         cout &lt;&lt; i &lt;&lt; ":\t" &lt;&lt; pi_frac &lt;&lt; endl
1411              &lt;&lt; "\t" &lt;&lt; pi_frac.evalf() &lt;&lt; endl;
1412     }
1413     return 0;
1414 }
1415 </programlisting>
1416 <para>When you run this program, it will type out:</para>
1417 <screen>
1418 2:      3804/1195
1419         3.1832635983263598326
1420 4:      5359397032/1706489875
1421         3.1405970293260603143
1422 6:      38279241713339684/12184551018734375
1423         3.141621029325034425
1424 8:      76528487109180192540976/24359780855939418203125
1425         3.141591772182177295
1426 10:     327853873402258685803048818236/104359128170408663038552734375
1427         3.1415926824043995174
1428 </screen>
1429 </example>
1430 </para>
1431
1432 </sect1>
1433
1434 </chapter>
1435
1436
1437 <chapter>
1438 <title>Extending GiNaC</title>
1439
1440 <para>Longish chapter follows here.</para>
1441
1442 </chapter>
1443
1444
1445 <chapter>
1446 <title>A Comparison with other CAS</title>
1447
1448 <para>This chapter will give you some information on how GiNaC
1449 compares to other, traditional Computer Algebra Systems, like
1450 <literal>Maple</literal>, <literal>Mathematica</literal> or
1451 <literal>Reduce</literal>, where it has advantages and disadvantages
1452 over these systems.</para>
1453
1454 <sect1><title>Advantages</title>
1455
1456 <para>GiNaC has several advantages over traditional Computer
1457 Algebra Systems, like 
1458
1459 <itemizedlist>
1460   <listitem>
1461     <para>familiar language: all common CAS implement their own
1462     proprietary grammar which you have to learn first (and maybe learn
1463     again when your vendor chooses to "enhance" it).  With GiNaC you
1464     can write your program in common <literal>C++</literal>, which is
1465     standardized.</para>
1466   </listitem>
1467   <listitem>
1468     <para>structured data types: you can build up structured data
1469     types using <literal>struct</literal>s or <literal>class</literal>es
1470     together with STL features instead of using unnamed lists of lists
1471     of lists.</para>
1472   </listitem>
1473   <listitem>
1474     <para>strongly typed: in CAS, you usually have only one kind of
1475     variables which can hold contents of an arbitrary type.  This
1476     4GL like feature is nice for novice programmers, but dangerous.
1477     </para>
1478   </listitem>
1479   <listitem>
1480     <para>development tools: powerful development tools exist for
1481     <literal>C++</literal>, like fancy editors (e.g. with automatic
1482     indentation and syntax highlighting), debuggers, visualization
1483     tools, documentation tools...</para>
1484   </listitem>
1485   <listitem>
1486     <para>modularization: <literal>C++</literal> programs can 
1487     easily be split into modules by separating interface and
1488     implementation.</para>
1489   </listitem>
1490   <listitem>
1491     <para>price: GiNaC is distributed under the GNU Public License
1492     which means that it is free and available with source code.  And
1493     there are excellent <literal>C++</literal>-compilers for free, too.
1494     </para>
1495   </listitem>
1496   <listitem>
1497     <para>extendable: you can add your own classes to GiNaC, thus
1498     extending it on a very low level.  Compare this to a traditional
1499     CAS that you can usually only extend on a high level by writing in
1500     the language defined by the parser.  In particular, it turns out
1501     to be almost impossible to fix bugs in a traditional system.
1502   </listitem>
1503   <listitem>
1504     <para>seemless integration: it is somewhere between difficult
1505     and impossible to call CAS functions from within a program 
1506     written in <literal>C++</literal> or any other programming 
1507     language and vice versa.  With GiNaC, your symbolic routines
1508     are part of your program.  You can easily call third party
1509     libraries, e.g. for numerical evaluation or graphical 
1510     interaction.  All other approaches are much more cumbersome: they
1511     range from simply ignoring the problem
1512     (i.e. <literal>Maple</literal>) to providing a
1513     method for "embedding" the system
1514     (i.e. <literal>Yacas</literal>).</para>
1515   </listitem>
1516   <listitem>
1517     <para>efficiency: often large parts of a program do not need
1518     symbolic calculations at all.  Why use large integers for loop
1519     variables or arbitrary precision arithmetics where double
1520     accuracy is sufficient?  For pure symbolic applications,
1521     GiNaC is comparable in speed with other CAS.
1522   </listitem>
1523 </itemizedlist>
1524 </para>
1525
1526 <sect1><title>Disadvantages</title>
1527
1528 <para>Of course it also has some disadvantages
1529
1530 <itemizedlist>
1531   <listitem>
1532     <para>not interactive: GiNaC programs have to be written in 
1533     an editor, compiled and executed. You cannot play with 
1534     expressions interactively.  However, such an extension is not
1535     inherently forbidden by design.  In fact, two interactive
1536     interfaces are possible: First, a simple shell that exposes GiNaC's
1537     types to a command line can readily be written (and has been
1538     written) and second, as a more consistent approach we plan
1539     an integration with the <literal>CINT</literal>
1540     <literal>C++</literal> interpreter.</para>
1541   </listitem>
1542   <listitem>
1543     <para>advanced features: GiNaC cannot compete with a program
1544     like <literal>Reduce</literal> which exists for more than
1545     30 years now or <literal>Maple</literal> which grows since 
1546     1981 by the work of dozens of programmers, with respect to
1547     mathematical features. Integration, factorization, non-trivial
1548     simplifications, limits etc. are missing in GiNaC (and are not
1549     planned for the near future).</para>
1550   </listitem>
1551   <listitem>
1552     <para>portability: While the GiNaC library itself is designed
1553     to avoid any platform dependent features (it should compile
1554     on any ANSI compliant <literal>C++</literal> compiler), the
1555     currently used version of the CLN library (fast large integer and
1556     arbitrary precision arithmetics) can be compiled only on systems
1557     with a recently new <literal>C++</literal> compiler from the
1558     GNU Compiler Collection (<literal>GCC</literal>).  GiNaC uses
1559     recent language features like explicit constructors, mutable
1560     members, RTTI, dynamic_casts and STL, so ANSI compliance is meant
1561     literally.  Recent <literal>GCC</literal> versions starting at
1562     2.95, although itself not yet ANSI compliant, support all needed
1563     features.
1564     </para>
1565   </listitem>
1566 </itemizedlist>
1567 </para>
1568
1569 <sect1><title>Why <literal>C++</literal>?</title>
1570
1571 <para>Why did we choose to implement GiNaC in <literal>C++</literal>
1572 instead of <literal>Java</literal> or any other language?
1573 <literal>C++</literal> is not perfect: type checking is not strict
1574 (casting is possible), separation between interface and implementation
1575 is not complete, object oriented design is not enforced.  The main
1576 reason is the often scolded feature of operator overloading in
1577 <literal>C++</literal>. While it may be true that operating on classes
1578 with a <literal>+</literal> operator is rarely meaningful, it is
1579 perfectly suited for algebraic expressions. Writing 3x+5y as
1580 <literal>3*x+5*y</literal> instead of
1581 <literal>x.times(3).plus(y.times(5))</literal> looks much more
1582 natural. Furthermore, the main developers are more familiar with
1583 <literal>C++</literal> than with any other programming
1584 language.</para>
1585
1586 </chapter>
1587
1588
1589 <bibliography>
1590 <bibliodiv>
1591
1592 <biblioentry>
1593   <bookbiblio>
1594     <title>ISO/IEC 14882:1998</title>
1595     <subtitle>Programming Languages: C++</subtitle>
1596   </bookbiblio>
1597 </biblioentry>
1598
1599 <bibliomixed>
1600   <title>CLN: A Class Library for Numbers</title>
1601   <authorgroup>
1602     <author>
1603       <firstname>Bruno</firstname><surname>Haible</surname>
1604       <affiliation><address><email>haible@ilog.fr</email></address></affiliation>
1605     </author>
1606   </authorgroup>
1607 </bibliomixed>
1608
1609 <biblioentry>
1610   <bookbiblio>
1611     <title>The C++ Programming Language</title>
1612     <authorgroup><author><firstname>Bjarne</firstname><surname>Stroustrup</surname></author></authorgroup>
1613     <edition>3</edition>
1614     <isbn>0-201-88954-4</isbn>
1615     <publisher><publishername>Addison Wesley</publishername></publisher>
1616   </bookbiblio>
1617 </biblioentry>
1618
1619 <biblioentry>
1620   <bookbiblio>
1621     <title>Algorithms for Computer Algebra</title>
1622     <authorgroup>
1623       <author><firstname>Keith</firstname><othername>O.</othername><surname>Geddes</surname></author>
1624       <author><firstname>Stephen</firstname><othername>R.</othername><surname>Czapor</surname></author>
1625       <author><firstname>George</firstname><surname>Labahn</surname></author>
1626     </authorgroup>
1627     <isbn>0-7923-9259-0</isbn>
1628     <pubdate>1992</pubdate>
1629     <publisher>
1630       <publishername>Kluwer Academic Publishers</publishername>
1631       <address><city>Norwell</city>, <state>Massachusetts</state></address>
1632     </publisher>
1633   </bookbiblio>
1634 </biblioentry>
1635
1636 <biblioentry>
1637   <bookbiblio>
1638     <title>Computer Algebra</title>
1639     <subtitle>Systems and Algorithms for Algebraic Computation</subtitle>
1640     <authorgroup>
1641       <author><firstname>J.</firstname><othername>H.</othername><surname>Davenport</surname></author>
1642       <author><firstname>Y.</firstname><surname>Siret</surname></author>
1643       <author><firstname>E.</firstname><surname>Tournier</surname></author>
1644     </authorgroup>
1645     <isbn>0-12-204230-1</isbn>
1646     <pubdate>1988</pubdate>
1647     <publisher>
1648       <publishername>Academic Press</publishername>
1649       <address><city>London</city></address>
1650     </publisher>
1651   </bookbiblio>
1652 </biblioentry>
1653
1654 </bibliodiv>
1655 </bibliography>
1656
1657 </book>