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