added tutorial section about indexed objects
[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 on
6 @afourpaper
7 @c For `info' only.
8 @paragraphindent 0
9 @c For TeX only.
10 @iftex
11 @c I hate putting "@noindent" in front of every paragraph.
12 @parindent=0pt
13 @end iftex
14 @c %**end of header
15
16 @include version.texi
17
18 @direntry
19 * ginac: (ginac).                   C++ library for symbolic computation.
20 @end direntry
21
22 @ifinfo
23 This is a tutorial that documents GiNaC @value{VERSION}, an open
24 framework for symbolic computation within the C++ programming language.
25
26 Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
27
28 Permission is granted to make and distribute verbatim copies of
29 this manual provided the copyright notice and this permission notice
30 are preserved on all copies.
31
32 @ignore
33 Permission is granted to process this file through TeX and print the
34 results, provided the printed document carries copying permission
35 notice identical to this one except for the removal of this paragraph
36
37 @end ignore
38 Permission is granted to copy and distribute modified versions of this
39 manual under the conditions for verbatim copying, provided that the entire
40 resulting derived work is distributed under the terms of a permission
41 notice identical to this one.
42 @end ifinfo
43
44 @finalout
45 @c finalout prevents ugly black rectangles on overfull hbox lines
46 @titlepage
47 @title GiNaC @value{VERSION}
48 @subtitle An open framework for symbolic computation within the C++ programming language
49 @subtitle @value{UPDATED}
50 @author The GiNaC Group:
51 @author Christian Bauer, Alexander Frink, Richard Kreckel
52
53 @page
54 @vskip 0pt plus 1filll
55 Copyright @copyright{} 1999-2001 Johannes Gutenberg University Mainz, Germany
56 @sp 2
57 Permission is granted to make and distribute verbatim copies of
58 this manual provided the copyright notice and this permission notice
59 are preserved on all copies.
60
61 Permission is granted to copy and distribute modified versions of this
62 manual under the conditions for verbatim copying, provided that the entire
63 resulting derived work is distributed under the terms of a permission
64 notice identical to this one.
65 @end titlepage
66
67 @page
68 @contents
69
70 @page
71
72
73 @node Top, Introduction, (dir), (dir)
74 @c    node-name, next, previous, up
75 @top GiNaC
76
77 This is a tutorial that documents GiNaC @value{VERSION}, an open
78 framework for symbolic computation within the C++ programming language.
79
80 @menu
81 * Introduction::                 GiNaC's purpose.
82 * A Tour of GiNaC::              A quick tour of the library.
83 * Installation::                 How to install the package.
84 * Basic Concepts::               Description of fundamental classes.
85 * Methods and Functions::        Algorithms for symbolic manipulations.
86 * Extending GiNaC::              How to extend the library.
87 * A Comparison With Other CAS::  Compares GiNaC to traditional CAS.
88 * Internal Structures::          Description of some internal structures.
89 * Package Tools::                Configuring packages to work with GiNaC.
90 * Bibliography::
91 * Concept Index::
92 @end menu
93
94
95 @node Introduction, A Tour of GiNaC, Top, Top
96 @c    node-name, next, previous, up
97 @chapter Introduction
98 @cindex history of GiNaC
99
100 The motivation behind GiNaC derives from the observation that most
101 present day computer algebra systems (CAS) are linguistically and
102 semantically impoverished.  Although they are quite powerful tools for
103 learning math and solving particular problems they lack modern
104 linguistical structures that allow for the creation of large-scale
105 projects.  GiNaC is an attempt to overcome this situation by extending a
106 well established and standardized computer language (C++) by some
107 fundamental symbolic capabilities, thus allowing for integrated systems
108 that embed symbolic manipulations together with more established areas
109 of computer science (like computation-intense numeric applications,
110 graphical interfaces, etc.) under one roof.
111
112 The particular problem that led to the writing of the GiNaC framework is
113 still a very active field of research, namely the calculation of higher
114 order corrections to elementary particle interactions.  There,
115 theoretical physicists are interested in matching present day theories
116 against experiments taking place at particle accelerators.  The
117 computations involved are so complex they call for a combined symbolical
118 and numerical approach.  This turned out to be quite difficult to
119 accomplish with the present day CAS we have worked with so far and so we
120 tried to fill the gap by writing GiNaC.  But of course its applications
121 are in no way restricted to theoretical physics.
122
123 This tutorial is intended for the novice user who is new to GiNaC but
124 already has some background in C++ programming.  However, since a
125 hand-made documentation like this one is difficult to keep in sync with
126 the development, the actual documentation is inside the sources in the
127 form of comments.  That documentation may be parsed by one of the many
128 Javadoc-like documentation systems.  If you fail at generating it you
129 may access it from @uref{http://www.ginac.de/reference/, the GiNaC home
130 page}.  It is an invaluable resource not only for the advanced user who
131 wishes to extend the system (or chase bugs) but for everybody who wants
132 to comprehend the inner workings of GiNaC.  This little tutorial on the
133 other hand only covers the basic things that are unlikely to change in
134 the near future.
135
136 @section License
137 The GiNaC framework for symbolic computation within the C++ programming
138 language is Copyright @copyright{} 1999-2001 Johannes Gutenberg
139 University Mainz, Germany.
140
141 This program is free software; you can redistribute it and/or
142 modify it under the terms of the GNU General Public License as
143 published by the Free Software Foundation; either version 2 of the
144 License, or (at your option) any later version.
145
146 This program is distributed in the hope that it will be useful, but
147 WITHOUT ANY WARRANTY; without even the implied warranty of
148 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
149 General Public License for more details.
150
151 You should have received a copy of the GNU General Public License
152 along with this program; see the file COPYING.  If not, write to the
153 Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
154 MA 02111-1307, USA.
155
156
157 @node A Tour of GiNaC, How to use it from within C++, Introduction, Top
158 @c    node-name, next, previous, up
159 @chapter A Tour of GiNaC
160
161 This quick tour of GiNaC wants to arise your interest in the
162 subsequent chapters by showing off a bit.  Please excuse us if it
163 leaves many open questions.
164
165 @menu
166 * How to use it from within C++::  Two simple examples.
167 * What it can do for you::         A Tour of GiNaC's features.
168 @end menu
169
170
171 @node How to use it from within C++, What it can do for you, A Tour of GiNaC, A Tour of GiNaC
172 @c    node-name, next, previous, up
173 @section How to use it from within C++
174
175 The GiNaC open framework for symbolic computation within the C++ programming
176 language does not try to define a language of its own as conventional
177 CAS do.  Instead, it extends the capabilities of C++ by symbolic
178 manipulations.  Here is how to generate and print a simple (and rather
179 pointless) bivariate polynomial with some large coefficients:
180
181 @example
182 #include <ginac/ginac.h>
183 using namespace GiNaC;
184
185 int main()
186 @{
187     symbol x("x"), y("y");
188     ex poly;
189
190     for (int i=0; i<3; ++i)
191         poly += factorial(i+16)*pow(x,i)*pow(y,2-i);
192
193     cout << poly << endl;
194     return 0;
195 @}
196 @end example
197
198 Assuming the file is called @file{hello.cc}, on our system we can compile
199 and run it like this:
200
201 @example
202 $ c++ hello.cc -o hello -lcln -lginac
203 $ ./hello
204 355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
205 @end example
206
207 (@xref{Package Tools}, for tools that help you when creating a software
208 package that uses GiNaC.)
209
210 @cindex Hermite polynomial
211 Next, there is a more meaningful C++ program that calls a function which
212 generates Hermite polynomials in a specified free variable.
213
214 @example
215 #include <ginac/ginac.h>
216 using namespace GiNaC;
217
218 ex HermitePoly(const symbol & x, int n)
219 @{
220     ex HKer=exp(-pow(x, 2));
221     // uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2)
222     return normal(pow(-1, n) * diff(HKer, x, n) / HKer);
223 @}
224
225 int main()
226 @{
227     symbol z("z");
228
229     for (int i=0; i<6; ++i)
230         cout << "H_" << i << "(z) == " << HermitePoly(z,i) << endl;
231
232     return 0;
233 @}
234 @end example
235
236 When run, this will type out
237
238 @example
239 H_0(z) == 1
240 H_1(z) == 2*z
241 H_2(z) == 4*z^2-2
242 H_3(z) == -12*z+8*z^3
243 H_4(z) == -48*z^2+16*z^4+12
244 H_5(z) == 120*z-160*z^3+32*z^5
245 @end example
246
247 This method of generating the coefficients is of course far from optimal
248 for production purposes.
249
250 In order to show some more examples of what GiNaC can do we will now use
251 the @command{ginsh}, a simple GiNaC interactive shell that provides a
252 convenient window into GiNaC's capabilities.
253
254
255 @node What it can do for you, Installation, How to use it from within C++, A Tour of GiNaC
256 @c    node-name, next, previous, up
257 @section What it can do for you
258
259 @cindex @command{ginsh}
260 After invoking @command{ginsh} one can test and experiment with GiNaC's
261 features much like in other Computer Algebra Systems except that it does
262 not provide programming constructs like loops or conditionals.  For a
263 concise description of the @command{ginsh} syntax we refer to its
264 accompanied man page. Suffice to say that assignments and comparisons in
265 @command{ginsh} are written as they are in C, i.e. @code{=} assigns and
266 @code{==} compares.
267
268 It can manipulate arbitrary precision integers in a very fast way.
269 Rational numbers are automatically converted to fractions of coprime
270 integers:
271
272 @example
273 > x=3^150;
274 369988485035126972924700782451696644186473100389722973815184405301748249
275 > y=3^149;
276 123329495011708990974900260817232214728824366796574324605061468433916083
277 > x/y;
278 3
279 > y/x;
280 1/3
281 @end example
282
283 Exact numbers are always retained as exact numbers and only evaluated as
284 floating point numbers if requested.  For instance, with numeric
285 radicals is dealt pretty much as with symbols.  Products of sums of them
286 can be expanded:
287
288 @example
289 > expand((1+a^(1/5)-a^(2/5))^3);
290 1+3*a+3*a^(1/5)-5*a^(3/5)-a^(6/5)
291 > expand((1+3^(1/5)-3^(2/5))^3);
292 10-5*3^(3/5)
293 > evalf((1+3^(1/5)-3^(2/5))^3);
294 0.33408977534118624228
295 @end example
296
297 The function @code{evalf} that was used above converts any number in
298 GiNaC's expressions into floating point numbers.  This can be done to
299 arbitrary predefined accuracy:
300
301 @example
302 > evalf(1/7);
303 0.14285714285714285714
304 > Digits=150;
305 150
306 > evalf(1/7);
307 0.1428571428571428571428571428571428571428571428571428571428571428571428
308 5714285714285714285714285714285714285
309 @end example
310
311 Exact numbers other than rationals that can be manipulated in GiNaC
312 include predefined constants like Archimedes' @code{Pi}.  They can both
313 be used in symbolic manipulations (as an exact number) as well as in
314 numeric expressions (as an inexact number):
315
316 @example
317 > a=Pi^2+x;
318 x+Pi^2
319 > evalf(a);
320 9.869604401089358619+x
321 > x=2;
322 2
323 > evalf(a);
324 11.869604401089358619
325 @end example
326
327 Built-in functions evaluate immediately to exact numbers if
328 this is possible.  Conversions that can be safely performed are done
329 immediately; conversions that are not generally valid are not done:
330
331 @example
332 > cos(42*Pi);
333 1
334 > cos(acos(x));
335 x
336 > acos(cos(x));
337 acos(cos(x))
338 @end example
339
340 (Note that converting the last input to @code{x} would allow one to
341 conclude that @code{42*Pi} is equal to @code{0}.)
342
343 Linear equation systems can be solved along with basic linear
344 algebra manipulations over symbolic expressions.  In C++ GiNaC offers
345 a matrix class for this purpose but we can see what it can do using
346 @command{ginsh}'s notation of double brackets to type them in:
347
348 @example
349 > lsolve(a+x*y==z,x);
350 y^(-1)*(z-a);
351 > lsolve([3*x+5*y == 7, -2*x+10*y == -5], [x, y]);
352 [x==19/8,y==-1/40]
353 > M = [[ [[1, 3]], [[-3, 2]] ]];
354 [[ [[1,3]], [[-3,2]] ]]
355 > determinant(M);
356 11
357 > charpoly(M,lambda);
358 lambda^2-3*lambda+11
359 @end example
360
361 Multivariate polynomials and rational functions may be expanded,
362 collected and normalized (i.e. converted to a ratio of two coprime 
363 polynomials):
364
365 @example
366 > a = x^4 + 2*x^2*y^2 + 4*x^3*y + 12*x*y^3 - 3*y^4;
367 -3*y^4+x^4+12*x*y^3+2*x^2*y^2+4*x^3*y
368 > b = x^2 + 4*x*y - y^2;
369 -y^2+x^2+4*x*y
370 > expand(a*b);
371 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
372 > collect(a*b,x);
373 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)
374 > normal(a/b);
375 3*y^2+x^2
376 @end example
377
378 You can differentiate functions and expand them as Taylor or Laurent
379 series in a very natural syntax (the second argument of @code{series} is
380 a relation defining the evaluation point, the third specifies the
381 order):
382
383 @cindex Zeta function
384 @example
385 > diff(tan(x),x);
386 tan(x)^2+1
387 > series(sin(x),x==0,4);
388 x-1/6*x^3+Order(x^4)
389 > series(1/tan(x),x==0,4);
390 x^(-1)-1/3*x+Order(x^2)
391 > series(tgamma(x),x==0,3);
392 x^(-1)-Euler+(1/12*Pi^2+1/2*Euler^2)*x+
393 (-1/3*zeta(3)-1/12*Pi^2*Euler-1/6*Euler^3)*x^2+Order(x^3)
394 > evalf(");
395 x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
396 -(0.90747907608088628905)*x^2+Order(x^3)
397 > series(tgamma(2*sin(x)-2),x==Pi/2,6);
398 -(x-1/2*Pi)^(-2)+(-1/12*Pi^2-1/2*Euler^2-1/240)*(x-1/2*Pi)^2
399 -Euler-1/12+Order((x-1/2*Pi)^3)
400 @end example
401
402 Here we have made use of the @command{ginsh}-command @code{"} to pop the
403 previously evaluated element from @command{ginsh}'s internal stack.
404
405 If you ever wanted to convert units in C or C++ and found this is
406 cumbersome, here is the solution.  Symbolic types can always be used as
407 tags for different types of objects.  Converting from wrong units to the
408 metric system is now easy:
409
410 @example
411 > in=.0254*m;
412 0.0254*m
413 > lb=.45359237*kg;
414 0.45359237*kg
415 > 200*lb/in^2;
416 140613.91592783185568*kg*m^(-2)
417 @end example
418
419
420 @node Installation, Prerequisites, What it can do for you, Top
421 @c    node-name, next, previous, up
422 @chapter Installation
423
424 @cindex CLN
425 GiNaC's installation follows the spirit of most GNU software. It is
426 easily installed on your system by three steps: configuration, build,
427 installation.
428
429 @menu
430 * Prerequisites::                Packages upon which GiNaC depends.
431 * Configuration::                How to configure GiNaC.
432 * Building GiNaC::               How to compile GiNaC.
433 * Installing GiNaC::             How to install GiNaC on your system.
434 @end menu
435
436
437 @node Prerequisites, Configuration, Installation, Installation
438 @c    node-name, next, previous, up
439 @section Prerequisites
440
441 In order to install GiNaC on your system, some prerequisites need to be
442 met.  First of all, you need to have a C++-compiler adhering to the
443 ANSI-standard @cite{ISO/IEC 14882:1998(E)}.  We used @acronym{GCC} for
444 development so if you have a different compiler you are on your own.
445 For the configuration to succeed you need a Posix compliant shell
446 installed in @file{/bin/sh}, GNU @command{bash} is fine.  Perl is needed
447 by the built process as well, since some of the source files are
448 automatically generated by Perl scripts.  Last but not least, Bruno
449 Haible's library @acronym{CLN} is extensively used and needs to be
450 installed on your system.  Please get it either from
451 @uref{ftp://ftp.santafe.edu/pub/gnu/}, from
452 @uref{ftp://ftpthep.physik.uni-mainz.de/pub/gnu/, GiNaC's FTP site} or
453 from @uref{ftp://ftp.ilog.fr/pub/Users/haible/gnu/, Bruno Haible's FTP
454 site} (it is covered by GPL) and install it prior to trying to install
455 GiNaC.  The configure script checks if it can find it and if it cannot
456 it will refuse to continue.
457
458
459 @node Configuration, Building GiNaC, Prerequisites, Installation
460 @c    node-name, next, previous, up
461 @section Configuration
462 @cindex configuration
463 @cindex Autoconf
464
465 To configure GiNaC means to prepare the source distribution for
466 building.  It is done via a shell script called @command{configure} that
467 is shipped with the sources and was originally generated by GNU
468 Autoconf.  Since a configure script generated by GNU Autoconf never
469 prompts, all customization must be done either via command line
470 parameters or environment variables.  It accepts a list of parameters,
471 the complete set of which can be listed by calling it with the
472 @option{--help} option.  The most important ones will be shortly
473 described in what follows:
474
475 @itemize @bullet
476
477 @item
478 @option{--disable-shared}: When given, this option switches off the
479 build of a shared library, i.e. a @file{.so} file.  This may be convenient
480 when developing because it considerably speeds up compilation.
481
482 @item
483 @option{--prefix=@var{PREFIX}}: The directory where the compiled library
484 and headers are installed. It defaults to @file{/usr/local} which means
485 that the library is installed in the directory @file{/usr/local/lib},
486 the header files in @file{/usr/local/include/ginac} and the documentation
487 (like this one) into @file{/usr/local/share/doc/GiNaC}.
488
489 @item
490 @option{--libdir=@var{LIBDIR}}: Use this option in case you want to have
491 the library installed in some other directory than
492 @file{@var{PREFIX}/lib/}.
493
494 @item
495 @option{--includedir=@var{INCLUDEDIR}}: Use this option in case you want
496 to have the header files installed in some other directory than
497 @file{@var{PREFIX}/include/ginac/}. For instance, if you specify
498 @option{--includedir=/usr/include} you will end up with the header files
499 sitting in the directory @file{/usr/include/ginac/}. Note that the
500 subdirectory @file{ginac} is enforced by this process in order to
501 keep the header files separated from others.  This avoids some
502 clashes and allows for an easier deinstallation of GiNaC. This ought
503 to be considered A Good Thing (tm).
504
505 @item
506 @option{--datadir=@var{DATADIR}}: This option may be given in case you
507 want to have the documentation installed in some other directory than
508 @file{@var{PREFIX}/share/doc/GiNaC/}.
509
510 @end itemize
511
512 In addition, you may specify some environment variables.
513 @env{CXX} holds the path and the name of the C++ compiler
514 in case you want to override the default in your path.  (The
515 @command{configure} script searches your path for @command{c++},
516 @command{g++}, @command{gcc}, @command{CC}, @command{cxx}
517 and @command{cc++} in that order.)  It may be very useful to
518 define some compiler flags with the @env{CXXFLAGS} environment
519 variable, like optimization, debugging information and warning
520 levels.  If omitted, it defaults to @option{-g -O2}.
521
522 The whole process is illustrated in the following two
523 examples. (Substitute @command{setenv @var{VARIABLE} @var{value}} for
524 @command{export @var{VARIABLE}=@var{value}} if the Berkeley C shell is
525 your login shell.)
526
527 Here is a simple configuration for a site-wide GiNaC library assuming
528 everything is in default paths:
529
530 @example
531 $ export CXXFLAGS="-Wall -O2"
532 $ ./configure
533 @end example
534
535 And here is a configuration for a private static GiNaC library with
536 several components sitting in custom places (site-wide @acronym{GCC} and
537 private @acronym{CLN}).  The compiler is pursuaded to be picky and full
538 assertions and debugging information are switched on:
539
540 @example
541 $ export CXX=/usr/local/gnu/bin/c++
542 $ export CPPFLAGS="$(CPPFLAGS) -I$(HOME)/include"
543 $ export CXXFLAGS="$(CXXFLAGS) -DDO_GINAC_ASSERT -ggdb -Wall -ansi -pedantic"
544 $ export LDFLAGS="$(LDFLAGS) -L$(HOME)/lib"
545 $ ./configure --disable-shared --prefix=$(HOME)
546 @end example
547
548
549 @node Building GiNaC, Installing GiNaC, Configuration, Installation
550 @c    node-name, next, previous, up
551 @section Building GiNaC
552 @cindex building GiNaC
553
554 After proper configuration you should just build the whole
555 library by typing
556 @example
557 $ make
558 @end example
559 at the command prompt and go for a cup of coffee.  The exact time it
560 takes to compile GiNaC depends not only on the speed of your machines
561 but also on other parameters, for instance what value for @env{CXXFLAGS}
562 you entered.  Optimization may be very time-consuming.
563
564 Just to make sure GiNaC works properly you may run a collection of
565 regression tests by typing
566
567 @example
568 $ make check
569 @end example
570
571 This will compile some sample programs, run them and check the output
572 for correctness.  The regression tests fall in three categories.  First,
573 the so called @emph{exams} are performed, simple tests where some
574 predefined input is evaluated (like a pupils' exam).  Second, the
575 @emph{checks} test the coherence of results among each other with
576 possible random input.  Third, some @emph{timings} are performed, which
577 benchmark some predefined problems with different sizes and display the
578 CPU time used in seconds.  Each individual test should return a message
579 @samp{passed}.  This is mostly intended to be a QA-check if something
580 was broken during development, not a sanity check of your system.  Some
581 of the tests in sections @emph{checks} and @emph{timings} may require
582 insane amounts of memory and CPU time.  Feel free to kill them if your
583 machine catches fire.  Another quite important intent is to allow people
584 to fiddle around with optimization.
585
586 Generally, the top-level Makefile runs recursively to the
587 subdirectories.  It is therfore safe to go into any subdirectory
588 (@code{doc/}, @code{ginsh/}, ...) and simply type @code{make}
589 @var{target} there in case something went wrong.
590
591
592 @node Installing GiNaC, Basic Concepts, Building GiNaC, Installation
593 @c    node-name, next, previous, up
594 @section Installing GiNaC
595 @cindex installation
596
597 To install GiNaC on your system, simply type
598
599 @example
600 $ make install
601 @end example
602
603 As described in the section about configuration the files will be
604 installed in the following directories (the directories will be created
605 if they don't already exist):
606
607 @itemize @bullet
608
609 @item
610 @file{libginac.a} will go into @file{@var{PREFIX}/lib/} (or
611 @file{@var{LIBDIR}}) which defaults to @file{/usr/local/lib/}.
612 So will @file{libginac.so} unless the configure script was
613 given the option @option{--disable-shared}.  The proper symlinks
614 will be established as well.
615
616 @item
617 All the header files will be installed into @file{@var{PREFIX}/include/ginac/}
618 (or @file{@var{INCLUDEDIR}/ginac/}, if specified).
619
620 @item
621 All documentation (HTML and Postscript) will be stuffed into
622 @file{@var{PREFIX}/share/doc/GiNaC/} (or
623 @file{@var{DATADIR}/doc/GiNaC/}, if @var{DATADIR} was specified).
624
625 @end itemize
626
627 For the sake of completeness we will list some other useful make
628 targets: @command{make clean} deletes all files generated by
629 @command{make}, i.e. all the object files.  In addition @command{make
630 distclean} removes all files generated by the configuration and
631 @command{make maintainer-clean} goes one step further and deletes files
632 that may require special tools to rebuild (like the @command{libtool}
633 for instance).  Finally @command{make uninstall} removes the installed
634 library, header files and documentation@footnote{Uninstallation does not
635 work after you have called @command{make distclean} since the
636 @file{Makefile} is itself generated by the configuration from
637 @file{Makefile.in} and hence deleted by @command{make distclean}.  There
638 are two obvious ways out of this dilemma.  First, you can run the
639 configuration again with the same @var{PREFIX} thus creating a
640 @file{Makefile} with a working @samp{uninstall} target.  Second, you can
641 do it by hand since you now know where all the files went during
642 installation.}.
643
644
645 @node Basic Concepts, Expressions, Installing GiNaC, Top
646 @c    node-name, next, previous, up
647 @chapter Basic Concepts
648
649 This chapter will describe the different fundamental objects that can be
650 handled by GiNaC.  But before doing so, it is worthwhile introducing you
651 to the more commonly used class of expressions, representing a flexible
652 meta-class for storing all mathematical objects.
653
654 @menu
655 * Expressions::                  The fundamental GiNaC class.
656 * The Class Hierarchy::          Overview of GiNaC's classes.
657 * Symbols::                      Symbolic objects.
658 * Numbers::                      Numerical objects.
659 * Constants::                    Pre-defined constants.
660 * Fundamental containers::       The power, add and mul classes.
661 * Lists::                        Lists of expressions.
662 * Mathematical functions::       Mathematical functions.
663 * Relations::                    Equality, Inequality and all that.
664 * Indexed objects::              Handling indexed quantities.
665 @end menu
666
667
668 @node Expressions, The Class Hierarchy, Basic Concepts, Basic Concepts
669 @c    node-name, next, previous, up
670 @section Expressions
671 @cindex expression (class @code{ex})
672 @cindex @code{has()}
673
674 The most common class of objects a user deals with is the expression
675 @code{ex}, representing a mathematical object like a variable, number,
676 function, sum, product, etc...  Expressions may be put together to form
677 new expressions, passed as arguments to functions, and so on.  Here is a
678 little collection of valid expressions:
679
680 @example
681 ex MyEx1 = 5;                       // simple number
682 ex MyEx2 = x + 2*y;                 // polynomial in x and y
683 ex MyEx3 = (x + 1)/(x - 1);         // rational expression
684 ex MyEx4 = sin(x + 2*y) + 3*z + 41; // containing a function
685 ex MyEx5 = MyEx4 + 1;               // similar to above
686 @end example
687
688 Expressions are handles to other more fundamental objects, that often
689 contain other expressions thus creating a tree of expressions
690 (@xref{Internal Structures}, for particular examples).  Most methods on
691 @code{ex} therefore run top-down through such an expression tree.  For
692 example, the method @code{has()} scans recursively for occurrences of
693 something inside an expression.  Thus, if you have declared @code{MyEx4}
694 as in the example above @code{MyEx4.has(y)} will find @code{y} inside
695 the argument of @code{sin} and hence return @code{true}.
696
697 The next sections will outline the general picture of GiNaC's class
698 hierarchy and describe the classes of objects that are handled by
699 @code{ex}.
700
701
702 @node The Class Hierarchy, Symbols, Expressions, Basic Concepts
703 @c    node-name, next, previous, up
704 @section The Class Hierarchy
705
706 GiNaC's class hierarchy consists of several classes representing
707 mathematical objects, all of which (except for @code{ex} and some
708 helpers) are internally derived from one abstract base class called
709 @code{basic}.  You do not have to deal with objects of class
710 @code{basic}, instead you'll be dealing with symbols, numbers,
711 containers of expressions and so on.
712
713 @cindex container
714 @cindex atom
715 To get an idea about what kinds of symbolic composits may be built we
716 have a look at the most important classes in the class hierarchy and
717 some of the relations among the classes:
718
719 @image{classhierarchy}
720
721 The abstract classes shown here (the ones without drop-shadow) are of no
722 interest for the user.  They are used internally in order to avoid code
723 duplication if two or more classes derived from them share certain
724 features.  An example is @code{expairseq}, a container for a sequence of
725 pairs each consisting of one expression and a number (@code{numeric}).
726 What @emph{is} visible to the user are the derived classes @code{add}
727 and @code{mul}, representing sums and products.  @xref{Internal
728 Structures}, where these two classes are described in more detail.  The
729 following table shortly summarizes what kinds of mathematical objects
730 are stored in the different classes:
731
732 @cartouche
733 @multitable @columnfractions .22 .78
734 @item @code{symbol} @tab Algebraic symbols @math{a}, @math{x}, @math{y}@dots{}
735 @item @code{constant} @tab Constants like 
736 @tex
737 $\pi$
738 @end tex
739 @ifnottex
740 @math{Pi}
741 @end ifnottex
742 @item @code{numeric} @tab All kinds of numbers, @math{42}, @math{7/3*I}, @math{3.14159}@dots{}
743 @item @code{add} @tab Sums like @math{x+y} or @math{a-(2*b)+3}
744 @item @code{mul} @tab Products like @math{x*y} or @math{2*a^2*(x+y+z)/b}
745 @item @code{power} @tab Exponentials such as @math{x^2}, @math{a^b}, 
746 @tex
747 $\sqrt{2}$
748 @end tex
749 @ifnottex
750 @code{sqrt(}@math{2}@code{)}
751 @end ifnottex
752 @dots{}
753 @item @code{pseries} @tab Power Series, e.g. @math{x-1/6*x^3+1/120*x^5+O(x^7)}
754 @item @code{function} @tab A symbolic function like @math{sin(2*x)}
755 @item @code{lst} @tab Lists of expressions [@math{x}, @math{2*y}, @math{3+z}]
756 @item @code{matrix} @tab @math{n}x@math{m} matrices of expressions
757 @item @code{relational} @tab A relation like the identity @math{x}@code{==}@math{y}
758 @item @code{indexed} @tab Indexed object like @math{A_ij}
759 @item @code{tensor} @tab Special tensor like the delta and metric tensors
760 @item @code{idx} @tab Index of an indexed object
761 @item @code{varidx} @tab Index with variance
762 @end multitable
763 @end cartouche
764
765 @node Symbols, Numbers, The Class Hierarchy, Basic Concepts
766 @c    node-name, next, previous, up
767 @section Symbols
768 @cindex @code{symbol} (class)
769 @cindex hierarchy of classes
770
771 @cindex atom
772 Symbols are for symbolic manipulation what atoms are for chemistry.  You
773 can declare objects of class @code{symbol} as any other object simply by
774 saying @code{symbol x,y;}.  There is, however, a catch in here having to
775 do with the fact that C++ is a compiled language.  The information about
776 the symbol's name is thrown away by the compiler but at a later stage
777 you may want to print expressions holding your symbols.  In order to
778 avoid confusion GiNaC's symbols are able to know their own name.  This
779 is accomplished by declaring its name for output at construction time in
780 the fashion @code{symbol x("x");}.  If you declare a symbol using the
781 default constructor (i.e. without string argument) the system will deal
782 out a unique name.  That name may not be suitable for printing but for
783 internal routines when no output is desired it is often enough.  We'll
784 come across examples of such symbols later in this tutorial.
785
786 This implies that the strings passed to symbols at construction time may
787 not be used for comparing two of them.  It is perfectly legitimate to
788 write @code{symbol x("x"),y("x");} but it is likely to lead into
789 trouble.  Here, @code{x} and @code{y} are different symbols and
790 statements like @code{x-y} will not be simplified to zero although the
791 output @code{x-x} looks funny.  Such output may also occur when there
792 are two different symbols in two scopes, for instance when you call a
793 function that declares a symbol with a name already existent in a symbol
794 in the calling function.  Again, comparing them (using @code{operator==}
795 for instance) will always reveal their difference.  Watch out, please.
796
797 @cindex @code{subs()}
798 Although symbols can be assigned expressions for internal reasons, you
799 should not do it (and we are not going to tell you how it is done).  If
800 you want to replace a symbol with something else in an expression, you
801 can use the expression's @code{.subs()} method (@xref{Substituting Symbols},
802 for more information).
803
804
805 @node Numbers, Constants, Symbols, Basic Concepts
806 @c    node-name, next, previous, up
807 @section Numbers
808 @cindex @code{numeric} (class)
809
810 @cindex GMP
811 @cindex CLN
812 @cindex rational
813 @cindex fraction
814 For storing numerical things, GiNaC uses Bruno Haible's library
815 @acronym{CLN}.  The classes therein serve as foundation classes for
816 GiNaC.  @acronym{CLN} stands for Class Library for Numbers or
817 alternatively for Common Lisp Numbers.  In order to find out more about
818 @acronym{CLN}'s internals the reader is refered to the documentation of
819 that library.  @inforef{Introduction, , cln}, for more
820 information. Suffice to say that it is by itself build on top of another
821 library, the GNU Multiple Precision library @acronym{GMP}, which is an
822 extremely fast library for arbitrary long integers and rationals as well
823 as arbitrary precision floating point numbers.  It is very commonly used
824 by several popular cryptographic applications.  @acronym{CLN} extends
825 @acronym{GMP} by several useful things: First, it introduces the complex
826 number field over either reals (i.e. floating point numbers with
827 arbitrary precision) or rationals.  Second, it automatically converts
828 rationals to integers if the denominator is unity and complex numbers to
829 real numbers if the imaginary part vanishes and also correctly treats
830 algebraic functions.  Third it provides good implementations of
831 state-of-the-art algorithms for all trigonometric and hyperbolic
832 functions as well as for calculation of some useful constants.
833
834 The user can construct an object of class @code{numeric} in several
835 ways.  The following example shows the four most important constructors.
836 It uses construction from C-integer, construction of fractions from two
837 integers, construction from C-float and construction from a string:
838
839 @example
840 #include <ginac/ginac.h>
841 using namespace GiNaC;
842
843 int main()
844 @{
845     numeric two(2);                       // exact integer 2
846     numeric r(2,3);                       // exact fraction 2/3
847     numeric e(2.71828);                   // floating point number
848     numeric p("3.1415926535897932385");   // floating point number
849     // Trott's constant in scientific notation:
850     numeric trott("1.0841015122311136151E-2");
851     
852     cout << two*p << endl;  // floating point 6.283...
853 @}
854 @end example
855
856 Note that all those constructors are @emph{explicit} which means you are
857 not allowed to write @code{numeric two=2;}.  This is because the basic
858 objects to be handled by GiNaC are the expressions @code{ex} and we want
859 to keep things simple and wish objects like @code{pow(x,2)} to be
860 handled the same way as @code{pow(x,a)}, which means that we need to
861 allow a general @code{ex} as base and exponent.  Therefore there is an
862 implicit constructor from C-integers directly to expressions handling
863 numerics at work in most of our examples.  This design really becomes
864 convenient when one declares own functions having more than one
865 parameter but it forbids using implicit constructors because that would
866 lead to compile-time ambiguities.
867
868 It may be tempting to construct numbers writing @code{numeric r(3/2)}.
869 This would, however, call C's built-in operator @code{/} for integers
870 first and result in a numeric holding a plain integer 1.  @strong{Never
871 use the operator @code{/} on integers} unless you know exactly what you
872 are doing!  Use the constructor from two integers instead, as shown in
873 the example above.  Writing @code{numeric(1)/2} may look funny but works
874 also.
875
876 @cindex @code{Digits}
877 @cindex accuracy
878 We have seen now the distinction between exact numbers and floating
879 point numbers.  Clearly, the user should never have to worry about
880 dynamically created exact numbers, since their `exactness' always
881 determines how they ought to be handled, i.e. how `long' they are.  The
882 situation is different for floating point numbers.  Their accuracy is
883 controlled by one @emph{global} variable, called @code{Digits}.  (For
884 those readers who know about Maple: it behaves very much like Maple's
885 @code{Digits}).  All objects of class numeric that are constructed from
886 then on will be stored with a precision matching that number of decimal
887 digits:
888
889 @example
890 #include <ginac/ginac.h>
891 using namespace GiNaC;
892
893 void foo()
894 @{
895     numeric three(3.0), one(1.0);
896     numeric x = one/three;
897
898     cout << "in " << Digits << " digits:" << endl;
899     cout << x << endl;
900     cout << Pi.evalf() << endl;
901 @}
902
903 int main()
904 @{
905     foo();
906     Digits = 60;
907     foo();
908     return 0;
909 @}
910 @end example
911
912 The above example prints the following output to screen:
913
914 @example
915 in 17 digits:
916 0.333333333333333333
917 3.14159265358979324
918 in 60 digits:
919 0.333333333333333333333333333333333333333333333333333333333333333333
920 3.14159265358979323846264338327950288419716939937510582097494459231
921 @end example
922
923 It should be clear that objects of class @code{numeric} should be used
924 for constructing numbers or for doing arithmetic with them.  The objects
925 one deals with most of the time are the polymorphic expressions @code{ex}.
926
927 @subsection Tests on numbers
928
929 Once you have declared some numbers, assigned them to expressions and
930 done some arithmetic with them it is frequently desired to retrieve some
931 kind of information from them like asking whether that number is
932 integer, rational, real or complex.  For those cases GiNaC provides
933 several useful methods.  (Internally, they fall back to invocations of
934 certain CLN functions.)
935
936 As an example, let's construct some rational number, multiply it with
937 some multiple of its denominator and test what comes out:
938
939 @example
940 #include <ginac/ginac.h>
941 using namespace GiNaC;
942
943 // some very important constants:
944 const numeric twentyone(21);
945 const numeric ten(10);
946 const numeric five(5);
947
948 int main()
949 @{
950     numeric answer = twentyone;
951
952     answer /= five;
953     cout << answer.is_integer() << endl;  // false, it's 21/5
954     answer *= ten;
955     cout << answer.is_integer() << endl;  // true, it's 42 now!
956 @}
957 @end example
958
959 Note that the variable @code{answer} is constructed here as an integer
960 by @code{numeric}'s copy constructor but in an intermediate step it
961 holds a rational number represented as integer numerator and integer
962 denominator.  When multiplied by 10, the denominator becomes unity and
963 the result is automatically converted to a pure integer again.
964 Internally, the underlying @acronym{CLN} is responsible for this
965 behaviour and we refer the reader to @acronym{CLN}'s documentation.
966 Suffice to say that the same behaviour applies to complex numbers as
967 well as return values of certain functions.  Complex numbers are
968 automatically converted to real numbers if the imaginary part becomes
969 zero.  The full set of tests that can be applied is listed in the
970 following table.
971
972 @cartouche
973 @multitable @columnfractions .30 .70
974 @item @strong{Method} @tab @strong{Returns true if the object is@dots{}}
975 @item @code{.is_zero()}
976 @tab @dots{}equal to zero
977 @item @code{.is_positive()}
978 @tab @dots{}not complex and greater than 0
979 @item @code{.is_integer()}
980 @tab @dots{}a (non-complex) integer
981 @item @code{.is_pos_integer()}
982 @tab @dots{}an integer and greater than 0
983 @item @code{.is_nonneg_integer()}
984 @tab @dots{}an integer and greater equal 0
985 @item @code{.is_even()}
986 @tab @dots{}an even integer
987 @item @code{.is_odd()}
988 @tab @dots{}an odd integer
989 @item @code{.is_prime()}
990 @tab @dots{}a prime integer (probabilistic primality test)
991 @item @code{.is_rational()}
992 @tab @dots{}an exact rational number (integers are rational, too)
993 @item @code{.is_real()}
994 @tab @dots{}a real integer, rational or float (i.e. is not complex)
995 @item @code{.is_cinteger()}
996 @tab @dots{}a (complex) integer (such as @math{2-3*I})
997 @item @code{.is_crational()}
998 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
999 @end multitable
1000 @end cartouche
1001
1002
1003 @node Constants, Fundamental containers, Numbers, Basic Concepts
1004 @c    node-name, next, previous, up
1005 @section Constants
1006 @cindex @code{constant} (class)
1007
1008 @cindex @code{Pi}
1009 @cindex @code{Catalan}
1010 @cindex @code{Euler}
1011 @cindex @code{evalf()}
1012 Constants behave pretty much like symbols except that they return some
1013 specific number when the method @code{.evalf()} is called.
1014
1015 The predefined known constants are:
1016
1017 @cartouche
1018 @multitable @columnfractions .14 .30 .56
1019 @item @strong{Name} @tab @strong{Common Name} @tab @strong{Numerical Value (to 35 digits)}
1020 @item @code{Pi}
1021 @tab Archimedes' constant
1022 @tab 3.14159265358979323846264338327950288
1023 @item @code{Catalan}
1024 @tab Catalan's constant
1025 @tab 0.91596559417721901505460351493238411
1026 @item @code{Euler}
1027 @tab Euler's (or Euler-Mascheroni) constant
1028 @tab 0.57721566490153286060651209008240243
1029 @end multitable
1030 @end cartouche
1031
1032
1033 @node Fundamental containers, Lists, Constants, Basic Concepts
1034 @c    node-name, next, previous, up
1035 @section Fundamental containers: the @code{power}, @code{add} and @code{mul} classes
1036 @cindex polynomial
1037 @cindex @code{add}
1038 @cindex @code{mul}
1039 @cindex @code{power}
1040
1041 Simple polynomial expressions are written down in GiNaC pretty much like
1042 in other CAS or like expressions involving numerical variables in C.
1043 The necessary operators @code{+}, @code{-}, @code{*} and @code{/} have
1044 been overloaded to achieve this goal.  When you run the following
1045 program, the constructor for an object of type @code{mul} is
1046 automatically called to hold the product of @code{a} and @code{b} and
1047 then the constructor for an object of type @code{add} is called to hold
1048 the sum of that @code{mul} object and the number one:
1049
1050 @example
1051 #include <ginac/ginac.h>
1052 using namespace GiNaC;
1053
1054 int main()
1055 @{
1056     symbol a("a"), b("b");
1057     ex MyTerm = 1+a*b;
1058     // ...
1059 @}
1060 @end example
1061
1062 @cindex @code{pow()}
1063 For exponentiation, you have already seen the somewhat clumsy (though C-ish)
1064 statement @code{pow(x,2);} to represent @code{x} squared.  This direct
1065 construction is necessary since we cannot safely overload the constructor
1066 @code{^} in C++ to construct a @code{power} object.  If we did, it would
1067 have several counterintuitive effects:
1068
1069 @itemize @bullet
1070 @item
1071 Due to C's operator precedence, @code{2*x^2} would be parsed as @code{(2*x)^2}.
1072 @item
1073 Due to the binding of the operator @code{^}, @code{x^a^b} would result in
1074 @code{(x^a)^b}. This would be confusing since most (though not all) other CAS
1075 interpret this as @code{x^(a^b)}.
1076 @item
1077 Also, expressions involving integer exponents are very frequently used,
1078 which makes it even more dangerous to overload @code{^} since it is then
1079 hard to distinguish between the semantics as exponentiation and the one
1080 for exclusive or.  (It would be embarassing to return @code{1} where one
1081 has requested @code{2^3}.)
1082 @end itemize
1083
1084 @cindex @command{ginsh}
1085 All effects are contrary to mathematical notation and differ from the
1086 way most other CAS handle exponentiation, therefore overloading @code{^}
1087 is ruled out for GiNaC's C++ part.  The situation is different in
1088 @command{ginsh}, there the exponentiation-@code{^} exists.  (Also note
1089 that the other frequently used exponentiation operator @code{**} does
1090 not exist at all in C++).
1091
1092 To be somewhat more precise, objects of the three classes described
1093 here, are all containers for other expressions.  An object of class
1094 @code{power} is best viewed as a container with two slots, one for the
1095 basis, one for the exponent.  All valid GiNaC expressions can be
1096 inserted.  However, basic transformations like simplifying
1097 @code{pow(pow(x,2),3)} to @code{x^6} automatically are only performed
1098 when this is mathematically possible.  If we replace the outer exponent
1099 three in the example by some symbols @code{a}, the simplification is not
1100 safe and will not be performed, since @code{a} might be @code{1/2} and
1101 @code{x} negative.
1102
1103 Objects of type @code{add} and @code{mul} are containers with an
1104 arbitrary number of slots for expressions to be inserted.  Again, simple
1105 and safe simplifications are carried out like transforming
1106 @code{3*x+4-x} to @code{2*x+4}.
1107
1108 The general rule is that when you construct such objects, GiNaC
1109 automatically creates them in canonical form, which might differ from
1110 the form you typed in your program.  This allows for rapid comparison of
1111 expressions, since after all @code{a-a} is simply zero.  Note, that the
1112 canonical form is not necessarily lexicographical ordering or in any way
1113 easily guessable.  It is only guaranteed that constructing the same
1114 expression twice, either implicitly or explicitly, results in the same
1115 canonical form.
1116
1117
1118 @node Lists, Mathematical functions, Fundamental containers, Basic Concepts
1119 @c    node-name, next, previous, up
1120 @section Lists of expressions
1121 @cindex @code{lst} (class)
1122 @cindex lists
1123 @cindex @code{nops()}
1124 @cindex @code{op()}
1125 @cindex @code{append()}
1126 @cindex @code{prepend()}
1127
1128 The GiNaC class @code{lst} serves for holding a list of arbitrary expressions.
1129 These are sometimes used to supply a variable number of arguments of the same
1130 type to GiNaC methods such as @code{subs()} and @code{to_rational()}, so you
1131 should have a basic understanding about them.
1132
1133 Lists of up to 15 expressions can be directly constructed from single
1134 expressions:
1135
1136 @example
1137 @{
1138     symbol x("x"), y("y");
1139     lst l(x, 2, y, x+y);
1140     // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y'
1141     // ...
1142 @end example
1143
1144 Use the @code{nops()} method to determine the size (number of expressions) of
1145 a list and the @code{op()} method to access individual elements:
1146
1147 @example
1148     // ...
1149     cout << l.nops() << endl;                   // prints '4'
1150     cout << l.op(2) << " " << l.op(0) << endl;  // prints 'y x'
1151     // ...
1152 @end example
1153
1154 Finally you can append or prepend an expression to a list with the
1155 @code{append()} and @code{prepend()} methods:
1156
1157 @example
1158     // ...
1159     l.append(4*x);   // l is now [x, 2, y, x+y, 4*x]
1160     l.prepend(0);    // l is now [0, x, 2, y, x+y, 4*x]
1161 @}
1162 @end example
1163
1164
1165 @node Mathematical functions, Relations, Lists, Basic Concepts
1166 @c    node-name, next, previous, up
1167 @section Mathematical functions
1168 @cindex @code{function} (class)
1169 @cindex trigonometric function
1170 @cindex hyperbolic function
1171
1172 There are quite a number of useful functions hard-wired into GiNaC.  For
1173 instance, all trigonometric and hyperbolic functions are implemented
1174 (@xref{Built-in Functions}, for a complete list).
1175
1176 These functions are all objects of class @code{function}.  They accept one
1177 or more expressions as arguments and return one expression.  If the arguments
1178 are not numerical, the evaluation of the function may be halted, as it
1179 does in the next example:
1180
1181 @cindex Gamma function
1182 @cindex @code{subs()}
1183 @example
1184 #include <ginac/ginac.h>
1185 using namespace GiNaC;
1186
1187 int main()
1188 @{
1189     symbol x("x"), y("y");
1190     
1191     ex foo = x+y/2;
1192     cout << "tgamma(" << foo << ") -> " << tgamma(foo) << endl;
1193     ex bar = foo.subs(y==1);
1194     cout << "tgamma(" << bar << ") -> " << tgamma(bar) << endl;
1195     ex foobar = bar.subs(x==7);
1196     cout << "tgamma(" << foobar << ") -> " << tgamma(foobar) << endl;
1197 @}
1198 @end example
1199
1200 This program shows how the function returns itself twice and finally an
1201 expression that may be really useful:
1202
1203 @example
1204 tgamma(x+(1/2)*y) -> tgamma(x+(1/2)*y)
1205 tgamma(x+1/2) -> tgamma(x+1/2)
1206 tgamma(15/2) -> (135135/128)*Pi^(1/2)
1207 @end example
1208
1209 Besides evaluation most of these functions allow differentiation, series
1210 expansion and so on.  Read the next chapter in order to learn more about
1211 this.
1212
1213
1214 @node Relations, Indexed objects, Mathematical functions, Basic Concepts
1215 @c    node-name, next, previous, up
1216 @section Relations
1217 @cindex @code{relational} (class)
1218
1219 Sometimes, a relation holding between two expressions must be stored
1220 somehow.  The class @code{relational} is a convenient container for such
1221 purposes.  A relation is by definition a container for two @code{ex} and
1222 a relation between them that signals equality, inequality and so on.
1223 They are created by simply using the C++ operators @code{==}, @code{!=},
1224 @code{<}, @code{<=}, @code{>} and @code{>=} between two expressions.
1225
1226 @xref{Mathematical functions}, for examples where various applications
1227 of the @code{.subs()} method show how objects of class relational are
1228 used as arguments.  There they provide an intuitive syntax for
1229 substitutions.  They are also used as arguments to the @code{ex::series}
1230 method, where the left hand side of the relation specifies the variable
1231 to expand in and the right hand side the expansion point.  They can also
1232 be used for creating systems of equations that are to be solved for
1233 unknown variables.  But the most common usage of objects of this class
1234 is rather inconspicuous in statements of the form @code{if
1235 (expand(pow(a+b,2))==a*a+2*a*b+b*b) @{...@}}.  Here, an implicit
1236 conversion from @code{relational} to @code{bool} takes place.  Note,
1237 however, that @code{==} here does not perform any simplifications, hence
1238 @code{expand()} must be called explicitly.
1239
1240
1241 @node Indexed objects, Methods and Functions, Relations, Basic Concepts
1242 @c    node-name, next, previous, up
1243 @section Indexed objects
1244
1245 GiNaC allows you to handle expressions containing general indexed objects in
1246 arbitrary spaces. It is also able to canonicalize and simplify such
1247 expressions and perform symbolic dummy index summations. There are a number
1248 of predefined indexed objects provided, like delta and metric tensors.
1249
1250 There are few restrictions placed on indexed objects and their indices and
1251 it is easy to construct nonsense expressions, but our intention is to
1252 provide a general framework that allows you to implement algorithms with
1253 indexed quantities, getting in the way as little as possible.
1254
1255 @cindex @code{idx} (class)
1256 @cindex @code{indexed} (class)
1257 @subsection Indexed quantities and their indices
1258
1259 Indexed expressions in GiNaC are constructed of two special types of objects,
1260 @dfn{index objects} and @dfn{indexed objects}.
1261
1262 @itemize @bullet
1263
1264 @item Index objects are of class @code{idx} or a subclass. Every index has
1265 a @dfn{value} and a @dfn{dimension} (which is the dimension of the space
1266 the index lives in) which can both be arbitrary expressions but are usually
1267 a number or a simple symbol. In addition, indices of class @code{varidx} have
1268 a @dfn{variance} (they can be co- or contravariant).
1269
1270 @item Indexed objects are of class @code{indexed} or a subclass. They
1271 contain a @dfn{base expression} (which is the expression being indexed), and
1272 one or more indices.
1273
1274 @end itemize
1275
1276 @strong{Note:} when printing expressions, covariant indices and indices
1277 without variance are denoted @samp{.i} while contravariant indices are denoted
1278 @samp{~i}. In the following, we are going to use that notation in the text
1279 so instead of @math{A^i_jk} we will write @samp{A~i.j.k}. Index dimensions
1280 are not visible in the output.
1281
1282 A simple example shall illustrate the concepts:
1283
1284 @example
1285 #include <ginac/ginac.h>
1286 using namespace GiNaC;
1287
1288 int main()
1289 @{
1290     symbol i_sym("i"), j_sym("j");
1291     idx i(i_sym, 3), j(j_sym, 3);
1292
1293     symbol A("A");
1294     cout << indexed(A, i, j) << endl;
1295      // -> A.i.j
1296     ...
1297 @end example
1298
1299 The @code{idx} constructor takes two arguments, the index value and the
1300 index dimension. First we define two index objects, @code{i} and @code{j},
1301 both with the numeric dimension 3. The value of the index @code{i} is the
1302 symbol @code{i_sym} (which prints as @samp{i}) and the value of the index
1303 @code{j} is the symbol @code{j_sym} (which prints as @samp{j}). Next we
1304 construct an expression containing one indexed object, @samp{A.i.j}. It has
1305 the symbol @code{A} as its base expression and the two indices @code{i} and
1306 @code{j}.
1307
1308 Note the difference between the indices @code{i} and @code{j} which are of
1309 class @code{idx}, and the index values which are the sybols @code{i_sym}
1310 and @code{j_sym}. The indices of indexed objects cannot directly be symbols
1311 or numbers but must be index objects. For example, the following is not
1312 correct and will raise an exception:
1313
1314 @example
1315 symbol i("i"), j("j");
1316 e = indexed(A, i, j); // ERROR: indices must be of type idx
1317 @end example
1318
1319 You can have multiple indexed objects in an expression, index values can
1320 be numeric, and index dimensions symbolic:
1321
1322 @example
1323     ...
1324     symbol B("B"), dim("dim");
1325     cout << 4 * indexed(A, i)
1326           + indexed(B, idx(j_sym, 4), idx(2, 3), idx(i_sym, dim)) << endl;
1327      // -> B.j.2.i+4*A.i
1328     ...
1329 @end example
1330
1331 @code{B} has a 4-dimensional symbolic index @samp{k}, a 3-dimensional numeric
1332 index of value 2, and a symbolic index @samp{i} with the symbolic dimension
1333 @samp{dim}. Note that GiNaC doesn't automatically notify you that the free
1334 indices of @samp{A} and @samp{B} in the sum don't match (you have to call
1335 @code{simplify_indexed()} for that, see below).
1336
1337 In fact, base expressions, index values and index dimensions can be
1338 arbitrary expressions:
1339
1340 @example
1341     ...
1342     cout << indexed(A+B, idx(2*i_sym+1, dim/2)) << endl;
1343      // -> (B+A).(1+2*i)
1344     ...
1345 @end example
1346
1347 It's also possible to construct nonsense like @samp{Pi.sin(x)}. You will not
1348 get an error message from this but you will probably not be able to do
1349 anything useful with it.
1350
1351 @cindex @code{get_value()}
1352 @cindex @code{get_dimension()}
1353 The methods
1354
1355 @example
1356 ex idx::get_value(void);
1357 ex idx::get_dimension(void);
1358 @end example
1359
1360 return the value and dimension of an @code{idx} object. If you have an index
1361 in an expression, such as returned by calling @code{.op()} on an indexed
1362 object, you can get a reference to the @code{idx} object with the function
1363 @code{ex_to_idx()} on the expression.
1364
1365 There are also the methods
1366
1367 @example
1368 bool idx::is_numeric(void);
1369 bool idx::is_symbolic(void);
1370 bool idx::is_dim_numeric(void);
1371 bool idx::is_dim_symbolic(void);
1372 @end example
1373
1374 for checking whether the value and dimension are numeric or symbolic
1375 (non-numeric). Using the @code{info()} method of an index (see @ref{Information
1376 About Expressions}) returns information about the index value.
1377
1378 @cindex @code{varidx} (class)
1379 If you need co- and contravariant indices, use the @code{varidx} class:
1380
1381 @example
1382     ...
1383     symbol mu_sym("mu"), nu_sym("nu");
1384     varidx mu(mu_sym, 4), nu(nu_sym, 4); // default is contravariant ~mu, ~nu
1385     varidx mu_co(mu_sym, 4, true);       // covariant index .mu
1386
1387     cout << indexed(A, mu, nu) << endl;
1388      // -> A~mu~nu
1389     cout << indexed(A, mu_co, nu) << endl;
1390      // -> A.mu~nu
1391     cout << indexed(A, mu.toggle_variance(), nu) << endl;
1392      // -> A.mu~nu
1393     ...
1394 @end example
1395
1396 A @code{varidx} is an @code{idx} with an additional flag that marks it as
1397 co- or contravariant. The default is a contravariant (upper) index, but
1398 this can be overridden by supplying a third argument to the @code{varidx}
1399 constructor. The two methods
1400
1401 @example
1402 bool varidx::is_covariant(void);
1403 bool varidx::is_contravariant(void);
1404 @end example
1405
1406 allow you to check the variance of a @code{varidx} object (use @code{ex_to_varidx()}
1407 to get the object reference from an expression). There's also the very useful
1408 method
1409
1410 @example
1411 ex varidx::toggle_variance(void);
1412 @end example
1413
1414 which makes a new index with the same value and dimension but the opposite
1415 variance. By using it you only have to define the index once.
1416
1417 @subsection Substituting indices
1418
1419 @cindex @code{subs()}
1420 Sometimes you will want to substitute one symbolic index with another
1421 symbolic or numeric index, for example when calculating one specific element
1422 of a tensor expression. This is done with the @code{.subs()} method, as it
1423 is done for symbols (see @ref{Substituting Symbols}).
1424
1425 You have two possibilities here. You can either substitute the whole index
1426 by another index or expression:
1427
1428 @example
1429     ...
1430     ex e = indexed(A, mu_co);
1431     cout << e << " becomes " << e.subs(mu_co == nu) << endl;
1432      // -> A.mu becomes A~nu
1433     cout << e << " becomes " << e.subs(mu_co == varidx(0, 4)) << endl;
1434      // -> A.mu becomes A~0
1435     cout << e << " becomes " << e.subs(mu_co == 0) << endl;
1436      // -> A.mu becomes A.0
1437     ...
1438 @end example
1439
1440 The third example shows that trying to replace an index with something that
1441 is not an index will substitute the index value instead.
1442
1443 Alternatively, you can substitute the @emph{symbol} of a symbolic index by
1444 another expression:
1445
1446 @example
1447     ...
1448     ex e = indexed(A, mu_co);
1449     cout << e << " becomes " << e.subs(mu_sym == nu_sym) << endl;
1450      // -> A.mu becomes A.nu
1451     cout << e << " becomes " << e.subs(mu_sym == 0) << endl;
1452      // -> A.mu becomes A.0
1453     ...
1454 @end example
1455
1456 As you see, with the second method only the value of the index will get
1457 substituted. Its other properties, including its dimension, remain unchanged.
1458 If you want to change the dimension of an index you have to substitute the
1459 whole index by another one with the new dimension.
1460
1461 Finally, substituting the base expression of an indexed object works as
1462 expected:
1463
1464 @example
1465     ...
1466     ex e = indexed(A, mu_co);
1467     cout << e << " becomes " << e.subs(A == A+B) << endl;
1468      // -> A.mu becomes (B+A).mu
1469     ...
1470 @end example
1471
1472 @subsection Symmetries
1473
1474 Indexed objects can be declared as being totally symmetric or antisymmetric
1475 with respect to their indices. In this case, GiNaC will automatically bring
1476 the indices into a canonical order which allows for some immediate
1477 simplifications:
1478
1479 @example
1480     ...
1481     cout << indexed(A, indexed::symmetric, i, j)
1482           + indexed(A, indexed::symmetric, j, i) << endl;
1483      // -> 2*A.j.i
1484     cout << indexed(B, indexed::antisymmetric, i, j)
1485           + indexed(B, indexed::antisymmetric, j, j) << endl;
1486      // -> -B.j.i
1487     cout << indexed(B, indexed::antisymmetric, i, j)
1488           + indexed(B, indexed::antisymmetric, j, i) << endl;
1489      // -> 0
1490     ...
1491 @end example
1492
1493 @cindex @code{get_free_indices()}
1494 @subsection Dummy indices
1495
1496 GiNaC treats certain symbolic index pairs as @dfn{dummy indices} meaning
1497 that a summation over the index range is implied. Symbolic indices which are
1498 not dummy indices are called @dfn{free indices}. Numeric indices are neither
1499 dummy nor free indices.
1500
1501 To be recognized as a dummy index pair, the two indices must be of the same
1502 class and dimension and their value must be the same single symbol (an index
1503 like @samp{2*n+1} is never a dummy index). If the indices are of class
1504 @code{varidx}, they must also be of opposite variance.
1505
1506 The method @code{.get_free_indices()} returns a vector containing the free
1507 indices of an expression. It also checks that the free indices of the terms
1508 of a sum are consistent:
1509
1510 @example
1511 @{
1512     symbol A("A"), B("B"), C("C");
1513
1514     symbol i_sym("i"), j_sym("j"), k_sym("k"), l_sym("l");
1515     idx i(i_sym, 3), j(j_sym, 3), k(k_sym, 3), l(l_sym, 3);
1516
1517     ex e = indexed(A, i, j) * indexed(B, j, k) + indexed(C, k, l, i, l);
1518     cout << exprseq(e.get_free_indices()) << endl;
1519      // -> (.i,.k)
1520      // 'j' and 'l' are dummy indices
1521
1522     symbol mu_sym("mu"), nu_sym("nu"), rho_sym("rho"), sigma_sym("sigma");
1523     varidx mu(mu_sym, 4), nu(nu_sym, 4), rho(rho_sym, 4), sigma(sigma_sym, 4);
1524
1525     e = indexed(A, mu, nu) * indexed(B, nu.toggle_variance(), rho)
1526       + indexed(C, mu, sigma, rho, sigma.toggle_variance());
1527     cout << exprseq(e.get_free_indices()) << endl;
1528      // -> (~mu,~rho)
1529      // 'nu' is a dummy index, but 'sigma' is not
1530
1531     e = indexed(A, mu, mu);
1532     cout << exprseq(e.get_free_indices()) << endl;
1533      // -> (~mu)
1534      // 'mu' is not a dummy index because it appears twice with the same
1535      // variance
1536
1537     e = indexed(A, mu, nu) + 42;
1538     cout << exprseq(e.get_free_indices()) << endl; // ERROR
1539      // this will throw an exception:
1540      // "add::get_free_indices: inconsistent indices in sum"
1541 @}
1542 @end example
1543
1544 @cindex @code{simplify_indexed()}
1545 @subsection Simplifying indexed expressions
1546
1547 In addition to the few automatic simplifications that GiNaC performs on
1548 indexed expressions (such as re-ordering the indices of symmetric tensors
1549 and calculating traces and convolutions of matrices and predefined tensors)
1550 there is the method
1551
1552 @example
1553 ex ex::simplify_indexed(void);
1554 ex ex::simplify_indexed(const scalar_products & sp);
1555 @end example
1556
1557 that performs some more expensive operations:
1558
1559 @itemize
1560 @item it checks the consistency of free indices in sums in the same way
1561   @code{get_free_indices()} does
1562 @item it (symbolically) calculates all possible dummy index summations/contractions
1563   with the predefined tensors (this will be explained in more detail in the
1564   next section)
1565 @item as a special case of dummy index summation, it can replace scalar products
1566   of two tensors with a user-defined value
1567 @end itemize
1568
1569 The last point is done with the help of the @code{scalar_products} class
1570 which is used to store scalar products with known values (this is not an
1571 arithmetic class, you just pass it to @code{simplify_indexed()}):
1572
1573 @example
1574 @{
1575     symbol A("A"), B("B"), C("C"), i_sym("i");
1576     idx i(i_sym, 3);
1577
1578     scalar_products sp;
1579     sp.add(A, B, 0); // A and B are orthogonal
1580     sp.add(A, C, 0); // A and C are orthogonal
1581     sp.add(A, A, 4); // A^2 = 4 (A has length 2)
1582
1583     e = indexed(A + B, i) * indexed(A + C, i);
1584     cout << e << endl;
1585      // -> (B+A).i*(A+C).i
1586
1587     cout << e.expand(expand_options::expand_indexed).simplify_indexed(sp)
1588          << endl;
1589      // -> 4+C.i*B.i
1590 @}
1591 @end example
1592
1593 The @code{scalar_products} object @code{sp} acts as a storage for the
1594 scalar products added to it with the @code{.add()} method. This method
1595 takes three arguments: the two expressions of which the scalar product is
1596 taken, and the expression to replace it with. After @code{sp.add(A, B, 0)},
1597 @code{simplify_indexed()} will replace all scalar products of indexed
1598 objects that have the symbols @code{A} and @code{B} as base expressions
1599 with the single value 0. The number, type and dimension of the indices
1600 doesn't matter; @samp{A~mu~nu*B.mu.nu} would also be replaced by 0.
1601
1602 @cindex @code{expand()}
1603 The example above also illustrates a feature of the @code{expand()} method:
1604 if passed the @code{expand_indexed} option it will distribute indices
1605 over sums, so @samp{(A+B).i} becomes @samp{A.i+B.i}.
1606
1607 @cindex @code{tensor} (class)
1608 @subsection Predefined tensors
1609
1610 Some frequently used special tensors such as the delta, epsilon and metric
1611 tensors are predefined in GiNaC. They have special properties when
1612 contracted with other tensor expressions and some of them have constant
1613 matrix representations (they will evaluate to a number when numeric
1614 indices are specified).
1615
1616 @cindex @code{delta_tensor()}
1617 @subsubsection Delta tensor
1618
1619 The delta tensor takes two indices, is symmetric and has the matrix
1620 representation @code{diag(1,1,1,...)}. It is constructed by the function
1621 @code{delta_tensor()}:
1622
1623 @example
1624 @{
1625     symbol A("A"), B("B");
1626
1627     idx i(symbol("i"), 3), j(symbol("j"), 3),
1628         k(symbol("k"), 3), l(symbol("l"), 3);
1629
1630     ex e = indexed(A, i, j) * indexed(B, k, l)
1631          * delta_tensor(i, k) * delta_tensor(j, l) << endl;
1632     cout << e.simplify_indexed() << endl;
1633      // -> B.i.j*A.i.j
1634
1635     cout << delta_tensor(i, i) << endl;
1636      // -> 3
1637 @}
1638 @end example
1639
1640 @cindex @code{metric_tensor()}
1641 @subsubsection General metric tensor
1642
1643 The function @code{metric_tensor()} creates a general metric tensor with
1644 two indices that can be used to raise/lower tensor indices. The metric
1645 tensor is denoted as @samp{g} in the output and if its indices are of
1646 mixed variance it is automatically replaced by a delta tensor:
1647
1648 @example
1649 @{
1650     symbol A("A");
1651
1652     varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
1653
1654     ex e = metric_tensor(mu, nu) * indexed(A, nu.toggle_variance(), rho);
1655     cout << e.simplify_indexed() << endl;
1656      // -> A~mu~rho
1657
1658     e = delta_tensor(mu, nu.toggle_variance()) * metric_tensor(nu, rho);
1659     cout << e.simplify_indexed() << endl;
1660      // -> g~mu~rho
1661
1662     e = metric_tensor(mu.toggle_variance(), nu.toggle_variance())
1663       * metric_tensor(nu, rho);
1664     cout << e.simplify_indexed() << endl;
1665      // -> delta.mu~rho
1666
1667     e = metric_tensor(nu.toggle_variance(), rho.toggle_variance())
1668       * metric_tensor(mu, nu) * (delta_tensor(mu.toggle_variance(), rho)
1669         + indexed(A, mu.toggle_variance(), rho));
1670     cout << e.simplify_indexed() << endl;
1671      // -> 4+A.rho~rho
1672 @}
1673 @end example
1674
1675 @cindex @code{lorentz_g()}
1676 @subsubsection Minkowski metric tensor
1677
1678 The Minkowski metric tensor is a special metric tensor with a constant
1679 matrix representation which is either @code{diag(1, -1, -1, ...)} (negative
1680 signature, the default) or @code{diag(-1, 1, 1, ...)} (positive signature).
1681 It is created with the function @code{lorentz_g()} (although it is output as
1682 @samp{eta}):
1683
1684 @example
1685 @{
1686     varidx mu(symbol("mu"), 4);
1687
1688     e = delta_tensor(varidx(0, 4), mu.toggle_variance())
1689       * lorentz_g(mu, varidx(0, 4));       // negative signature
1690     cout << e.simplify_indexed() << endl;
1691      // -> 1
1692
1693     e = delta_tensor(varidx(0, 4), mu.toggle_variance())
1694       * lorentz_g(mu, varidx(0, 4), true); // positive signature
1695     cout << e.simplify_indexed() << endl;
1696      // -> -1
1697 @}
1698 @end example
1699
1700 @subsubsection Epsilon tensor
1701
1702 The epsilon tensor is totally antisymmetric, its number of indices is equal
1703 to the dimension of the index space (the indices must all be of the same
1704 numeric dimension), and @samp{eps.1.2.3...} (resp. @samp{eps~0~1~2...}) is
1705 defined to be 1. Its behaviour with indices that have a variance also
1706 depends on the signature of the metric. Epsilon tensors are output as
1707 @samp{eps}.
1708
1709 There are three functions defined to create epsilon tensors in 2, 3 and 4
1710 dimensions:
1711
1712 @example
1713 ex epsilon_tensor(const ex & i1, const ex & i2);
1714 ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
1715 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig = false);
1716 @end example
1717
1718 The first two functions create an epsilon tensor in 2 or 3 Euclidean
1719 dimensions, the last function creates an epsilon tensor in a 4-dimensional
1720 Minkowski space (the last @code{bool} argument specifies whether the metric
1721 has negative or positive signature, as in the case of the Minkowski metric
1722 tensor).
1723
1724 @subsection Linear algebra
1725
1726 The @code{matrix} class can be used with indices to do some simple linear
1727 algebra (products of vectors and matrices, traces and scalar products):
1728
1729 @example
1730 @{
1731     idx i(symbol("i"), 2), j(symbol("j"), 2);
1732     symbol x("x"), y("y");
1733
1734     matrix A(2, 2), X(2, 1);
1735     A.set(0, 0, 1); A.set(0, 1, 2);
1736     A.set(1, 0, 3); A.set(1, 1, 4);
1737     X.set(0, 0, x); X.set(1, 0, y);
1738
1739     cout << indexed(A, i, i) << endl;
1740      // -> 5
1741
1742     ex e = indexed(A, i, j) * indexed(X, j);
1743     cout << e.simplify_indexed() << endl;
1744      // -> [[ [[2*y+x]], [[4*y+3*x]] ]].i
1745
1746     e = indexed(A, i, j) * indexed(X, i);
1747     cout << e.simplify_indexed() << endl;
1748      // -> [[ [[3*y+x,4*y+2*x]] ]].j
1749 @}
1750 @end example
1751
1752 You can of course obtain the same results with the @code{matrix::mul()}
1753 and @code{matrix::trace()} methods but with indices you don't have to
1754 worry about transposing matrices.
1755
1756 Matrix indices always start at 0 and their dimension must match the number
1757 of rows/columns of the matrix. Matrices with one row or one column are
1758 vectors and can have one or two indices (it doesn't matter whether it's a
1759 row or a columnt vector). Other matrices must have two indices.
1760
1761 You should be careful when using indices with variance on matrices. GiNaC
1762 doesn't look at the variance and doesn't know that @samp{F~mu~nu} and
1763 @samp{F.mu.nu} are different matrices. In this case you should use only
1764 one form for @samp{F} and explicitly multiply it with a matrix representation
1765 of the metric tensor.
1766
1767
1768 @node Methods and Functions, Information About Expressions, Indexed objects, Top
1769 @c    node-name, next, previous, up
1770 @chapter Methods and Functions
1771 @cindex polynomial
1772
1773 In this chapter the most important algorithms provided by GiNaC will be
1774 described.  Some of them are implemented as functions on expressions,
1775 others are implemented as methods provided by expression objects.  If
1776 they are methods, there exists a wrapper function around it, so you can
1777 alternatively call it in a functional way as shown in the simple
1778 example:
1779
1780 @example
1781 #include <ginac/ginac.h>
1782 using namespace GiNaC;
1783
1784 int main()
1785 @{
1786     ex x = numeric(1.0);
1787     
1788     cout << "As method:   " << sin(x).evalf() << endl;
1789     cout << "As function: " << evalf(sin(x)) << endl;
1790 @}
1791 @end example
1792
1793 @cindex @code{subs()}
1794 The general rule is that wherever methods accept one or more parameters
1795 (@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
1796 wrapper accepts is the same but preceded by the object to act on
1797 (@var{object}, @var{arg1}, @var{arg2}, @dots{}).  This approach is the
1798 most natural one in an OO model but it may lead to confusion for MapleV
1799 users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
1800 would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
1801 @code{A} and @code{x}).  On the other hand, since MapleV returns 3 on
1802 @code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
1803 coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
1804 here.  Also, users of MuPAD will in most cases feel more comfortable
1805 with GiNaC's convention.  All function wrappers are implemented
1806 as simple inline functions which just call the corresponding method and
1807 are only provided for users uncomfortable with OO who are dead set to
1808 avoid method invocations.  Generally, nested function wrappers are much
1809 harder to read than a sequence of methods and should therefore be
1810 avoided if possible.  On the other hand, not everything in GiNaC is a
1811 method on class @code{ex} and sometimes calling a function cannot be
1812 avoided.
1813
1814 @menu
1815 * Information About Expressions::
1816 * Substituting Symbols::
1817 * Polynomial Arithmetic::           Working with polynomials.
1818 * Rational Expressions::            Working with rational functions.
1819 * Symbolic Differentiation::
1820 * Series Expansion::                Taylor and Laurent expansion.
1821 * Built-in Functions::              List of predefined mathematical functions.
1822 * Input/Output::                    Input and output of expressions.
1823 @end menu
1824
1825
1826 @node Information About Expressions, Substituting Symbols, Methods and Functions, Methods and Functions
1827 @c    node-name, next, previous, up
1828 @section Getting information about expressions
1829
1830 @subsection Checking expression types
1831 @cindex @code{is_ex_of_type()}
1832 @cindex @code{ex_to_numeric()}
1833 @cindex @code{ex_to_@dots{}}
1834 @cindex @code{Converting ex to other classes}
1835 @cindex @code{info()}
1836
1837 Sometimes it's useful to check whether a given expression is a plain number,
1838 a sum, a polynomial with integer coefficients, or of some other specific type.
1839 GiNaC provides two functions for this (the first one is actually a macro):
1840
1841 @example
1842 bool is_ex_of_type(const ex & e, TYPENAME t);
1843 bool ex::info(unsigned flag);
1844 @end example
1845
1846 When the test made by @code{is_ex_of_type()} returns true, it is safe to
1847 call one of the functions @code{ex_to_@dots{}}, where @code{@dots{}} is
1848 one of the class names (@xref{The Class Hierarchy}, for a list of all
1849 classes). For example, assuming @code{e} is an @code{ex}:
1850
1851 @example
1852 @{
1853     @dots{}
1854     if (is_ex_of_type(e, numeric))
1855         numeric n = ex_to_numeric(e);
1856     @dots{}
1857 @}
1858 @end example
1859
1860 @code{is_ex_of_type()} allows you to check whether the top-level object of
1861 an expression @samp{e} is an instance of the GiNaC class @samp{t}
1862 (@xref{The Class Hierarchy}, for a list of all classes). This is most useful,
1863 e.g., for checking whether an expression is a number, a sum, or a product:
1864
1865 @example
1866 @{
1867     symbol x("x");
1868     ex e1 = 42;
1869     ex e2 = 4*x - 3;
1870     is_ex_of_type(e1, numeric);  // true
1871     is_ex_of_type(e2, numeric);  // false
1872     is_ex_of_type(e1, add);      // false
1873     is_ex_of_type(e2, add);      // true
1874     is_ex_of_type(e1, mul);      // false
1875     is_ex_of_type(e2, mul);      // false
1876 @}
1877 @end example
1878
1879 The @code{info()} method is used for checking certain attributes of
1880 expressions. The possible values for the @code{flag} argument are defined
1881 in @file{ginac/flags.h}, the most important being explained in the following
1882 table:
1883
1884 @cartouche
1885 @multitable @columnfractions .30 .70
1886 @item @strong{Flag} @tab @strong{Returns true if the object is@dots{}}
1887 @item @code{numeric}
1888 @tab @dots{}a number (same as @code{is_ex_of_type(..., numeric)})
1889 @item @code{real}
1890 @tab @dots{}a real integer, rational or float (i.e. is not complex)
1891 @item @code{rational}
1892 @tab @dots{}an exact rational number (integers are rational, too)
1893 @item @code{integer}
1894 @tab @dots{}a (non-complex) integer
1895 @item @code{crational}
1896 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
1897 @item @code{cinteger}
1898 @tab @dots{}a (complex) integer (such as @math{2-3*I})
1899 @item @code{positive}
1900 @tab @dots{}not complex and greater than 0
1901 @item @code{negative}
1902 @tab @dots{}not complex and less than 0
1903 @item @code{nonnegative}
1904 @tab @dots{}not complex and greater than or equal to 0
1905 @item @code{posint}
1906 @tab @dots{}an integer greater than 0
1907 @item @code{negint}
1908 @tab @dots{}an integer less than 0
1909 @item @code{nonnegint}
1910 @tab @dots{}an integer greater than or equal to 0
1911 @item @code{even}
1912 @tab @dots{}an even integer
1913 @item @code{odd}
1914 @tab @dots{}an odd integer
1915 @item @code{prime}
1916 @tab @dots{}a prime integer (probabilistic primality test)
1917 @item @code{relation}
1918 @tab @dots{}a relation (same as @code{is_ex_of_type(..., relational)})
1919 @item @code{relation_equal}
1920 @tab @dots{}a @code{==} relation
1921 @item @code{relation_not_equal}
1922 @tab @dots{}a @code{!=} relation
1923 @item @code{relation_less}
1924 @tab @dots{}a @code{<} relation
1925 @item @code{relation_less_or_equal}
1926 @tab @dots{}a @code{<=} relation
1927 @item @code{relation_greater}
1928 @tab @dots{}a @code{>} relation
1929 @item @code{relation_greater_or_equal}
1930 @tab @dots{}a @code{>=} relation
1931 @item @code{symbol}
1932 @tab @dots{}a symbol (same as @code{is_ex_of_type(..., symbol)})
1933 @item @code{list}
1934 @tab @dots{}a list (same as @code{is_ex_of_type(..., lst)})
1935 @item @code{polynomial}
1936 @tab @dots{}a polynomial (i.e. only consists of sums and products of numbers and symbols with positive integer powers)
1937 @item @code{integer_polynomial}
1938 @tab @dots{}a polynomial with (non-complex) integer coefficients
1939 @item @code{cinteger_polynomial}
1940 @tab @dots{}a polynomial with (possibly complex) integer coefficients (such as @math{2-3*I})
1941 @item @code{rational_polynomial}
1942 @tab @dots{}a polynomial with (non-complex) rational coefficients
1943 @item @code{crational_polynomial}
1944 @tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I})
1945 @item @code{rational_function}
1946 @tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)})
1947 @item @code{algebraic}
1948 @tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1})
1949 @end multitable
1950 @end cartouche
1951
1952
1953 @subsection Accessing subexpressions
1954 @cindex @code{nops()}
1955 @cindex @code{op()}
1956 @cindex @code{has()}
1957 @cindex container
1958 @cindex @code{relational} (class)
1959
1960 GiNaC provides the two methods
1961
1962 @example
1963 unsigned ex::nops();
1964 ex ex::op(unsigned i);
1965 @end example
1966
1967 for accessing the subexpressions in the container-like GiNaC classes like
1968 @code{add}, @code{mul}, @code{lst}, and @code{function}. @code{nops()}
1969 determines the number of subexpressions (@samp{operands}) contained, while
1970 @code{op()} returns the @code{i}-th (0..@code{nops()-1}) subexpression.
1971 In the case of a @code{power} object, @code{op(0)} will return the basis
1972 and @code{op(1)} the exponent. For @code{indexed} objects, @code{op(0)}
1973 is the base expression and @code{op(i)}, @math{i>0} are the indices.
1974
1975 The left-hand and right-hand side expressions of objects of class
1976 @code{relational} (and only of these) can also be accessed with the methods
1977
1978 @example
1979 ex ex::lhs();
1980 ex ex::rhs();
1981 @end example
1982
1983 Finally, the method
1984
1985 @example
1986 bool ex::has(const ex & other);
1987 @end example
1988
1989 checks whether an expression contains the given subexpression @code{other}.
1990 This only works reliably if @code{other} is of an atomic class such as a
1991 @code{numeric} or a @code{symbol}. It is, e.g., not possible to verify that
1992 @code{a+b+c} contains @code{a+c} (or @code{a+b}) as a subexpression.
1993
1994
1995 @subsection Comparing expressions
1996 @cindex @code{is_equal()}
1997 @cindex @code{is_zero()}
1998
1999 Expressions can be compared with the usual C++ relational operators like
2000 @code{==}, @code{>}, and @code{<} but if the expressions contain symbols,
2001 the result is usually not determinable and the result will be @code{false},
2002 except in the case of the @code{!=} operator. You should also be aware that
2003 GiNaC will only do the most trivial test for equality (subtracting both
2004 expressions), so something like @code{(pow(x,2)+x)/x==x+1} will return
2005 @code{false}.
2006
2007 Actually, if you construct an expression like @code{a == b}, this will be
2008 represented by an object of the @code{relational} class (@xref{Relations}.)
2009 which is not evaluated until (explicitly or implicitely) cast to a @code{bool}.
2010
2011 There are also two methods
2012
2013 @example
2014 bool ex::is_equal(const ex & other);
2015 bool ex::is_zero();
2016 @end example
2017
2018 for checking whether one expression is equal to another, or equal to zero,
2019 respectively.
2020
2021 @strong{Warning:} You will also find an @code{ex::compare()} method in the
2022 GiNaC header files. This method is however only to be used internally by
2023 GiNaC to establish a canonical sort order for terms, and using it to compare
2024 expressions will give very surprising results.
2025
2026
2027 @node Substituting Symbols, Polynomial Arithmetic, Information About Expressions, Methods and Functions
2028 @c    node-name, next, previous, up
2029 @section Substituting symbols
2030 @cindex @code{subs()}
2031
2032 Symbols can be replaced with expressions via the @code{.subs()} method:
2033
2034 @example
2035 ex ex::subs(const ex & e);
2036 ex ex::subs(const lst & syms, const lst & repls);
2037 @end example
2038
2039 In the first form, @code{subs()} accepts a relational of the form
2040 @samp{symbol == expression} or a @code{lst} of such relationals. E.g.
2041
2042 @example
2043 @{
2044     symbol x("x"), y("y");
2045     ex e1 = 2*x^2-4*x+3;
2046     cout << "e1(7) = " << e1.subs(x == 7) << endl;
2047     ex e2 = x*y + x;
2048     cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
2049 @}
2050 @end example
2051
2052 will print @samp{73} and @samp{-10}, respectively.
2053
2054 If you specify multiple substitutions, they are performed in parallel, so e.g.
2055 @code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
2056
2057 The second form of @code{subs()} takes two lists, one for the symbols and
2058 one for the expressions to be substituted (both lists must contain the same
2059 number of elements). Using this form, you would write @code{subs(lst(x, y), lst(y, x))}
2060 to exchange @samp{x} and @samp{y}.
2061
2062
2063 @node Polynomial Arithmetic, Rational Expressions, Substituting Symbols, Methods and Functions
2064 @c    node-name, next, previous, up
2065 @section Polynomial arithmetic
2066
2067 @subsection Expanding and collecting
2068 @cindex @code{expand()}
2069 @cindex @code{collect()}
2070
2071 A polynomial in one or more variables has many equivalent
2072 representations.  Some useful ones serve a specific purpose.  Consider
2073 for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
2074 21*y*z + 4*z^2} (written down here in output-style).  It is equivalent
2075 to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}.  Other
2076 representations are the recursive ones where one collects for exponents
2077 in one of the three variable.  Since the factors are themselves
2078 polynomials in the remaining two variables the procedure can be
2079 repeated.  In our expample, two possibilities would be @math{(4*y + z)*x
2080 + 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
2081 x*z}.
2082
2083 To bring an expression into expanded form, its method
2084
2085 @example
2086 ex ex::expand();
2087 @end example
2088
2089 may be called.  In our example above, this corresponds to @math{4*x*y +
2090 x*z + 20*y^2 + 21*y*z + 4*z^2}.  Again, since the canonical form in
2091 GiNaC is not easily guessable you should be prepared to see different
2092 orderings of terms in such sums!
2093
2094 Another useful representation of multivariate polynomials is as a
2095 univariate polynomial in one of the variables with the coefficients
2096 being polynomials in the remaining variables.  The method
2097 @code{collect()} accomplishes this task:
2098
2099 @example
2100 ex ex::collect(const symbol & s);
2101 @end example
2102
2103 Note that the original polynomial needs to be in expanded form in order
2104 to be able to find the coefficients properly.
2105
2106 @subsection Degree and coefficients
2107 @cindex @code{degree()}
2108 @cindex @code{ldegree()}
2109 @cindex @code{coeff()}
2110
2111 The degree and low degree of a polynomial can be obtained using the two
2112 methods
2113
2114 @example
2115 int ex::degree(const symbol & s);
2116 int ex::ldegree(const symbol & s);
2117 @end example
2118
2119 which also work reliably on non-expanded input polynomials (they even work
2120 on rational functions, returning the asymptotic degree). To extract
2121 a coefficient with a certain power from an expanded polynomial you use
2122
2123 @example
2124 ex ex::coeff(const symbol & s, int n);
2125 @end example
2126
2127 You can also obtain the leading and trailing coefficients with the methods
2128
2129 @example
2130 ex ex::lcoeff(const symbol & s);
2131 ex ex::tcoeff(const symbol & s);
2132 @end example
2133
2134 which are equivalent to @code{coeff(s, degree(s))} and @code{coeff(s, ldegree(s))},
2135 respectively.
2136
2137 An application is illustrated in the next example, where a multivariate
2138 polynomial is analyzed:
2139
2140 @example
2141 #include <ginac/ginac.h>
2142 using namespace GiNaC;
2143
2144 int main()
2145 @{
2146     symbol x("x"), y("y");
2147     ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
2148                  - pow(x+y,2) + 2*pow(y+2,2) - 8;
2149     ex Poly = PolyInp.expand();
2150     
2151     for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) @{
2152         cout << "The x^" << i << "-coefficient is "
2153              << Poly.coeff(x,i) << endl;
2154     @}
2155     cout << "As polynomial in y: " 
2156          << Poly.collect(y) << endl;
2157 @}
2158 @end example
2159
2160 When run, it returns an output in the following fashion:
2161
2162 @example
2163 The x^0-coefficient is y^2+11*y
2164 The x^1-coefficient is 5*y^2-2*y
2165 The x^2-coefficient is -1
2166 The x^3-coefficient is 4*y
2167 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
2168 @end example
2169
2170 As always, the exact output may vary between different versions of GiNaC
2171 or even from run to run since the internal canonical ordering is not
2172 within the user's sphere of influence.
2173
2174
2175 @subsection Polynomial division
2176 @cindex polynomial division
2177 @cindex quotient
2178 @cindex remainder
2179 @cindex pseudo-remainder
2180 @cindex @code{quo()}
2181 @cindex @code{rem()}
2182 @cindex @code{prem()}
2183 @cindex @code{divide()}
2184
2185 The two functions
2186
2187 @example
2188 ex quo(const ex & a, const ex & b, const symbol & x);
2189 ex rem(const ex & a, const ex & b, const symbol & x);
2190 @end example
2191
2192 compute the quotient and remainder of univariate polynomials in the variable
2193 @samp{x}. The results satisfy @math{a = b*quo(a, b, x) + rem(a, b, x)}.
2194
2195 The additional function
2196
2197 @example
2198 ex prem(const ex & a, const ex & b, const symbol & x);
2199 @end example
2200
2201 computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
2202 @math{c*a = b*q + prem(a, b, x)}, where @math{c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1)}.
2203
2204 Exact division of multivariate polynomials is performed by the function
2205
2206 @example
2207 bool divide(const ex & a, const ex & b, ex & q);
2208 @end example
2209
2210 If @samp{b} divides @samp{a} over the rationals, this function returns @code{true}
2211 and returns the quotient in the variable @code{q}. Otherwise it returns @code{false}
2212 in which case the value of @code{q} is undefined.
2213
2214
2215 @subsection Unit, content and primitive part
2216 @cindex @code{unit()}
2217 @cindex @code{content()}
2218 @cindex @code{primpart()}
2219
2220 The methods
2221
2222 @example
2223 ex ex::unit(const symbol & x);
2224 ex ex::content(const symbol & x);
2225 ex ex::primpart(const symbol & x);
2226 @end example
2227
2228 return the unit part, content part, and primitive polynomial of a multivariate
2229 polynomial with respect to the variable @samp{x} (the unit part being the sign
2230 of the leading coefficient, the content part being the GCD of the coefficients,
2231 and the primitive polynomial being the input polynomial divided by the unit and
2232 content parts). The product of unit, content, and primitive part is the
2233 original polynomial.
2234
2235
2236 @subsection GCD and LCM
2237 @cindex GCD
2238 @cindex LCM
2239 @cindex @code{gcd()}
2240 @cindex @code{lcm()}
2241
2242 The functions for polynomial greatest common divisor and least common
2243 multiple have the synopsis
2244
2245 @example
2246 ex gcd(const ex & a, const ex & b);
2247 ex lcm(const ex & a, const ex & b);
2248 @end example
2249
2250 The functions @code{gcd()} and @code{lcm()} accept two expressions
2251 @code{a} and @code{b} as arguments and return a new expression, their
2252 greatest common divisor or least common multiple, respectively.  If the
2253 polynomials @code{a} and @code{b} are coprime @code{gcd(a,b)} returns 1
2254 and @code{lcm(a,b)} returns the product of @code{a} and @code{b}.
2255
2256 @example
2257 #include <ginac/ginac.h>
2258 using namespace GiNaC;
2259
2260 int main()
2261 @{
2262     symbol x("x"), y("y"), z("z");
2263     ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
2264     ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
2265
2266     ex P_gcd = gcd(P_a, P_b);
2267     // x + 5*y + 4*z
2268     ex P_lcm = lcm(P_a, P_b);
2269     // 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
2270 @}
2271 @end example
2272
2273
2274 @node Rational Expressions, Symbolic Differentiation, Polynomial Arithmetic, Methods and Functions
2275 @c    node-name, next, previous, up
2276 @section Rational expressions
2277
2278 @subsection The @code{normal} method
2279 @cindex @code{normal()}
2280 @cindex simplification
2281 @cindex temporary replacement
2282
2283 Some basic form of simplification of expressions is called for frequently.
2284 GiNaC provides the method @code{.normal()}, which converts a rational function
2285 into an equivalent rational function of the form @samp{numerator/denominator}
2286 where numerator and denominator are coprime.  If the input expression is already
2287 a fraction, it just finds the GCD of numerator and denominator and cancels it,
2288 otherwise it performs fraction addition and multiplication.
2289
2290 @code{.normal()} can also be used on expressions which are not rational functions
2291 as it will replace all non-rational objects (like functions or non-integer
2292 powers) by temporary symbols to bring the expression to the domain of rational
2293 functions before performing the normalization, and re-substituting these
2294 symbols afterwards. This algorithm is also available as a separate method
2295 @code{.to_rational()}, described below.
2296
2297 This means that both expressions @code{t1} and @code{t2} are indeed
2298 simplified in this little program:
2299
2300 @example
2301 #include <ginac/ginac.h>
2302 using namespace GiNaC;
2303
2304 int main()
2305 @{
2306     symbol x("x");
2307     ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
2308     ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
2309     cout << "t1 is " << t1.normal() << endl;
2310     cout << "t2 is " << t2.normal() << endl;
2311 @}
2312 @end example
2313
2314 Of course this works for multivariate polynomials too, so the ratio of
2315 the sample-polynomials from the section about GCD and LCM above would be
2316 normalized to @code{P_a/P_b} = @code{(4*y+z)/(y+3*z)}.
2317
2318
2319 @subsection Numerator and denominator
2320 @cindex numerator
2321 @cindex denominator
2322 @cindex @code{numer()}
2323 @cindex @code{denom()}
2324
2325 The numerator and denominator of an expression can be obtained with
2326
2327 @example
2328 ex ex::numer();
2329 ex ex::denom();
2330 @end example
2331
2332 These functions will first normalize the expression as described above and
2333 then return the numerator or denominator, respectively.
2334
2335
2336 @subsection Converting to a rational expression
2337 @cindex @code{to_rational()}
2338
2339 Some of the methods described so far only work on polynomials or rational
2340 functions. GiNaC provides a way to extend the domain of these functions to
2341 general expressions by using the temporary replacement algorithm described
2342 above. You do this by calling
2343
2344 @example
2345 ex ex::to_rational(lst &l);
2346 @end example
2347
2348 on the expression to be converted. The supplied @code{lst} will be filled
2349 with the generated temporary symbols and their replacement expressions in
2350 a format that can be used directly for the @code{subs()} method. It can also
2351 already contain a list of replacements from an earlier application of
2352 @code{.to_rational()}, so it's possible to use it on multiple expressions
2353 and get consistent results.
2354
2355 For example,
2356
2357 @example
2358 @{
2359     symbol x("x");
2360     ex a = pow(sin(x), 2) - pow(cos(x), 2);
2361     ex b = sin(x) + cos(x);
2362     ex q;
2363     lst l;
2364     divide(a.to_rational(l), b.to_rational(l), q);
2365     cout << q.subs(l) << endl;
2366 @}
2367 @end example
2368
2369 will print @samp{sin(x)-cos(x)}.
2370
2371
2372 @node Symbolic Differentiation, Series Expansion, Rational Expressions, Methods and Functions
2373 @c    node-name, next, previous, up
2374 @section Symbolic differentiation
2375 @cindex differentiation
2376 @cindex @code{diff()}
2377 @cindex chain rule
2378 @cindex product rule
2379
2380 GiNaC's objects know how to differentiate themselves.  Thus, a
2381 polynomial (class @code{add}) knows that its derivative is the sum of
2382 the derivatives of all the monomials:
2383
2384 @example
2385 #include <ginac/ginac.h>
2386 using namespace GiNaC;
2387
2388 int main()
2389 @{
2390     symbol x("x"), y("y"), z("z");
2391     ex P = pow(x, 5) + pow(x, 2) + y;
2392
2393     cout << P.diff(x,2) << endl;  // 20*x^3 + 2
2394     cout << P.diff(y) << endl;    // 1
2395     cout << P.diff(z) << endl;    // 0
2396 @}
2397 @end example
2398
2399 If a second integer parameter @var{n} is given, the @code{diff} method
2400 returns the @var{n}th derivative.
2401
2402 If @emph{every} object and every function is told what its derivative
2403 is, all derivatives of composed objects can be calculated using the
2404 chain rule and the product rule.  Consider, for instance the expression
2405 @code{1/cosh(x)}.  Since the derivative of @code{cosh(x)} is
2406 @code{sinh(x)} and the derivative of @code{pow(x,-1)} is
2407 @code{-pow(x,-2)}, GiNaC can readily compute the composition.  It turns
2408 out that the composition is the generating function for Euler Numbers,
2409 i.e. the so called @var{n}th Euler number is the coefficient of
2410 @code{x^n/n!} in the expansion of @code{1/cosh(x)}.  We may use this
2411 identity to code a function that generates Euler numbers in just three
2412 lines:
2413
2414 @cindex Euler numbers
2415 @example
2416 #include <ginac/ginac.h>
2417 using namespace GiNaC;
2418
2419 ex EulerNumber(unsigned n)
2420 @{
2421     symbol x;
2422     const ex generator = pow(cosh(x),-1);
2423     return generator.diff(x,n).subs(x==0);
2424 @}
2425
2426 int main()
2427 @{
2428     for (unsigned i=0; i<11; i+=2)
2429         cout << EulerNumber(i) << endl;
2430     return 0;
2431 @}
2432 @end example
2433
2434 When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
2435 @code{-61}, @code{1385}, @code{-50521}.  We increment the loop variable
2436 @code{i} by two since all odd Euler numbers vanish anyways.
2437
2438
2439 @node Series Expansion, Built-in Functions, Symbolic Differentiation, Methods and Functions
2440 @c    node-name, next, previous, up
2441 @section Series expansion
2442 @cindex @code{series()}
2443 @cindex Taylor expansion
2444 @cindex Laurent expansion
2445 @cindex @code{pseries} (class)
2446
2447 Expressions know how to expand themselves as a Taylor series or (more
2448 generally) a Laurent series.  As in most conventional Computer Algebra
2449 Systems, no distinction is made between those two.  There is a class of
2450 its own for storing such series (@code{class pseries}) and a built-in
2451 function (called @code{Order}) for storing the order term of the series.
2452 As a consequence, if you want to work with series, i.e. multiply two
2453 series, you need to call the method @code{ex::series} again to convert
2454 it to a series object with the usual structure (expansion plus order
2455 term).  A sample application from special relativity could read:
2456
2457 @example
2458 #include <ginac/ginac.h>
2459 using namespace GiNaC;
2460
2461 int main()
2462 @{
2463     symbol v("v"), c("c");
2464     
2465     ex gamma = 1/sqrt(1 - pow(v/c,2));
2466     ex mass_nonrel = gamma.series(v==0, 10);
2467     
2468     cout << "the relativistic mass increase with v is " << endl
2469          << mass_nonrel << endl;
2470     
2471     cout << "the inverse square of this series is " << endl
2472          << pow(mass_nonrel,-2).series(v==0, 10) << endl;
2473 @}
2474 @end example
2475
2476 Only calling the series method makes the last output simplify to
2477 @math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
2478 series raised to the power @math{-2}.
2479
2480 @cindex M@'echain's formula
2481 As another instructive application, let us calculate the numerical 
2482 value of Archimedes' constant
2483 @tex
2484 $\pi$
2485 @end tex
2486 (for which there already exists the built-in constant @code{Pi}) 
2487 using M@'echain's amazing formula
2488 @tex
2489 $\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
2490 @end tex
2491 @ifnottex
2492 @math{Pi==16*atan(1/5)-4*atan(1/239)}.
2493 @end ifnottex
2494 We may expand the arcus tangent around @code{0} and insert the fractions
2495 @code{1/5} and @code{1/239}.  But, as we have seen, a series in GiNaC
2496 carries an order term with it and the question arises what the system is
2497 supposed to do when the fractions are plugged into that order term.  The
2498 solution is to use the function @code{series_to_poly()} to simply strip
2499 the order term off:
2500
2501 @example
2502 #include <ginac/ginac.h>
2503 using namespace GiNaC;
2504
2505 ex mechain_pi(int degr)
2506 @{
2507     symbol x;
2508     ex pi_expansion = series_to_poly(atan(x).series(x,degr));
2509     ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
2510                    -4*pi_expansion.subs(x==numeric(1,239));
2511     return pi_approx;
2512 @}
2513
2514 int main()
2515 @{
2516     ex pi_frac;
2517     for (int i=2; i<12; i+=2) @{
2518         pi_frac = mechain_pi(i);
2519         cout << i << ":\t" << pi_frac << endl
2520              << "\t" << pi_frac.evalf() << endl;
2521     @}
2522     return 0;
2523 @}
2524 @end example
2525
2526 Note how we just called @code{.series(x,degr)} instead of
2527 @code{.series(x==0,degr)}.  This is a simple shortcut for @code{ex}'s
2528 method @code{series()}: if the first argument is a symbol the expression
2529 is expanded in that symbol around point @code{0}.  When you run this
2530 program, it will type out:
2531
2532 @example
2533 2:      3804/1195
2534         3.1832635983263598326
2535 4:      5359397032/1706489875
2536         3.1405970293260603143
2537 6:      38279241713339684/12184551018734375
2538         3.141621029325034425
2539 8:      76528487109180192540976/24359780855939418203125
2540         3.141591772182177295
2541 10:     327853873402258685803048818236/104359128170408663038552734375
2542         3.1415926824043995174
2543 @end example
2544
2545
2546 @node Built-in Functions, Input/Output, Series Expansion, Methods and Functions
2547 @c    node-name, next, previous, up
2548 @section Predefined mathematical functions
2549
2550 GiNaC contains the following predefined mathematical functions:
2551
2552 @cartouche
2553 @multitable @columnfractions .30 .70
2554 @item @strong{Name} @tab @strong{Function}
2555 @item @code{abs(x)}
2556 @tab absolute value
2557 @item @code{csgn(x)}
2558 @tab complex sign
2559 @item @code{sqrt(x)}
2560 @tab square root (not a GiNaC function proper but equivalent to @code{pow(x, numeric(1, 2)})
2561 @item @code{sin(x)}
2562 @tab sine
2563 @item @code{cos(x)}
2564 @tab cosine
2565 @item @code{tan(x)}
2566 @tab tangent
2567 @item @code{asin(x)}
2568 @tab inverse sine
2569 @item @code{acos(x)}
2570 @tab inverse cosine
2571 @item @code{atan(x)}
2572 @tab inverse tangent
2573 @item @code{atan2(y, x)}
2574 @tab inverse tangent with two arguments
2575 @item @code{sinh(x)}
2576 @tab hyperbolic sine
2577 @item @code{cosh(x)}
2578 @tab hyperbolic cosine
2579 @item @code{tanh(x)}
2580 @tab hyperbolic tangent
2581 @item @code{asinh(x)}
2582 @tab inverse hyperbolic sine
2583 @item @code{acosh(x)}
2584 @tab inverse hyperbolic cosine
2585 @item @code{atanh(x)}
2586 @tab inverse hyperbolic tangent
2587 @item @code{exp(x)}
2588 @tab exponential function
2589 @item @code{log(x)}
2590 @tab natural logarithm
2591 @item @code{Li2(x)}
2592 @tab Dilogarithm
2593 @item @code{zeta(x)}
2594 @tab Riemann's zeta function
2595 @item @code{zeta(n, x)}
2596 @tab derivatives of Riemann's zeta function
2597 @item @code{tgamma(x)}
2598 @tab Gamma function
2599 @item @code{lgamma(x)}
2600 @tab logarithm of Gamma function
2601 @item @code{beta(x, y)}
2602 @tab Beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
2603 @item @code{psi(x)}
2604 @tab psi (digamma) function
2605 @item @code{psi(n, x)}
2606 @tab derivatives of psi function (polygamma functions)
2607 @item @code{factorial(n)}
2608 @tab factorial function
2609 @item @code{binomial(n, m)}
2610 @tab binomial coefficients
2611 @item @code{Order(x)}
2612 @tab order term function in truncated power series
2613 @item @code{Derivative(x, l)}
2614 @tab inert partial differentiation operator (used internally)
2615 @end multitable
2616 @end cartouche
2617
2618 @cindex branch cut
2619 For functions that have a branch cut in the complex plane GiNaC follows
2620 the conventions for C++ as defined in the ANSI standard as far as
2621 possible.  In particular: the natural logarithm (@code{log}) and the
2622 square root (@code{sqrt}) both have their branch cuts running along the
2623 negative real axis where the points on the axis itself belong to the
2624 upper part (i.e. continuous with quadrant II).  The inverse
2625 trigonometric and hyperbolic functions are not defined for complex
2626 arguments by the C++ standard, however.  In GiNaC we follow the
2627 conventions used by CLN, which in turn follow the carefully designed
2628 definitions in the Common Lisp standard.  It should be noted that this
2629 convention is identical to the one used by the C99 standard and by most
2630 serious CAS.  It is to be expected that future revisions of the C++
2631 standard incorporate these functions in the complex domain in a manner
2632 compatible with C99.
2633
2634
2635 @node Input/Output, Extending GiNaC, Built-in Functions, Methods and Functions
2636 @c    node-name, next, previous, up
2637 @section Input and output of expressions
2638 @cindex I/O
2639
2640 @subsection Expression output
2641 @cindex printing
2642 @cindex output of expressions
2643
2644 The easiest way to print an expression is to write it to a stream:
2645
2646 @example
2647 @{
2648     symbol x("x");
2649     ex e = 4.5+pow(x,2)*3/2;
2650     cout << e << endl;    // prints '4.5+3/2*x^2'
2651     // ...
2652 @end example
2653
2654 The output format is identical to the @command{ginsh} input syntax and
2655 to that used by most computer algebra systems, but not directly pastable
2656 into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
2657 is printed as @samp{x^2}).
2658
2659 To print an expression in a way that can be directly used in a C or C++
2660 program, you use the method
2661
2662 @example
2663 void ex::printcsrc(ostream & os, unsigned type, const char *name);
2664 @end example
2665
2666 This outputs a line in the form of a variable definition @code{<type> <name> = <expression>}.
2667 The possible types are defined in @file{ginac/flags.h} (@code{csrc_types})
2668 and mostly affect the way in which floating point numbers are written:
2669
2670 @example
2671     // ...
2672     e.printcsrc(cout, csrc_types::ctype_float, "f");
2673     e.printcsrc(cout, csrc_types::ctype_double, "d");
2674     e.printcsrc(cout, csrc_types::ctype_cl_N, "n");
2675     // ...
2676 @end example
2677
2678 The above example will produce (note the @code{x^2} being converted to @code{x*x}):
2679
2680 @example
2681 float f = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
2682 double d = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
2683 cl_N n = (cl_F("3.0")/cl_F("2.0"))*(x*x)+cl_F("4.5");
2684 @end example
2685
2686 Finally, there are the two methods @code{printraw()} and @code{printtree()} intended for GiNaC
2687 developers, that provide a dump of the internal structure of an expression for
2688 debugging purposes:
2689
2690 @example
2691     // ...
2692     e.printraw(cout); cout << endl << endl;
2693     e.printtree(cout);
2694 @}
2695 @end example
2696
2697 produces
2698
2699 @example
2700 ex(+((power(ex(symbol(name=x,serial=1,hash=150875740,flags=11)),ex(numeric(2)),hash=2,flags=3),numeric(3/2)),,hash=0,flags=3))
2701
2702 type=Q25GiNaC3add, hash=0 (0x0), flags=3, nops=2
2703     power: hash=2 (0x2), flags=3
2704         x (symbol): serial=1, hash=150875740 (0x8fe2e5c), flags=11
2705         2 (numeric): hash=2147483714 (0x80000042), flags=11
2706     3/2 (numeric): hash=2147483745 (0x80000061), flags=11
2707     -----
2708     overall_coeff
2709     4.5L0 (numeric): hash=2147483723 (0x8000004b), flags=11
2710     =====
2711 @end example
2712
2713 The @code{printtree()} method is also available in @command{ginsh} as the
2714 @code{print()} function.
2715
2716
2717 @subsection Expression input
2718 @cindex input of expressions
2719
2720 GiNaC provides no way to directly read an expression from a stream because
2721 you will usually want the user to be able to enter something like @samp{2*x+sin(y)}
2722 and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
2723 @code{y} you defined in your program and there is no way to specify the
2724 desired symbols to the @code{>>} stream input operator.
2725
2726 Instead, GiNaC lets you construct an expression from a string, specifying the
2727 list of symbols to be used:
2728
2729 @example
2730 @{
2731     symbol x("x"), y("y");
2732     ex e("2*x+sin(y)", lst(x, y));
2733 @}
2734 @end example
2735
2736 The input syntax is the same as that used by @command{ginsh} and the stream
2737 output operator @code{<<}. The symbols in the string are matched by name to
2738 the symbols in the list and if GiNaC encounters a symbol not specified in
2739 the list it will throw an exception.
2740
2741 With this constructor, it's also easy to implement interactive GiNaC programs:
2742
2743 @example
2744 #include <iostream>
2745 #include <string>
2746 #include <stdexcept>
2747 #include <ginac/ginac.h>
2748 using namespace GiNaC;
2749
2750 int main()
2751 @{
2752      symbol x("x");
2753      string s;
2754
2755      cout << "Enter an expression containing 'x': ";
2756      getline(cin, s);
2757
2758      try @{
2759          ex e(s, lst(x));
2760          cout << "The derivative of " << e << " with respect to x is ";
2761          cout << e.diff(x) << ".\n";
2762      @} catch (exception &p) @{
2763          cerr << p.what() << endl;
2764      @}
2765 @}
2766 @end example
2767
2768
2769 @subsection Archiving
2770 @cindex @code{archive} (class)
2771 @cindex archiving
2772
2773 GiNaC allows creating @dfn{archives} of expressions which can be stored
2774 to or retrieved from files. To create an archive, you declare an object
2775 of class @code{archive} and archive expressions in it, giving each
2776 expression a unique name:
2777
2778 @example
2779 #include <ginac/ginac.h>
2780 #include <fstream>
2781 using namespace GiNaC;
2782
2783 int main()
2784 @{
2785     symbol x("x"), y("y"), z("z");
2786
2787     ex foo = sin(x + 2*y) + 3*z + 41;
2788     ex bar = foo + 1;
2789
2790     archive a;
2791     a.archive_ex(foo, "foo");
2792     a.archive_ex(bar, "the second one");
2793     // ...
2794 @end example
2795
2796 The archive can then be written to a file:
2797
2798 @example
2799     // ...
2800     ofstream out("foobar.gar");
2801     out << a;
2802     out.close();
2803     // ...
2804 @end example
2805
2806 The file @file{foobar.gar} contains all information that is needed to
2807 reconstruct the expressions @code{foo} and @code{bar}.
2808
2809 @cindex @command{viewgar}
2810 The tool @command{viewgar} that comes with GiNaC can be used to view
2811 the contents of GiNaC archive files:
2812
2813 @example
2814 $ viewgar foobar.gar
2815 foo = 41+sin(x+2*y)+3*z
2816 the second one = 42+sin(x+2*y)+3*z
2817 @end example
2818
2819 The point of writing archive files is of course that they can later be
2820 read in again:
2821
2822 @example
2823     // ...
2824     archive a2;
2825     ifstream in("foobar.gar");
2826     in >> a2;
2827     // ...
2828 @end example
2829
2830 And the stored expressions can be retrieved by their name:
2831
2832 @example
2833     // ...
2834     lst syms(x, y);
2835
2836     ex ex1 = a2.unarchive_ex(syms, "foo");
2837     ex ex2 = a2.unarchive_ex(syms, "the second one");
2838
2839     cout << ex1 << endl;              // prints "41+sin(x+2*y)+3*z"
2840     cout << ex2 << endl;              // prints "42+sin(x+2*y)+3*z"
2841     cout << ex1.subs(x == 2) << endl; // prints "41+sin(2+2*y)+3*z"
2842 @}
2843 @end example
2844
2845 Note that you have to supply a list of the symbols which are to be inserted
2846 in the expressions. Symbols in archives are stored by their name only and
2847 if you don't specify which symbols you have, unarchiving the expression will
2848 create new symbols with that name. E.g. if you hadn't included @code{x} in
2849 the @code{syms} list above, the @code{ex1.subs(x == 2)} statement would
2850 have had no effect because the @code{x} in @code{ex1} would have been a
2851 different symbol than the @code{x} which was defined at the beginning of
2852 the program, altough both would appear as @samp{x} when printed.
2853
2854
2855
2856 @node Extending GiNaC, What does not belong into GiNaC, Input/Output, Top
2857 @c    node-name, next, previous, up
2858 @chapter Extending GiNaC
2859
2860 By reading so far you should have gotten a fairly good understanding of
2861 GiNaC's design-patterns.  From here on you should start reading the
2862 sources.  All we can do now is issue some recommendations how to tackle
2863 GiNaC's many loose ends in order to fulfill everybody's dreams.  If you
2864 develop some useful extension please don't hesitate to contact the GiNaC
2865 authors---they will happily incorporate them into future versions.
2866
2867 @menu
2868 * What does not belong into GiNaC::  What to avoid.
2869 * Symbolic functions::               Implementing symbolic functions.
2870 * Adding classes::                   Defining new algebraic classes.
2871 @end menu
2872
2873
2874 @node What does not belong into GiNaC, Symbolic functions, Extending GiNaC, Extending GiNaC
2875 @c    node-name, next, previous, up
2876 @section What doesn't belong into GiNaC
2877
2878 @cindex @command{ginsh}
2879 First of all, GiNaC's name must be read literally.  It is designed to be
2880 a library for use within C++.  The tiny @command{ginsh} accompanying
2881 GiNaC makes this even more clear: it doesn't even attempt to provide a
2882 language.  There are no loops or conditional expressions in
2883 @command{ginsh}, it is merely a window into the library for the
2884 programmer to test stuff (or to show off).  Still, the design of a
2885 complete CAS with a language of its own, graphical capabilites and all
2886 this on top of GiNaC is possible and is without doubt a nice project for
2887 the future.
2888
2889 There are many built-in functions in GiNaC that do not know how to
2890 evaluate themselves numerically to a precision declared at runtime
2891 (using @code{Digits}).  Some may be evaluated at certain points, but not
2892 generally.  This ought to be fixed.  However, doing numerical
2893 computations with GiNaC's quite abstract classes is doomed to be
2894 inefficient.  For this purpose, the underlying foundation classes
2895 provided by @acronym{CLN} are much better suited.
2896
2897
2898 @node Symbolic functions, Adding classes, What does not belong into GiNaC, Extending GiNaC
2899 @c    node-name, next, previous, up
2900 @section Symbolic functions
2901
2902 The easiest and most instructive way to start with is probably to
2903 implement your own function.  GiNaC's functions are objects of class
2904 @code{function}.  The preprocessor is then used to convert the function
2905 names to objects with a corresponding serial number that is used
2906 internally to identify them.  You usually need not worry about this
2907 number.  New functions may be inserted into the system via a kind of
2908 `registry'.  It is your responsibility to care for some functions that
2909 are called when the user invokes certain methods.  These are usual
2910 C++-functions accepting a number of @code{ex} as arguments and returning
2911 one @code{ex}.  As an example, if we have a look at a simplified
2912 implementation of the cosine trigonometric function, we first need a
2913 function that is called when one wishes to @code{eval} it.  It could
2914 look something like this:
2915
2916 @example
2917 static ex cos_eval_method(const ex & x)
2918 @{
2919     // if (!x%(2*Pi)) return 1
2920     // if (!x%Pi) return -1
2921     // if (!x%Pi/2) return 0
2922     // care for other cases...
2923     return cos(x).hold();
2924 @}
2925 @end example
2926
2927 @cindex @code{hold()}
2928 @cindex evaluation
2929 The last line returns @code{cos(x)} if we don't know what else to do and
2930 stops a potential recursive evaluation by saying @code{.hold()}, which
2931 sets a flag to the expression signaling that it has been evaluated.  We
2932 should also implement a method for numerical evaluation and since we are
2933 lazy we sweep the problem under the rug by calling someone else's
2934 function that does so, in this case the one in class @code{numeric}:
2935
2936 @example
2937 static ex cos_evalf(const ex & x)
2938 @{
2939     return cos(ex_to_numeric(x));
2940 @}
2941 @end example
2942
2943 Differentiation will surely turn up and so we need to tell @code{cos}
2944 what the first derivative is (higher derivatives (@code{.diff(x,3)} for
2945 instance are then handled automatically by @code{basic::diff} and
2946 @code{ex::diff}):
2947
2948 @example
2949 static ex cos_deriv(const ex & x, unsigned diff_param)
2950 @{
2951     return -sin(x);
2952 @}
2953 @end example
2954
2955 @cindex product rule
2956 The second parameter is obligatory but uninteresting at this point.  It
2957 specifies which parameter to differentiate in a partial derivative in
2958 case the function has more than one parameter and its main application
2959 is for correct handling of the chain rule.  For Taylor expansion, it is
2960 enough to know how to differentiate.  But if the function you want to
2961 implement does have a pole somewhere in the complex plane, you need to
2962 write another method for Laurent expansion around that point.
2963
2964 Now that all the ingredients for @code{cos} have been set up, we need
2965 to tell the system about it.  This is done by a macro and we are not
2966 going to descibe how it expands, please consult your preprocessor if you
2967 are curious:
2968
2969 @example
2970 REGISTER_FUNCTION(cos, eval_func(cos_eval).
2971                        evalf_func(cos_evalf).
2972                        derivative_func(cos_deriv));
2973 @end example
2974
2975 The first argument is the function's name used for calling it and for
2976 output.  The second binds the corresponding methods as options to this
2977 object.  Options are separated by a dot and can be given in an arbitrary
2978 order.  GiNaC functions understand several more options which are always
2979 specified as @code{.option(params)}, for example a method for series
2980 expansion @code{.series_func(cos_series)}.  Again, if no series
2981 expansion method is given, GiNaC defaults to simple Taylor expansion,
2982 which is correct if there are no poles involved as is the case for the
2983 @code{cos} function.  The way GiNaC handles poles in case there are any
2984 is best understood by studying one of the examples, like the Gamma
2985 (@code{tgamma}) function for instance.  (In essence the function first
2986 checks if there is a pole at the evaluation point and falls back to
2987 Taylor expansion if there isn't.  Then, the pole is regularized by some
2988 suitable transformation.)  Also, the new function needs to be declared
2989 somewhere.  This may also be done by a convenient preprocessor macro:
2990
2991 @example
2992 DECLARE_FUNCTION_1P(cos)
2993 @end example
2994
2995 The suffix @code{_1P} stands for @emph{one parameter}.  Of course, this
2996 implementation of @code{cos} is very incomplete and lacks several safety
2997 mechanisms.  Please, have a look at the real implementation in GiNaC.
2998 (By the way: in case you are worrying about all the macros above we can
2999 assure you that functions are GiNaC's most macro-intense classes.  We
3000 have done our best to avoid macros where we can.)
3001
3002
3003 @node Adding classes, A Comparison With Other CAS, Symbolic functions, Extending GiNaC
3004 @c    node-name, next, previous, up
3005 @section Adding classes
3006
3007 If you are doing some very specialized things with GiNaC you may find that
3008 you have to implement your own algebraic classes to fit your needs. This
3009 section will explain how to do this by giving the example of a simple
3010 'string' class. After reading this section you will know how to properly
3011 declare a GiNaC class and what the minimum required member functions are
3012 that you have to implement. We only cover the implementation of a 'leaf'
3013 class here (i.e. one that doesn't contain subexpressions). Creating a
3014 container class like, for example, a class representing tensor products is
3015 more involved but this section should give you enough information so you can
3016 consult the source to GiNaC's predefined classes if you want to implement
3017 something more complicated.
3018
3019 @subsection GiNaC's run-time type information system
3020
3021 @cindex hierarchy of classes
3022 @cindex RTTI
3023 All algebraic classes (that is, all classes that can appear in expressions)
3024 in GiNaC are direct or indirect subclasses of the class @code{basic}. So a
3025 @code{basic *} (which is essentially what an @code{ex} is) represents a
3026 generic pointer to an algebraic class. Occasionally it is necessary to find
3027 out what the class of an object pointed to by a @code{basic *} really is.
3028 Also, for the unarchiving of expressions it must be possible to find the
3029 @code{unarchive()} function of a class given the class name (as a string). A
3030 system that provides this kind of information is called a run-time type
3031 information (RTTI) system. The C++ language provides such a thing (see the
3032 standard header file @file{<typeinfo>}) but for efficiency reasons GiNaC
3033 implements its own, simpler RTTI.
3034
3035 The RTTI in GiNaC is based on two mechanisms:
3036
3037 @itemize @bullet
3038
3039 @item
3040 The @code{basic} class declares a member variable @code{tinfo_key} which
3041 holds an unsigned integer that identifies the object's class. These numbers
3042 are defined in the @file{tinfos.h} header file for the built-in GiNaC
3043 classes. They all start with @code{TINFO_}.
3044
3045 @item
3046 By means of some clever tricks with static members, GiNaC maintains a list
3047 of information for all classes derived from @code{basic}. The information
3048 available includes the class names, the @code{tinfo_key}s, and pointers
3049 to the unarchiving functions. This class registry is defined in the
3050 @file{registrar.h} header file.
3051
3052 @end itemize
3053
3054 The disadvantage of this proprietary RTTI implementation is that there's
3055 a little more to do when implementing new classes (C++'s RTTI works more
3056 or less automatic) but don't worry, most of the work is simplified by
3057 macros.
3058
3059 @subsection A minimalistic example
3060
3061 Now we will start implementing a new class @code{mystring} that allows
3062 placing character strings in algebraic expressions (this is not very useful,
3063 but it's just an example). This class will be a direct subclass of
3064 @code{basic}. You can use this sample implementation as a starting point
3065 for your own classes.
3066
3067 The code snippets given here assume that you have included some header files
3068 as follows:
3069
3070 @example
3071 #include <iostream>
3072 #include <string>   
3073 #include <stdexcept>
3074 using namespace std;
3075
3076 #include <ginac/ginac.h>
3077 using namespace GiNaC;
3078 @end example
3079
3080 The first thing we have to do is to define a @code{tinfo_key} for our new
3081 class. This can be any arbitrary unsigned number that is not already taken
3082 by one of the existing classes but it's better to come up with something
3083 that is unlikely to clash with keys that might be added in the future. The
3084 numbers in @file{tinfos.h} are modeled somewhat after the class hierarchy
3085 which is not a requirement but we are going to stick with this scheme:
3086
3087 @example
3088 const unsigned TINFO_mystring = 0x42420001U;
3089 @end example
3090
3091 Now we can write down the class declaration. The class stores a C++
3092 @code{string} and the user shall be able to construct a @code{mystring}
3093 object from a C or C++ string:
3094
3095 @example
3096 class mystring : public basic
3097 @{
3098     GINAC_DECLARE_REGISTERED_CLASS(mystring, basic)
3099   
3100 public:
3101     mystring(const string &s);
3102     mystring(const char *s);
3103
3104 private:
3105     string str;
3106 @};
3107
3108 GIANC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
3109 @end example
3110
3111 The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
3112 macros are defined in @file{registrar.h}. They take the name of the class
3113 and its direct superclass as arguments and insert all required declarations
3114 for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
3115 the first line after the opening brace of the class definition. The
3116 @code{GINAC_IMPLEMENT_REGISTERED_CLASS} may appear anywhere else in the
3117 source (at global scope, of course, not inside a function).
3118
3119 @code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
3120 declarations of the default and copy constructor, the destructor, the
3121 assignment operator and a couple of other functions that are required. It
3122 also defines a type @code{inherited} which refers to the superclass so you
3123 don't have to modify your code every time you shuffle around the class
3124 hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} implements the copy
3125 constructor, the destructor and the assignment operator.
3126
3127 Now there are nine member functions we have to implement to get a working
3128 class:
3129
3130 @itemize
3131
3132 @item
3133 @code{mystring()}, the default constructor.
3134
3135 @item
3136 @code{void destroy(bool call_parent)}, which is used in the destructor and the
3137 assignment operator to free dynamically allocated members. The @code{call_parent}
3138 specifies whether the @code{destroy()} function of the superclass is to be
3139 called also.
3140
3141 @item
3142 @code{void copy(const mystring &other)}, which is used in the copy constructor
3143 and assignment operator to copy the member variables over from another
3144 object of the same class.
3145
3146 @item
3147 @code{void archive(archive_node &n)}, the archiving function. This stores all
3148 information needed to reconstruct an object of this class inside an
3149 @code{archive_node}.
3150
3151 @item
3152 @code{mystring(const archive_node &n, const lst &sym_lst)}, the unarchiving
3153 constructor. This constructs an instance of the class from the information
3154 found in an @code{archive_node}.
3155
3156 @item
3157 @code{ex unarchive(const archive_node &n, const lst &sym_lst)}, the static
3158 unarchiving function. It constructs a new instance by calling the unarchiving
3159 constructor.
3160
3161 @item
3162 @code{int compare_same_type(const basic &other)}, which is used internally
3163 by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
3164 -1, depending on the relative order of this object and the @code{other}
3165 object. If it returns 0, the objects are considered equal.
3166 @strong{Note:} This has nothing to do with the (numeric) ordering
3167 relationship expressed by @code{<}, @code{>=} etc (which cannot be defined
3168 for non-numeric classes). For example, @code{numeric(1).compare_same_type(numeric(2))}
3169 may return +1 even though 1 is clearly smaller than 2. Every GiNaC class
3170 must provide a @code{compare_same_type()} function, even those representing
3171 objects for which no reasonable algebraic ordering relationship can be
3172 defined.
3173
3174 @item
3175 And, of course, @code{mystring(const string &s)} and @code{mystring(const char *s)}
3176 which are the two constructors we declared.
3177
3178 @end itemize
3179
3180 Let's proceed step-by-step. The default constructor looks like this:
3181
3182 @example
3183 mystring::mystring() : inherited(TINFO_mystring)
3184 @{
3185     // dynamically allocate resources here if required
3186 @}
3187 @end example
3188
3189 The golden rule is that in all constructors you have to set the
3190 @code{tinfo_key} member to the @code{TINFO_*} value of your class. Otherwise
3191 it will be set by the constructor of the superclass and all hell will break
3192 loose in the RTTI. For your convenience, the @code{basic} class provides
3193 a constructor that takes a @code{tinfo_key} value, which we are using here
3194 (remember that in our case @code{inherited = basic}). If the superclass
3195 didn't have such a constructor, we would have to set the @code{tinfo_key}
3196 to the right value manually.
3197
3198 In the default constructor you should set all other member variables to
3199 reasonable default values (we don't need that here since our @code{str}
3200 member gets set to an empty string automatically). The constructor(s) are of
3201 course also the right place to allocate any dynamic resources you require.
3202
3203 Next, the @code{destroy()} function:
3204
3205 @example
3206 void mystring::destroy(bool call_parent)
3207 @{
3208     // free dynamically allocated resources here if required
3209     if (call_parent)
3210         inherited::destroy(call_parent);
3211 @}
3212 @end example
3213
3214 This function is where we free all dynamically allocated resources. We don't
3215 have any so we're not doing anything here, but if we had, for example, used
3216 a C-style @code{char *} to store our string, this would be the place to
3217 @code{delete[]} the string storage. If @code{call_parent} is true, we have
3218 to call the @code{destroy()} function of the superclass after we're done
3219 (to mimic C++'s automatic invocation of superclass destructors where
3220 @code{destroy()} is called from outside a destructor).
3221
3222 The @code{copy()} function just copies over the member variables from
3223 another object:
3224
3225 @example
3226 void mystring::copy(const mystring &other)
3227 @{
3228     inherited::copy(other);
3229     str = other.str;
3230 @}
3231 @end example
3232
3233 We can simply overwrite the member variables here. There's no need to worry
3234 about dynamically allocated storage. The assignment operator (which is
3235 automatically defined by @code{GINAC_IMPLEMENT_REGISTERED_CLASS}, as you
3236 recall) calls @code{destroy()} before it calls @code{copy()}. You have to
3237 explicitly call the @code{copy()} function of the superclass here so
3238 all the member variables will get copied.
3239
3240 Next are the three functions for archiving. You have to implement them even
3241 if you don't plan to use archives, but the minimum required implementation
3242 is really simple. First, the archiving function:
3243
3244 @example
3245 void mystring::archive(archive_node &n) const
3246 @{
3247     inherited::archive(n);
3248     n.add_string("string", str);
3249 @}
3250 @end example
3251
3252 The only thing that is really required is calling the @code{archive()}
3253 function of the superclass. Optionally, you can store all information you
3254 deem necessary for representing the object into the passed
3255 @code{archive_node}. We are just storing our string here. For more
3256 information on how the archiving works, consult the @file{archive.h} header
3257 file.
3258
3259 The unarchiving constructor is basically the inverse of the archiving
3260 function:
3261
3262 @example
3263 mystring::mystring(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
3264 @{
3265     n.find_string("string", str);
3266 @}
3267 @end example
3268
3269 If you don't need archiving, just leave this function empty (but you must
3270 invoke the unarchiving constructor of the superclass). Note that we don't
3271 have to set the @code{tinfo_key} here because it is done automatically
3272 by the unarchiving constructor of the @code{basic} class.
3273
3274 Finally, the unarchiving function:
3275
3276 @example
3277 ex mystring::unarchive(const archive_node &n, const lst &sym_lst)
3278 @{
3279     return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
3280 @}
3281 @end example
3282
3283 You don't have to understand how exactly this works. Just copy these four
3284 lines into your code literally (replacing the class name, of course). It
3285 calls the unarchiving constructor of the class and unless you are doing
3286 something very special (like matching @code{archive_node}s to global
3287 objects) you don't need a different implementation. For those who are
3288 interested: setting the @code{dynallocated} flag puts the object under
3289 the control of GiNaC's garbage collection. It will get deleted automatically
3290 once it is no longer referenced.
3291
3292 Our @code{compare_same_type()} function uses a provided function to compare
3293 the string members:
3294
3295 @example
3296 int mystring::compare_same_type(const basic &other) const
3297 @{
3298     const mystring &o = static_cast<const mystring &>(other);
3299     int cmpval = str.compare(o.str);
3300     if (cmpval == 0)
3301         return 0;
3302     else if (cmpval < 0)
3303         return -1;
3304     else
3305         return 1;
3306 @}
3307 @end example
3308
3309 Although this function takes a @code{basic &}, it will always be a reference
3310 to an object of exactly the same class (objects of different classes are not
3311 comparable), so the cast is safe. If this function returns 0, the two objects
3312 are considered equal (in the sense that @math{A-B=0}), so you should compare
3313 all relevant member variables.
3314
3315 Now the only thing missing is our two new constructors:
3316
3317 @example
3318 mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s)
3319 @{
3320     // dynamically allocate resources here if required
3321 @}
3322
3323 mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s)
3324 @{
3325     // dynamically allocate resources here if required
3326 @}
3327 @end example
3328
3329 No surprises here. We set the @code{str} member from the argument and
3330 remember to pass the right @code{tinfo_key} to the @code{basic} constructor.
3331
3332 That's it! We now have a minimal working GiNaC class that can store
3333 strings in algebraic expressions. Let's confirm that the RTTI works:
3334
3335 @example
3336 ex e = mystring("Hello, world!");
3337 cout << is_ex_of_type(e, mystring) << endl;
3338  // -> 1 (true)
3339
3340 cout << e.bp->class_name() << endl;
3341  // -> mystring
3342 @end example
3343
3344 Obviously it does. Let's see what the expression @code{e} looks like:
3345
3346 @example
3347 cout << e << endl;
3348  // -> [mystring object]
3349 @end example
3350
3351 Hm, not exactly what we expect, but of course the @code{mystring} class
3352 doesn't yet know how to print itself. This is done in the @code{print()}
3353 member function. Let's say that we wanted to print the string surrounded
3354 by double quotes:
3355
3356 @example
3357 class mystring : public basic
3358 @{
3359     ...
3360 public:
3361     void print(ostream &os, unsigned upper_precedence) const;
3362     ...
3363 @};
3364
3365 void mystring::print(ostream &os, unsigned upper_precedence) const
3366 @{
3367     os << '\"' << str << '\"';
3368 @}
3369 @end example
3370
3371 The @code{upper_precedence} argument is only required for container classes
3372 to correctly parenthesize the output. Let's try again to print the expression:
3373
3374 @example
3375 cout << e << endl;
3376  // -> "Hello, world!"
3377 @end example
3378
3379 Much better. The @code{mystring} class can be used in arbitrary expressions:
3380
3381 @example
3382 e += mystring("GiNaC rulez"); 
3383 cout << e << endl;
3384  // -> "GiNaC rulez"+"Hello, world!"
3385 @end example
3386
3387 (note that GiNaC's automatic term reordering is in effect here), or even
3388
3389 @example
3390 e = pow(mystring("One string"), 2*sin(Pi-mystring("Another string")));
3391 cout << e << endl;
3392  // -> "One string"^(2*sin(-"Another string"+Pi))
3393 @end example
3394
3395 Whether this makes sense is debatable but remember that this is only an
3396 example. At least it allows you to implement your own symbolic algorithms
3397 for your objects.
3398
3399 Note that GiNaC's algebraic rules remain unchanged:
3400
3401 @example
3402 e = mystring("Wow") * mystring("Wow");
3403 cout << e << endl;
3404  // -> "Wow"^2
3405
3406 e = pow(mystring("First")-mystring("Second"), 2);
3407 cout << e.expand() << endl;
3408  // -> -2*"First"*"Second"+"First"^2+"Second"^2
3409 @end example
3410
3411 There's no way to, for example, make GiNaC's @code{add} class perform string
3412 concatenation. You would have to implement this yourself.
3413
3414 @subsection Automatic evaluation
3415
3416 @cindex @code{hold()}
3417 @cindex evaluation
3418 When dealing with objects that are just a little more complicated than the
3419 simple string objects we have implemented, chances are that you will want to
3420 have some automatic simplifications or canonicalizations performed on them.
3421 This is done in the evaluation member function @code{eval()}. Let's say that
3422 we wanted all strings automatically converted to lowercase with
3423 non-alphabetic characters stripped, and empty strings removed:
3424
3425 @example
3426 class mystring : public basic
3427 @{
3428     ...
3429 public:
3430     ex eval(int level = 0) const;
3431     ...
3432 @};
3433
3434 ex mystring::eval(int level) const
3435 @{
3436     string new_str;
3437     for (int i=0; i<str.length(); i++) @{
3438         char c = str[i];
3439         if (c >= 'A' && c <= 'Z') 
3440             new_str += tolower(c);
3441         else if (c >= 'a' && c <= 'z')
3442             new_str += c;
3443     @}
3444
3445     if (new_str.length() == 0)
3446         return _ex0();
3447     else
3448         return mystring(new_str).hold();
3449 @}
3450 @end example
3451
3452 The @code{level} argument is used to limit the recursion depth of the
3453 evaluation. We don't have any subexpressions in the @code{mystring} class
3454 so we are not concerned with this. If we had, we would call the @code{eval()}
3455 functions of the subexpressions with @code{level - 1} as the argument if
3456 @code{level != 1}. The @code{hold()} member function sets a flag in the
3457 object that prevents further evaluation. Otherwise we might end up in an
3458 endless loop. When you want to return the object unmodified, use
3459 @code{return this->hold();}.
3460
3461 Let's confirm that it works:
3462
3463 @example
3464 ex e = mystring("Hello, world!") + mystring("!?#");
3465 cout << e << endl;
3466  // -> "helloworld"
3467
3468 e = mystring("Wow!") + mystring("WOW") + mystring(" W ** o ** W");  
3469 cout << e << endl;
3470  // -> 3*"wow"
3471 @end example
3472
3473 @subsection Other member functions
3474
3475 We have implemented only a small set of member functions to make the class
3476 work in the GiNaC framework. For a real algebraic class, there are probably
3477 some more functions that you will want to re-implement, such as
3478 @code{evalf()}, @code{series()} or @code{op()}. Have a look at @file{basic.h}
3479 or the header file of the class you want to make a subclass of to see
3480 what's there. You can, of course, also add your own new member functions.
3481 In this case you will probably want to define a little helper function like
3482
3483 @example
3484 inline const mystring &ex_to_mystring(const ex &e)
3485 @{
3486     return static_cast<const mystring &>(*e.bp);
3487 @}
3488 @end example
3489
3490 that let's you get at the object inside an expression (after you have verified
3491 that the type is correct) so you can call member functions that are specific
3492 to the class.
3493
3494 That's it. May the source be with you!
3495
3496
3497 @node A Comparison With Other CAS, Advantages, Adding classes, Top
3498 @c    node-name, next, previous, up
3499 @chapter A Comparison With Other CAS
3500 @cindex advocacy
3501
3502 This chapter will give you some information on how GiNaC compares to
3503 other, traditional Computer Algebra Systems, like @emph{Maple},
3504 @emph{Mathematica} or @emph{Reduce}, where it has advantages and
3505 disadvantages over these systems.
3506
3507 @menu
3508 * Advantages::                       Stengths of the GiNaC approach.
3509 * Disadvantages::                    Weaknesses of the GiNaC approach.
3510 * Why C++?::                         Attractiveness of C++.
3511 @end menu
3512
3513 @node Advantages, Disadvantages, A Comparison With Other CAS, A Comparison With Other CAS
3514 @c    node-name, next, previous, up
3515 @section Advantages
3516
3517 GiNaC has several advantages over traditional Computer
3518 Algebra Systems, like 
3519
3520 @itemize @bullet
3521
3522 @item
3523 familiar language: all common CAS implement their own proprietary
3524 grammar which you have to learn first (and maybe learn again when your
3525 vendor decides to `enhance' it).  With GiNaC you can write your program
3526 in common C++, which is standardized.
3527
3528 @cindex STL
3529 @item
3530 structured data types: you can build up structured data types using
3531 @code{struct}s or @code{class}es together with STL features instead of
3532 using unnamed lists of lists of lists.
3533
3534 @item
3535 strongly typed: in CAS, you usually have only one kind of variables
3536 which can hold contents of an arbitrary type.  This 4GL like feature is
3537 nice for novice programmers, but dangerous.
3538     
3539 @item
3540 development tools: powerful development tools exist for C++, like fancy
3541 editors (e.g. with automatic indentation and syntax highlighting),
3542 debuggers, visualization tools, documentation generators...
3543
3544 @item
3545 modularization: C++ programs can easily be split into modules by
3546 separating interface and implementation.
3547
3548 @item
3549 price: GiNaC is distributed under the GNU Public License which means
3550 that it is free and available with source code.  And there are excellent
3551 C++-compilers for free, too.
3552     
3553 @item
3554 extendable: you can add your own classes to GiNaC, thus extending it on
3555 a very low level.  Compare this to a traditional CAS that you can
3556 usually only extend on a high level by writing in the language defined
3557 by the parser.  In particular, it turns out to be almost impossible to
3558 fix bugs in a traditional system.
3559
3560 @item
3561 multiple interfaces: Though real GiNaC programs have to be written in
3562 some editor, then be compiled, linked and executed, there are more ways
3563 to work with the GiNaC engine.  Many people want to play with
3564 expressions interactively, as in traditional CASs.  Currently, two such
3565 windows into GiNaC have been implemented and many more are possible: the
3566 tiny @command{ginsh} that is part of the distribution exposes GiNaC's
3567 types to a command line and second, as a more consistent approach, an
3568 interactive interface to the @acronym{Cint} C++ interpreter has been put
3569 together (called @acronym{GiNaC-cint}) that allows an interactive
3570 scripting interface consistent with the C++ language.
3571
3572 @item
3573 seemless integration: it is somewhere between difficult and impossible
3574 to call CAS functions from within a program written in C++ or any other
3575 programming language and vice versa.  With GiNaC, your symbolic routines
3576 are part of your program.  You can easily call third party libraries,
3577 e.g. for numerical evaluation or graphical interaction.  All other
3578 approaches are much more cumbersome: they range from simply ignoring the
3579 problem (i.e. @emph{Maple}) to providing a method for `embedding' the
3580 system (i.e. @emph{Yacas}).
3581
3582 @item
3583 efficiency: often large parts of a program do not need symbolic
3584 calculations at all.  Why use large integers for loop variables or
3585 arbitrary precision arithmetics where @code{int} and @code{double} are
3586 sufficient?  For pure symbolic applications, GiNaC is comparable in
3587 speed with other CAS.
3588
3589 @end itemize
3590
3591
3592 @node Disadvantages, Why C++?, Advantages, A Comparison With Other CAS
3593 @c    node-name, next, previous, up
3594 @section Disadvantages
3595
3596 Of course it also has some disadvantages:
3597
3598 @itemize @bullet
3599
3600 @item
3601 advanced features: GiNaC cannot compete with a program like
3602 @emph{Reduce} which exists for more than 30 years now or @emph{Maple}
3603 which grows since 1981 by the work of dozens of programmers, with
3604 respect to mathematical features.  Integration, factorization,
3605 non-trivial simplifications, limits etc. are missing in GiNaC (and are
3606 not planned for the near future).
3607
3608 @item
3609 portability: While the GiNaC library itself is designed to avoid any
3610 platform dependent features (it should compile on any ANSI compliant C++
3611 compiler), the currently used version of the CLN library (fast large
3612 integer and arbitrary precision arithmetics) can be compiled only on
3613 systems with a recently new C++ compiler from the GNU Compiler
3614 Collection (@acronym{GCC}).@footnote{This is because CLN uses
3615 PROVIDE/REQUIRE like macros to let the compiler gather all static
3616 initializations, which works for GNU C++ only.}  GiNaC uses recent
3617 language features like explicit constructors, mutable members, RTTI,
3618 @code{dynamic_cast}s and STL, so ANSI compliance is meant literally.
3619 Recent @acronym{GCC} versions starting at 2.95, although itself not yet
3620 ANSI compliant, support all needed features.
3621     
3622 @end itemize
3623
3624
3625 @node Why C++?, Internal Structures, Disadvantages, A Comparison With Other CAS
3626 @c    node-name, next, previous, up
3627 @section Why C++?
3628
3629 Why did we choose to implement GiNaC in C++ instead of Java or any other
3630 language?  C++ is not perfect: type checking is not strict (casting is
3631 possible), separation between interface and implementation is not
3632 complete, object oriented design is not enforced.  The main reason is
3633 the often scolded feature of operator overloading in C++.  While it may
3634 be true that operating on classes with a @code{+} operator is rarely
3635 meaningful, it is perfectly suited for algebraic expressions.  Writing
3636 @math{3x+5y} as @code{3*x+5*y} instead of
3637 @code{x.times(3).plus(y.times(5))} looks much more natural.
3638 Furthermore, the main developers are more familiar with C++ than with
3639 any other programming language.
3640
3641
3642 @node Internal Structures, Expressions are reference counted, Why C++? , Top
3643 @c    node-name, next, previous, up
3644 @appendix Internal Structures
3645
3646 @menu
3647 * Expressions are reference counted::
3648 * Internal representation of products and sums::
3649 @end menu
3650
3651 @node Expressions are reference counted, Internal representation of products and sums, Internal Structures, Internal Structures
3652 @c    node-name, next, previous, up
3653 @appendixsection Expressions are reference counted
3654
3655 @cindex reference counting
3656 @cindex copy-on-write
3657 @cindex garbage collection
3658 An expression is extremely light-weight since internally it works like a
3659 handle to the actual representation and really holds nothing more than a
3660 pointer to some other object. What this means in practice is that
3661 whenever you create two @code{ex} and set the second equal to the first
3662 no copying process is involved. Instead, the copying takes place as soon
3663 as you try to change the second.  Consider the simple sequence of code:
3664
3665 @example
3666 #include <ginac/ginac.h>
3667 using namespace GiNaC;
3668
3669 int main()
3670 @{
3671     symbol x("x"), y("y"), z("z");
3672     ex e1, e2;
3673
3674     e1 = sin(x + 2*y) + 3*z + 41;
3675     e2 = e1;                // e2 points to same object as e1
3676     cout << e2 << endl;     // prints sin(x+2*y)+3*z+41
3677     e2 += 1;                // e2 is copied into a new object
3678     cout << e2 << endl;     // prints sin(x+2*y)+3*z+42
3679 @}
3680 @end example
3681
3682 The line @code{e2 = e1;} creates a second expression pointing to the
3683 object held already by @code{e1}.  The time involved for this operation
3684 is therefore constant, no matter how large @code{e1} was.  Actual
3685 copying, however, must take place in the line @code{e2 += 1;} because
3686 @code{e1} and @code{e2} are not handles for the same object any more.
3687 This concept is called @dfn{copy-on-write semantics}.  It increases
3688 performance considerably whenever one object occurs multiple times and
3689 represents a simple garbage collection scheme because when an @code{ex}
3690 runs out of scope its destructor checks whether other expressions handle
3691 the object it points to too and deletes the object from memory if that
3692 turns out not to be the case.  A slightly less trivial example of
3693 differentiation using the chain-rule should make clear how powerful this
3694 can be:
3695
3696 @example
3697 #include <ginac/ginac.h>
3698 using namespace GiNaC;
3699
3700 int main()
3701 @{
3702     symbol x("x"), y("y");
3703
3704     ex e1 = x + 3*y;
3705     ex e2 = pow(e1, 3);
3706     ex e3 = diff(sin(e2), x);   // first derivative of sin(e2) by x
3707     cout << e1 << endl          // prints x+3*y
3708          << e2 << endl          // prints (x+3*y)^3
3709          << e3 << endl;         // prints 3*(x+3*y)^2*cos((x+3*y)^3)
3710 @}
3711 @end example
3712
3713 Here, @code{e1} will actually be referenced three times while @code{e2}
3714 will be referenced two times.  When the power of an expression is built,
3715 that expression needs not be copied.  Likewise, since the derivative of
3716 a power of an expression can be easily expressed in terms of that
3717 expression, no copying of @code{e1} is involved when @code{e3} is
3718 constructed.  So, when @code{e3} is constructed it will print as
3719 @code{3*(x+3*y)^2*cos((x+3*y)^3)} but the argument of @code{cos()} only
3720 holds a reference to @code{e2} and the factor in front is just
3721 @code{3*e1^2}.
3722
3723 As a user of GiNaC, you cannot see this mechanism of copy-on-write