simplify_indexed() recognizes linear combinations of matrices with numeric
[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 symmetric metric
1644 tensor with two indices that can be used to raise/lower tensor indices. The
1645 metric 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 (linear combinations and products of vectors and matrices, traces
1728 and scalar products):
1729
1730 @example
1731 @{
1732     idx i(symbol("i"), 2), j(symbol("j"), 2);
1733     symbol x("x"), y("y");
1734
1735     matrix A(2, 2, lst(1, 2, 3, 4)), X(2, 1, lst(x, y));
1736
1737     cout << indexed(A, i, i) << endl;
1738      // -> 5
1739
1740     ex e = indexed(A, i, j) * indexed(X, j);
1741     cout << e.simplify_indexed() << endl;
1742      // -> [[ [[2*y+x]], [[4*y+3*x]] ]].i
1743
1744     e = indexed(A, i, j) * indexed(X, i) + indexed(X, j) * 2;
1745     cout << e.simplify_indexed() << endl;
1746      // -> [[ [[3*y+3*x,6*y+2*x]] ]].j
1747 @}
1748 @end example
1749
1750 You can of course obtain the same results with the @code{matrix::add()},
1751 @code{matrix::mul()} and @code{matrix::trace()} methods but with indices you
1752 don't have to worry about transposing matrices.
1753
1754 Matrix indices always start at 0 and their dimension must match the number
1755 of rows/columns of the matrix. Matrices with one row or one column are
1756 vectors and can have one or two indices (it doesn't matter whether it's a
1757 row or a column vector). Other matrices must have two indices.
1758
1759 You should be careful when using indices with variance on matrices. GiNaC
1760 doesn't look at the variance and doesn't know that @samp{F~mu~nu} and
1761 @samp{F.mu.nu} are different matrices. In this case you should use only
1762 one form for @samp{F} and explicitly multiply it with a matrix representation
1763 of the metric tensor.
1764
1765
1766 @node Methods and Functions, Information About Expressions, Indexed objects, Top
1767 @c    node-name, next, previous, up
1768 @chapter Methods and Functions
1769 @cindex polynomial
1770
1771 In this chapter the most important algorithms provided by GiNaC will be
1772 described.  Some of them are implemented as functions on expressions,
1773 others are implemented as methods provided by expression objects.  If
1774 they are methods, there exists a wrapper function around it, so you can
1775 alternatively call it in a functional way as shown in the simple
1776 example:
1777
1778 @example
1779 #include <ginac/ginac.h>
1780 using namespace GiNaC;
1781
1782 int main()
1783 @{
1784     ex x = numeric(1.0);
1785     
1786     cout << "As method:   " << sin(x).evalf() << endl;
1787     cout << "As function: " << evalf(sin(x)) << endl;
1788 @}
1789 @end example
1790
1791 @cindex @code{subs()}
1792 The general rule is that wherever methods accept one or more parameters
1793 (@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
1794 wrapper accepts is the same but preceded by the object to act on
1795 (@var{object}, @var{arg1}, @var{arg2}, @dots{}).  This approach is the
1796 most natural one in an OO model but it may lead to confusion for MapleV
1797 users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
1798 would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
1799 @code{A} and @code{x}).  On the other hand, since MapleV returns 3 on
1800 @code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
1801 coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
1802 here.  Also, users of MuPAD will in most cases feel more comfortable
1803 with GiNaC's convention.  All function wrappers are implemented
1804 as simple inline functions which just call the corresponding method and
1805 are only provided for users uncomfortable with OO who are dead set to
1806 avoid method invocations.  Generally, nested function wrappers are much
1807 harder to read than a sequence of methods and should therefore be
1808 avoided if possible.  On the other hand, not everything in GiNaC is a
1809 method on class @code{ex} and sometimes calling a function cannot be
1810 avoided.
1811
1812 @menu
1813 * Information About Expressions::
1814 * Substituting Symbols::
1815 * Polynomial Arithmetic::           Working with polynomials.
1816 * Rational Expressions::            Working with rational functions.
1817 * Symbolic Differentiation::
1818 * Series Expansion::                Taylor and Laurent expansion.
1819 * Built-in Functions::              List of predefined mathematical functions.
1820 * Input/Output::                    Input and output of expressions.
1821 @end menu
1822
1823
1824 @node Information About Expressions, Substituting Symbols, Methods and Functions, Methods and Functions
1825 @c    node-name, next, previous, up
1826 @section Getting information about expressions
1827
1828 @subsection Checking expression types
1829 @cindex @code{is_ex_of_type()}
1830 @cindex @code{ex_to_numeric()}
1831 @cindex @code{ex_to_@dots{}}
1832 @cindex @code{Converting ex to other classes}
1833 @cindex @code{info()}
1834
1835 Sometimes it's useful to check whether a given expression is a plain number,
1836 a sum, a polynomial with integer coefficients, or of some other specific type.
1837 GiNaC provides two functions for this (the first one is actually a macro):
1838
1839 @example
1840 bool is_ex_of_type(const ex & e, TYPENAME t);
1841 bool ex::info(unsigned flag);
1842 @end example
1843
1844 When the test made by @code{is_ex_of_type()} returns true, it is safe to
1845 call one of the functions @code{ex_to_@dots{}}, where @code{@dots{}} is
1846 one of the class names (@xref{The Class Hierarchy}, for a list of all
1847 classes). For example, assuming @code{e} is an @code{ex}:
1848
1849 @example
1850 @{
1851     @dots{}
1852     if (is_ex_of_type(e, numeric))
1853         numeric n = ex_to_numeric(e);
1854     @dots{}
1855 @}
1856 @end example
1857
1858 @code{is_ex_of_type()} allows you to check whether the top-level object of
1859 an expression @samp{e} is an instance of the GiNaC class @samp{t}
1860 (@xref{The Class Hierarchy}, for a list of all classes). This is most useful,
1861 e.g., for checking whether an expression is a number, a sum, or a product:
1862
1863 @example
1864 @{
1865     symbol x("x");
1866     ex e1 = 42;
1867     ex e2 = 4*x - 3;
1868     is_ex_of_type(e1, numeric);  // true
1869     is_ex_of_type(e2, numeric);  // false
1870     is_ex_of_type(e1, add);      // false
1871     is_ex_of_type(e2, add);      // true
1872     is_ex_of_type(e1, mul);      // false
1873     is_ex_of_type(e2, mul);      // false
1874 @}
1875 @end example
1876
1877 The @code{info()} method is used for checking certain attributes of
1878 expressions. The possible values for the @code{flag} argument are defined
1879 in @file{ginac/flags.h}, the most important being explained in the following
1880 table:
1881
1882 @cartouche
1883 @multitable @columnfractions .30 .70
1884 @item @strong{Flag} @tab @strong{Returns true if the object is@dots{}}
1885 @item @code{numeric}
1886 @tab @dots{}a number (same as @code{is_ex_of_type(..., numeric)})
1887 @item @code{real}
1888 @tab @dots{}a real integer, rational or float (i.e. is not complex)
1889 @item @code{rational}
1890 @tab @dots{}an exact rational number (integers are rational, too)
1891 @item @code{integer}
1892 @tab @dots{}a (non-complex) integer
1893 @item @code{crational}
1894 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
1895 @item @code{cinteger}
1896 @tab @dots{}a (complex) integer (such as @math{2-3*I})
1897 @item @code{positive}
1898 @tab @dots{}not complex and greater than 0
1899 @item @code{negative}
1900 @tab @dots{}not complex and less than 0
1901 @item @code{nonnegative}
1902 @tab @dots{}not complex and greater than or equal to 0
1903 @item @code{posint}
1904 @tab @dots{}an integer greater than 0
1905 @item @code{negint}
1906 @tab @dots{}an integer less than 0
1907 @item @code{nonnegint}
1908 @tab @dots{}an integer greater than or equal to 0
1909 @item @code{even}
1910 @tab @dots{}an even integer
1911 @item @code{odd}
1912 @tab @dots{}an odd integer
1913 @item @code{prime}
1914 @tab @dots{}a prime integer (probabilistic primality test)
1915 @item @code{relation}
1916 @tab @dots{}a relation (same as @code{is_ex_of_type(..., relational)})
1917 @item @code{relation_equal}
1918 @tab @dots{}a @code{==} relation
1919 @item @code{relation_not_equal}
1920 @tab @dots{}a @code{!=} relation
1921 @item @code{relation_less}
1922 @tab @dots{}a @code{<} relation
1923 @item @code{relation_less_or_equal}
1924 @tab @dots{}a @code{<=} relation
1925 @item @code{relation_greater}
1926 @tab @dots{}a @code{>} relation
1927 @item @code{relation_greater_or_equal}
1928 @tab @dots{}a @code{>=} relation
1929 @item @code{symbol}
1930 @tab @dots{}a symbol (same as @code{is_ex_of_type(..., symbol)})
1931 @item @code{list}
1932 @tab @dots{}a list (same as @code{is_ex_of_type(..., lst)})
1933 @item @code{polynomial}
1934 @tab @dots{}a polynomial (i.e. only consists of sums and products of numbers and symbols with positive integer powers)
1935 @item @code{integer_polynomial}
1936 @tab @dots{}a polynomial with (non-complex) integer coefficients
1937 @item @code{cinteger_polynomial}
1938 @tab @dots{}a polynomial with (possibly complex) integer coefficients (such as @math{2-3*I})
1939 @item @code{rational_polynomial}
1940 @tab @dots{}a polynomial with (non-complex) rational coefficients
1941 @item @code{crational_polynomial}
1942 @tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I})
1943 @item @code{rational_function}
1944 @tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)})
1945 @item @code{algebraic}
1946 @tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1})
1947 @end multitable
1948 @end cartouche
1949
1950
1951 @subsection Accessing subexpressions
1952 @cindex @code{nops()}
1953 @cindex @code{op()}
1954 @cindex @code{has()}
1955 @cindex container
1956 @cindex @code{relational} (class)
1957
1958 GiNaC provides the two methods
1959
1960 @example
1961 unsigned ex::nops();
1962 ex ex::op(unsigned i);
1963 @end example
1964
1965 for accessing the subexpressions in the container-like GiNaC classes like
1966 @code{add}, @code{mul}, @code{lst}, and @code{function}. @code{nops()}
1967 determines the number of subexpressions (@samp{operands}) contained, while
1968 @code{op()} returns the @code{i}-th (0..@code{nops()-1}) subexpression.
1969 In the case of a @code{power} object, @code{op(0)} will return the basis
1970 and @code{op(1)} the exponent. For @code{indexed} objects, @code{op(0)}
1971 is the base expression and @code{op(i)}, @math{i>0} are the indices.
1972
1973 The left-hand and right-hand side expressions of objects of class
1974 @code{relational} (and only of these) can also be accessed with the methods
1975
1976 @example
1977 ex ex::lhs();
1978 ex ex::rhs();
1979 @end example
1980
1981 Finally, the method
1982
1983 @example
1984 bool ex::has(const ex & other);
1985 @end example
1986
1987 checks whether an expression contains the given subexpression @code{other}.
1988 This only works reliably if @code{other} is of an atomic class such as a
1989 @code{numeric} or a @code{symbol}. It is, e.g., not possible to verify that
1990 @code{a+b+c} contains @code{a+c} (or @code{a+b}) as a subexpression.
1991
1992
1993 @subsection Comparing expressions
1994 @cindex @code{is_equal()}
1995 @cindex @code{is_zero()}
1996
1997 Expressions can be compared with the usual C++ relational operators like
1998 @code{==}, @code{>}, and @code{<} but if the expressions contain symbols,
1999 the result is usually not determinable and the result will be @code{false},
2000 except in the case of the @code{!=} operator. You should also be aware that
2001 GiNaC will only do the most trivial test for equality (subtracting both
2002 expressions), so something like @code{(pow(x,2)+x)/x==x+1} will return
2003 @code{false}.
2004
2005 Actually, if you construct an expression like @code{a == b}, this will be
2006 represented by an object of the @code{relational} class (@xref{Relations}.)
2007 which is not evaluated until (explicitly or implicitely) cast to a @code{bool}.
2008
2009 There are also two methods
2010
2011 @example
2012 bool ex::is_equal(const ex & other);
2013 bool ex::is_zero();
2014 @end example
2015
2016 for checking whether one expression is equal to another, or equal to zero,
2017 respectively.
2018
2019 @strong{Warning:} You will also find an @code{ex::compare()} method in the
2020 GiNaC header files. This method is however only to be used internally by
2021 GiNaC to establish a canonical sort order for terms, and using it to compare
2022 expressions will give very surprising results.
2023
2024
2025 @node Substituting Symbols, Polynomial Arithmetic, Information About Expressions, Methods and Functions
2026 @c    node-name, next, previous, up
2027 @section Substituting symbols
2028 @cindex @code{subs()}
2029
2030 Symbols can be replaced with expressions via the @code{.subs()} method:
2031
2032 @example
2033 ex ex::subs(const ex & e);
2034 ex ex::subs(const lst & syms, const lst & repls);
2035 @end example
2036
2037 In the first form, @code{subs()} accepts a relational of the form
2038 @samp{symbol == expression} or a @code{lst} of such relationals. E.g.
2039
2040 @example
2041 @{
2042     symbol x("x"), y("y");
2043     ex e1 = 2*x^2-4*x+3;
2044     cout << "e1(7) = " << e1.subs(x == 7) << endl;
2045     ex e2 = x*y + x;
2046     cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
2047 @}
2048 @end example
2049
2050 will print @samp{73} and @samp{-10}, respectively.
2051
2052 If you specify multiple substitutions, they are performed in parallel, so e.g.
2053 @code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
2054
2055 The second form of @code{subs()} takes two lists, one for the symbols and
2056 one for the expressions to be substituted (both lists must contain the same
2057 number of elements). Using this form, you would write @code{subs(lst(x, y), lst(y, x))}
2058 to exchange @samp{x} and @samp{y}.
2059
2060
2061 @node Polynomial Arithmetic, Rational Expressions, Substituting Symbols, Methods and Functions
2062 @c    node-name, next, previous, up
2063 @section Polynomial arithmetic
2064
2065 @subsection Expanding and collecting
2066 @cindex @code{expand()}
2067 @cindex @code{collect()}
2068
2069 A polynomial in one or more variables has many equivalent
2070 representations.  Some useful ones serve a specific purpose.  Consider
2071 for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
2072 21*y*z + 4*z^2} (written down here in output-style).  It is equivalent
2073 to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}.  Other
2074 representations are the recursive ones where one collects for exponents
2075 in one of the three variable.  Since the factors are themselves
2076 polynomials in the remaining two variables the procedure can be
2077 repeated.  In our expample, two possibilities would be @math{(4*y + z)*x
2078 + 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
2079 x*z}.
2080
2081 To bring an expression into expanded form, its method
2082
2083 @example
2084 ex ex::expand();
2085 @end example
2086
2087 may be called.  In our example above, this corresponds to @math{4*x*y +
2088 x*z + 20*y^2 + 21*y*z + 4*z^2}.  Again, since the canonical form in
2089 GiNaC is not easily guessable you should be prepared to see different
2090 orderings of terms in such sums!
2091
2092 Another useful representation of multivariate polynomials is as a
2093 univariate polynomial in one of the variables with the coefficients
2094 being polynomials in the remaining variables.  The method
2095 @code{collect()} accomplishes this task:
2096
2097 @example
2098 ex ex::collect(const symbol & s);
2099 @end example
2100
2101 Note that the original polynomial needs to be in expanded form in order
2102 to be able to find the coefficients properly.
2103
2104 @subsection Degree and coefficients
2105 @cindex @code{degree()}
2106 @cindex @code{ldegree()}
2107 @cindex @code{coeff()}
2108
2109 The degree and low degree of a polynomial can be obtained using the two
2110 methods
2111
2112 @example
2113 int ex::degree(const symbol & s);
2114 int ex::ldegree(const symbol & s);
2115 @end example
2116
2117 which also work reliably on non-expanded input polynomials (they even work
2118 on rational functions, returning the asymptotic degree). To extract
2119 a coefficient with a certain power from an expanded polynomial you use
2120
2121 @example
2122 ex ex::coeff(const symbol & s, int n);
2123 @end example
2124
2125 You can also obtain the leading and trailing coefficients with the methods
2126
2127 @example
2128 ex ex::lcoeff(const symbol & s);
2129 ex ex::tcoeff(const symbol & s);
2130 @end example
2131
2132 which are equivalent to @code{coeff(s, degree(s))} and @code{coeff(s, ldegree(s))},
2133 respectively.
2134
2135 An application is illustrated in the next example, where a multivariate
2136 polynomial is analyzed:
2137
2138 @example
2139 #include <ginac/ginac.h>
2140 using namespace GiNaC;
2141
2142 int main()
2143 @{
2144     symbol x("x"), y("y");
2145     ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
2146                  - pow(x+y,2) + 2*pow(y+2,2) - 8;
2147     ex Poly = PolyInp.expand();
2148     
2149     for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) @{
2150         cout << "The x^" << i << "-coefficient is "
2151              << Poly.coeff(x,i) << endl;
2152     @}
2153     cout << "As polynomial in y: " 
2154          << Poly.collect(y) << endl;
2155 @}
2156 @end example
2157
2158 When run, it returns an output in the following fashion:
2159
2160 @example
2161 The x^0-coefficient is y^2+11*y
2162 The x^1-coefficient is 5*y^2-2*y
2163 The x^2-coefficient is -1
2164 The x^3-coefficient is 4*y
2165 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
2166 @end example
2167
2168 As always, the exact output may vary between different versions of GiNaC
2169 or even from run to run since the internal canonical ordering is not
2170 within the user's sphere of influence.
2171
2172
2173 @subsection Polynomial division
2174 @cindex polynomial division
2175 @cindex quotient
2176 @cindex remainder
2177 @cindex pseudo-remainder
2178 @cindex @code{quo()}
2179 @cindex @code{rem()}
2180 @cindex @code{prem()}
2181 @cindex @code{divide()}
2182
2183 The two functions
2184
2185 @example
2186 ex quo(const ex & a, const ex & b, const symbol & x);
2187 ex rem(const ex & a, const ex & b, const symbol & x);
2188 @end example
2189
2190 compute the quotient and remainder of univariate polynomials in the variable
2191 @samp{x}. The results satisfy @math{a = b*quo(a, b, x) + rem(a, b, x)}.
2192
2193 The additional function
2194
2195 @example
2196 ex prem(const ex & a, const ex & b, const symbol & x);
2197 @end example
2198
2199 computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
2200 @math{c*a = b*q + prem(a, b, x)}, where @math{c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1)}.
2201
2202 Exact division of multivariate polynomials is performed by the function
2203
2204 @example
2205 bool divide(const ex & a, const ex & b, ex & q);
2206 @end example
2207
2208 If @samp{b} divides @samp{a} over the rationals, this function returns @code{true}
2209 and returns the quotient in the variable @code{q}. Otherwise it returns @code{false}
2210 in which case the value of @code{q} is undefined.
2211
2212
2213 @subsection Unit, content and primitive part
2214 @cindex @code{unit()}
2215 @cindex @code{content()}
2216 @cindex @code{primpart()}
2217
2218 The methods
2219
2220 @example
2221 ex ex::unit(const symbol & x);
2222 ex ex::content(const symbol & x);
2223 ex ex::primpart(const symbol & x);
2224 @end example
2225
2226 return the unit part, content part, and primitive polynomial of a multivariate
2227 polynomial with respect to the variable @samp{x} (the unit part being the sign
2228 of the leading coefficient, the content part being the GCD of the coefficients,
2229 and the primitive polynomial being the input polynomial divided by the unit and
2230 content parts). The product of unit, content, and primitive part is the
2231 original polynomial.
2232
2233
2234 @subsection GCD and LCM
2235 @cindex GCD
2236 @cindex LCM
2237 @cindex @code{gcd()}
2238 @cindex @code{lcm()}
2239
2240 The functions for polynomial greatest common divisor and least common
2241 multiple have the synopsis
2242
2243 @example
2244 ex gcd(const ex & a, const ex & b);
2245 ex lcm(const ex & a, const ex & b);
2246 @end example
2247
2248 The functions @code{gcd()} and @code{lcm()} accept two expressions
2249 @code{a} and @code{b} as arguments and return a new expression, their
2250 greatest common divisor or least common multiple, respectively.  If the
2251 polynomials @code{a} and @code{b} are coprime @code{gcd(a,b)} returns 1
2252 and @code{lcm(a,b)} returns the product of @code{a} and @code{b}.
2253
2254 @example
2255 #include <ginac/ginac.h>
2256 using namespace GiNaC;
2257
2258 int main()
2259 @{
2260     symbol x("x"), y("y"), z("z");
2261     ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
2262     ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
2263
2264     ex P_gcd = gcd(P_a, P_b);
2265     // x + 5*y + 4*z
2266     ex P_lcm = lcm(P_a, P_b);
2267     // 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
2268 @}
2269 @end example
2270
2271
2272 @node Rational Expressions, Symbolic Differentiation, Polynomial Arithmetic, Methods and Functions
2273 @c    node-name, next, previous, up
2274 @section Rational expressions
2275
2276 @subsection The @code{normal} method
2277 @cindex @code{normal()}
2278 @cindex simplification
2279 @cindex temporary replacement
2280
2281 Some basic form of simplification of expressions is called for frequently.
2282 GiNaC provides the method @code{.normal()}, which converts a rational function
2283 into an equivalent rational function of the form @samp{numerator/denominator}
2284 where numerator and denominator are coprime.  If the input expression is already
2285 a fraction, it just finds the GCD of numerator and denominator and cancels it,
2286 otherwise it performs fraction addition and multiplication.
2287
2288 @code{.normal()} can also be used on expressions which are not rational functions
2289 as it will replace all non-rational objects (like functions or non-integer
2290 powers) by temporary symbols to bring the expression to the domain of rational
2291 functions before performing the normalization, and re-substituting these
2292 symbols afterwards. This algorithm is also available as a separate method
2293 @code{.to_rational()}, described below.
2294
2295 This means that both expressions @code{t1} and @code{t2} are indeed
2296 simplified in this little program:
2297
2298 @example
2299 #include <ginac/ginac.h>
2300 using namespace GiNaC;
2301
2302 int main()
2303 @{
2304     symbol x("x");
2305     ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
2306     ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
2307     cout << "t1 is " << t1.normal() << endl;
2308     cout << "t2 is " << t2.normal() << endl;
2309 @}
2310 @end example
2311
2312 Of course this works for multivariate polynomials too, so the ratio of
2313 the sample-polynomials from the section about GCD and LCM above would be
2314 normalized to @code{P_a/P_b} = @code{(4*y+z)/(y+3*z)}.
2315
2316
2317 @subsection Numerator and denominator
2318 @cindex numerator
2319 @cindex denominator
2320 @cindex @code{numer()}
2321 @cindex @code{denom()}
2322
2323 The numerator and denominator of an expression can be obtained with
2324
2325 @example
2326 ex ex::numer();
2327 ex ex::denom();
2328 @end example
2329
2330 These functions will first normalize the expression as described above and
2331 then return the numerator or denominator, respectively.
2332
2333
2334 @subsection Converting to a rational expression
2335 @cindex @code{to_rational()}
2336
2337 Some of the methods described so far only work on polynomials or rational
2338 functions. GiNaC provides a way to extend the domain of these functions to
2339 general expressions by using the temporary replacement algorithm described
2340 above. You do this by calling
2341
2342 @example
2343 ex ex::to_rational(lst &l);
2344 @end example
2345
2346 on the expression to be converted. The supplied @code{lst} will be filled
2347 with the generated temporary symbols and their replacement expressions in
2348 a format that can be used directly for the @code{subs()} method. It can also
2349 already contain a list of replacements from an earlier application of
2350 @code{.to_rational()}, so it's possible to use it on multiple expressions
2351 and get consistent results.
2352
2353 For example,
2354
2355 @example
2356 @{
2357     symbol x("x");
2358     ex a = pow(sin(x), 2) - pow(cos(x), 2);
2359     ex b = sin(x) + cos(x);
2360     ex q;
2361     lst l;
2362     divide(a.to_rational(l), b.to_rational(l), q);
2363     cout << q.subs(l) << endl;
2364 @}
2365 @end example
2366
2367 will print @samp{sin(x)-cos(x)}.
2368
2369
2370 @node Symbolic Differentiation, Series Expansion, Rational Expressions, Methods and Functions
2371 @c    node-name, next, previous, up
2372 @section Symbolic differentiation
2373 @cindex differentiation
2374 @cindex @code{diff()}
2375 @cindex chain rule
2376 @cindex product rule
2377
2378 GiNaC's objects know how to differentiate themselves.  Thus, a
2379 polynomial (class @code{add}) knows that its derivative is the sum of
2380 the derivatives of all the monomials:
2381
2382 @example
2383 #include <ginac/ginac.h>
2384 using namespace GiNaC;
2385
2386 int main()
2387 @{
2388     symbol x("x"), y("y"), z("z");
2389     ex P = pow(x, 5) + pow(x, 2) + y;
2390
2391     cout << P.diff(x,2) << endl;  // 20*x^3 + 2
2392     cout << P.diff(y) << endl;    // 1
2393     cout << P.diff(z) << endl;    // 0
2394 @}
2395 @end example
2396
2397 If a second integer parameter @var{n} is given, the @code{diff} method
2398 returns the @var{n}th derivative.
2399
2400 If @emph{every} object and every function is told what its derivative
2401 is, all derivatives of composed objects can be calculated using the
2402 chain rule and the product rule.  Consider, for instance the expression
2403 @code{1/cosh(x)}.  Since the derivative of @code{cosh(x)} is
2404 @code{sinh(x)} and the derivative of @code{pow(x,-1)} is
2405 @code{-pow(x,-2)}, GiNaC can readily compute the composition.  It turns
2406 out that the composition is the generating function for Euler Numbers,
2407 i.e. the so called @var{n}th Euler number is the coefficient of
2408 @code{x^n/n!} in the expansion of @code{1/cosh(x)}.  We may use this
2409 identity to code a function that generates Euler numbers in just three
2410 lines:
2411
2412 @cindex Euler numbers
2413 @example
2414 #include <ginac/ginac.h>
2415 using namespace GiNaC;
2416
2417 ex EulerNumber(unsigned n)
2418 @{
2419     symbol x;
2420     const ex generator = pow(cosh(x),-1);
2421     return generator.diff(x,n).subs(x==0);
2422 @}
2423
2424 int main()
2425 @{
2426     for (unsigned i=0; i<11; i+=2)
2427         cout << EulerNumber(i) << endl;
2428     return 0;
2429 @}
2430 @end example
2431
2432 When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
2433 @code{-61}, @code{1385}, @code{-50521}.  We increment the loop variable
2434 @code{i} by two since all odd Euler numbers vanish anyways.
2435
2436
2437 @node Series Expansion, Built-in Functions, Symbolic Differentiation, Methods and Functions
2438 @c    node-name, next, previous, up
2439 @section Series expansion
2440 @cindex @code{series()}
2441 @cindex Taylor expansion
2442 @cindex Laurent expansion
2443 @cindex @code{pseries} (class)
2444
2445 Expressions know how to expand themselves as a Taylor series or (more
2446 generally) a Laurent series.  As in most conventional Computer Algebra
2447 Systems, no distinction is made between those two.  There is a class of
2448 its own for storing such series (@code{class pseries}) and a built-in
2449 function (called @code{Order}) for storing the order term of the series.
2450 As a consequence, if you want to work with series, i.e. multiply two
2451 series, you need to call the method @code{ex::series} again to convert
2452 it to a series object with the usual structure (expansion plus order
2453 term).  A sample application from special relativity could read:
2454
2455 @example
2456 #include <ginac/ginac.h>
2457 using namespace GiNaC;
2458
2459 int main()
2460 @{
2461     symbol v("v"), c("c");
2462     
2463     ex gamma = 1/sqrt(1 - pow(v/c,2));
2464     ex mass_nonrel = gamma.series(v==0, 10);
2465     
2466     cout << "the relativistic mass increase with v is " << endl
2467          << mass_nonrel << endl;
2468     
2469     cout << "the inverse square of this series is " << endl
2470          << pow(mass_nonrel,-2).series(v==0, 10) << endl;
2471 @}
2472 @end example
2473
2474 Only calling the series method makes the last output simplify to
2475 @math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
2476 series raised to the power @math{-2}.
2477
2478 @cindex M@'echain's formula
2479 As another instructive application, let us calculate the numerical 
2480 value of Archimedes' constant
2481 @tex
2482 $\pi$
2483 @end tex
2484 (for which there already exists the built-in constant @code{Pi}) 
2485 using M@'echain's amazing formula
2486 @tex
2487 $\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
2488 @end tex
2489 @ifnottex
2490 @math{Pi==16*atan(1/5)-4*atan(1/239)}.
2491 @end ifnottex
2492 We may expand the arcus tangent around @code{0} and insert the fractions
2493 @code{1/5} and @code{1/239}.  But, as we have seen, a series in GiNaC
2494 carries an order term with it and the question arises what the system is
2495 supposed to do when the fractions are plugged into that order term.  The
2496 solution is to use the function @code{series_to_poly()} to simply strip
2497 the order term off:
2498
2499 @example
2500 #include <ginac/ginac.h>
2501 using namespace GiNaC;
2502
2503 ex mechain_pi(int degr)
2504 @{
2505     symbol x;
2506     ex pi_expansion = series_to_poly(atan(x).series(x,degr));
2507     ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
2508                    -4*pi_expansion.subs(x==numeric(1,239));
2509     return pi_approx;
2510 @}
2511
2512 int main()
2513 @{
2514     ex pi_frac;
2515     for (int i=2; i<12; i+=2) @{
2516         pi_frac = mechain_pi(i);
2517         cout << i << ":\t" << pi_frac << endl
2518              << "\t" << pi_frac.evalf() << endl;
2519     @}
2520     return 0;
2521 @}
2522 @end example
2523
2524 Note how we just called @code{.series(x,degr)} instead of
2525 @code{.series(x==0,degr)}.  This is a simple shortcut for @code{ex}'s
2526 method @code{series()}: if the first argument is a symbol the expression
2527 is expanded in that symbol around point @code{0}.  When you run this
2528 program, it will type out:
2529
2530 @example
2531 2:      3804/1195
2532         3.1832635983263598326
2533 4:      5359397032/1706489875
2534         3.1405970293260603143
2535 6:      38279241713339684/12184551018734375
2536         3.141621029325034425
2537 8:      76528487109180192540976/24359780855939418203125
2538         3.141591772182177295
2539 10:     327853873402258685803048818236/104359128170408663038552734375
2540         3.1415926824043995174
2541 @end example
2542
2543
2544 @node Built-in Functions, Input/Output, Series Expansion, Methods and Functions
2545 @c    node-name, next, previous, up
2546 @section Predefined mathematical functions
2547
2548 GiNaC contains the following predefined mathematical functions:
2549
2550 @cartouche
2551 @multitable @columnfractions .30 .70
2552 @item @strong{Name} @tab @strong{Function}
2553 @item @code{abs(x)}
2554 @tab absolute value
2555 @item @code{csgn(x)}
2556 @tab complex sign
2557 @item @code{sqrt(x)}
2558 @tab square root (not a GiNaC function proper but equivalent to @code{pow(x, numeric(1, 2)})
2559 @item @code{sin(x)}
2560 @tab sine
2561 @item @code{cos(x)}
2562 @tab cosine
2563 @item @code{tan(x)}
2564 @tab tangent
2565 @item @code{asin(x)}
2566 @tab inverse sine
2567 @item @code{acos(x)}
2568 @tab inverse cosine
2569 @item @code{atan(x)}
2570 @tab inverse tangent
2571 @item @code{atan2(y, x)}
2572 @tab inverse tangent with two arguments
2573 @item @code{sinh(x)}
2574 @tab hyperbolic sine
2575 @item @code{cosh(x)}
2576 @tab hyperbolic cosine
2577 @item @code{tanh(x)}
2578 @tab hyperbolic tangent
2579 @item @code{asinh(x)}
2580 @tab inverse hyperbolic sine
2581 @item @code{acosh(x)}
2582 @tab inverse hyperbolic cosine
2583 @item @code{atanh(x)}
2584 @tab inverse hyperbolic tangent
2585 @item @code{exp(x)}
2586 @tab exponential function
2587 @item @code{log(x)}
2588 @tab natural logarithm
2589 @item @code{Li2(x)}
2590 @tab Dilogarithm
2591 @item @code{zeta(x)}
2592 @tab Riemann's zeta function
2593 @item @code{zeta(n, x)}
2594 @tab derivatives of Riemann's zeta function
2595 @item @code{tgamma(x)}
2596 @tab Gamma function
2597 @item @code{lgamma(x)}
2598 @tab logarithm of Gamma function
2599 @item @code{beta(x, y)}
2600 @tab Beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
2601 @item @code{psi(x)}
2602 @tab psi (digamma) function
2603 @item @code{psi(n, x)}
2604 @tab derivatives of psi function (polygamma functions)
2605 @item @code{factorial(n)}
2606 @tab factorial function
2607 @item @code{binomial(n, m)}
2608 @tab binomial coefficients
2609 @item @code{Order(x)}
2610 @tab order term function in truncated power series
2611 @item @code{Derivative(x, l)}
2612 @tab inert partial differentiation operator (used internally)
2613 @end multitable
2614 @end cartouche
2615
2616 @cindex branch cut
2617 For functions that have a branch cut in the complex plane GiNaC follows
2618 the conventions for C++ as defined in the ANSI standard as far as
2619 possible.  In particular: the natural logarithm (@code{log}) and the
2620 square root (@code{sqrt}) both have their branch cuts running along the
2621 negative real axis where the points on the axis itself belong to the
2622 upper part (i.e. continuous with quadrant II).  The inverse
2623 trigonometric and hyperbolic functions are not defined for complex
2624 arguments by the C++ standard, however.  In GiNaC we follow the
2625 conventions used by CLN, which in turn follow the carefully designed
2626 definitions in the Common Lisp standard.  It should be noted that this
2627 convention is identical to the one used by the C99 standard and by most
2628 serious CAS.  It is to be expected that future revisions of the C++
2629 standard incorporate these functions in the complex domain in a manner
2630 compatible with C99.
2631
2632
2633 @node Input/Output, Extending GiNaC, Built-in Functions, Methods and Functions
2634 @c    node-name, next, previous, up
2635 @section Input and output of expressions
2636 @cindex I/O
2637
2638 @subsection Expression output
2639 @cindex printing
2640 @cindex output of expressions
2641
2642 The easiest way to print an expression is to write it to a stream:
2643
2644 @example
2645 @{
2646     symbol x("x");
2647     ex e = 4.5+pow(x,2)*3/2;
2648     cout << e << endl;    // prints '4.5+3/2*x^2'
2649     // ...
2650 @end example
2651
2652 The output format is identical to the @command{ginsh} input syntax and
2653 to that used by most computer algebra systems, but not directly pastable
2654 into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
2655 is printed as @samp{x^2}).
2656
2657 To print an expression in a way that can be directly used in a C or C++
2658 program, you use the method
2659
2660 @example
2661 void ex::printcsrc(ostream & os, unsigned type, const char *name);
2662 @end example
2663
2664 This outputs a line in the form of a variable definition @code{<type> <name> = <expression>}.
2665 The possible types are defined in @file{ginac/flags.h} (@code{csrc_types})
2666 and mostly affect the way in which floating point numbers are written:
2667
2668 @example
2669     // ...
2670     e.printcsrc(cout, csrc_types::ctype_float, "f");
2671     e.printcsrc(cout, csrc_types::ctype_double, "d");
2672     e.printcsrc(cout, csrc_types::ctype_cl_N, "n");
2673     // ...
2674 @end example
2675
2676 The above example will produce (note the @code{x^2} being converted to @code{x*x}):
2677
2678 @example
2679 float f = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
2680 double d = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
2681 cl_N n = (cl_F("3.0")/cl_F("2.0"))*(x*x)+cl_F("4.5");
2682 @end example
2683
2684 Finally, there are the two methods @code{printraw()} and @code{printtree()} intended for GiNaC
2685 developers, that provide a dump of the internal structure of an expression for
2686 debugging purposes:
2687
2688 @example
2689     // ...
2690     e.printraw(cout); cout << endl << endl;
2691     e.printtree(cout);
2692 @}
2693 @end example
2694
2695 produces
2696
2697 @example
2698 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))
2699
2700 type=Q25GiNaC3add, hash=0 (0x0), flags=3, nops=2
2701     power: hash=2 (0x2), flags=3
2702         x (symbol): serial=1, hash=150875740 (0x8fe2e5c), flags=11
2703         2 (numeric): hash=2147483714 (0x80000042), flags=11
2704     3/2 (numeric): hash=2147483745 (0x80000061), flags=11
2705     -----
2706     overall_coeff
2707     4.5L0 (numeric): hash=2147483723 (0x8000004b), flags=11
2708     =====
2709 @end example
2710
2711 The @code{printtree()} method is also available in @command{ginsh} as the
2712 @code{print()} function.
2713
2714
2715 @subsection Expression input
2716 @cindex input of expressions
2717
2718 GiNaC provides no way to directly read an expression from a stream because
2719 you will usually want the user to be able to enter something like @samp{2*x+sin(y)}
2720 and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
2721 @code{y} you defined in your program and there is no way to specify the
2722 desired symbols to the @code{>>} stream input operator.
2723
2724 Instead, GiNaC lets you construct an expression from a string, specifying the
2725 list of symbols to be used:
2726
2727 @example
2728 @{
2729     symbol x("x"), y("y");
2730     ex e("2*x+sin(y)", lst(x, y));
2731 @}
2732 @end example
2733
2734 The input syntax is the same as that used by @command{ginsh} and the stream
2735 output operator @code{<<}. The symbols in the string are matched by name to
2736 the symbols in the list and if GiNaC encounters a symbol not specified in
2737 the list it will throw an exception.
2738
2739 With this constructor, it's also easy to implement interactive GiNaC programs:
2740
2741 @example
2742 #include <iostream>
2743 #include <string>
2744 #include <stdexcept>
2745 #include <ginac/ginac.h>
2746 using namespace GiNaC;
2747
2748 int main()
2749 @{
2750      symbol x("x");
2751      string s;
2752
2753      cout << "Enter an expression containing 'x': ";
2754      getline(cin, s);
2755
2756      try @{
2757          ex e(s, lst(x));
2758          cout << "The derivative of " << e << " with respect to x is ";
2759          cout << e.diff(x) << ".\n";
2760      @} catch (exception &p) @{
2761          cerr << p.what() << endl;
2762      @}
2763 @}
2764 @end example
2765
2766
2767 @subsection Archiving
2768 @cindex @code{archive} (class)
2769 @cindex archiving
2770
2771 GiNaC allows creating @dfn{archives} of expressions which can be stored
2772 to or retrieved from files. To create an archive, you declare an object
2773 of class @code{archive} and archive expressions in it, giving each
2774 expression a unique name:
2775
2776 @example
2777 #include <ginac/ginac.h>
2778 #include <fstream>
2779 using namespace GiNaC;
2780
2781 int main()
2782 @{
2783     symbol x("x"), y("y"), z("z");
2784
2785     ex foo = sin(x + 2*y) + 3*z + 41;
2786     ex bar = foo + 1;
2787
2788     archive a;
2789     a.archive_ex(foo, "foo");
2790     a.archive_ex(bar, "the second one");
2791     // ...
2792 @end example
2793
2794 The archive can then be written to a file:
2795
2796 @example
2797     // ...
2798     ofstream out("foobar.gar");
2799     out << a;
2800     out.close();
2801     // ...
2802 @end example
2803
2804 The file @file{foobar.gar} contains all information that is needed to
2805 reconstruct the expressions @code{foo} and @code{bar}.
2806
2807 @cindex @command{viewgar}
2808 The tool @command{viewgar} that comes with GiNaC can be used to view
2809 the contents of GiNaC archive files:
2810
2811 @example
2812 $ viewgar foobar.gar
2813 foo = 41+sin(x+2*y)+3*z
2814 the second one = 42+sin(x+2*y)+3*z
2815 @end example
2816
2817 The point of writing archive files is of course that they can later be
2818 read in again:
2819
2820 @example
2821     // ...
2822     archive a2;
2823     ifstream in("foobar.gar");
2824     in >> a2;
2825     // ...
2826 @end example
2827
2828 And the stored expressions can be retrieved by their name:
2829
2830 @example
2831     // ...
2832     lst syms(x, y);
2833
2834     ex ex1 = a2.unarchive_ex(syms, "foo");
2835     ex ex2 = a2.unarchive_ex(syms, "the second one");
2836
2837     cout << ex1 << endl;              // prints "41+sin(x+2*y)+3*z"
2838     cout << ex2 << endl;              // prints "42+sin(x+2*y)+3*z"
2839     cout << ex1.subs(x == 2) << endl; // prints "41+sin(2+2*y)+3*z"
2840 @}
2841 @end example
2842
2843 Note that you have to supply a list of the symbols which are to be inserted
2844 in the expressions. Symbols in archives are stored by their name only and
2845 if you don't specify which symbols you have, unarchiving the expression will
2846 create new symbols with that name. E.g. if you hadn't included @code{x} in
2847 the @code{syms} list above, the @code{ex1.subs(x == 2)} statement would
2848 have had no effect because the @code{x} in @code{ex1} would have been a
2849 different symbol than the @code{x} which was defined at the beginning of
2850 the program, altough both would appear as @samp{x} when printed.
2851
2852
2853
2854 @node Extending GiNaC, What does not belong into GiNaC, Input/Output, Top
2855 @c    node-name, next, previous, up
2856 @chapter Extending GiNaC
2857
2858 By reading so far you should have gotten a fairly good understanding of
2859 GiNaC's design-patterns.  From here on you should start reading the
2860 sources.  All we can do now is issue some recommendations how to tackle
2861 GiNaC's many loose ends in order to fulfill everybody's dreams.  If you
2862 develop some useful extension please don't hesitate to contact the GiNaC
2863 authors---they will happily incorporate them into future versions.
2864
2865 @menu
2866 * What does not belong into GiNaC::  What to avoid.
2867 * Symbolic functions::               Implementing symbolic functions.
2868 * Adding classes::                   Defining new algebraic classes.
2869 @end menu
2870
2871
2872 @node What does not belong into GiNaC, Symbolic functions, Extending GiNaC, Extending GiNaC
2873 @c    node-name, next, previous, up
2874 @section What doesn't belong into GiNaC
2875
2876 @cindex @command{ginsh}
2877 First of all, GiNaC's name must be read literally.  It is designed to be
2878 a library for use within C++.  The tiny @command{ginsh} accompanying
2879 GiNaC makes this even more clear: it doesn't even attempt to provide a
2880 language.  There are no loops or conditional expressions in
2881 @command{ginsh}, it is merely a window into the library for the
2882 programmer to test stuff (or to show off).  Still, the design of a
2883 complete CAS with a language of its own, graphical capabilites and all
2884 this on top of GiNaC is possible and is without doubt a nice project for
2885 the future.
2886
2887 There are many built-in functions in GiNaC that do not know how to
2888 evaluate themselves numerically to a precision declared at runtime
2889 (using @code{Digits}).  Some may be evaluated at certain points, but not
2890 generally.  This ought to be fixed.  However, doing numerical
2891 computations with GiNaC's quite abstract classes is doomed to be
2892 inefficient.  For this purpose, the underlying foundation classes
2893 provided by @acronym{CLN} are much better suited.
2894
2895
2896 @node Symbolic functions, Adding classes, What does not belong into GiNaC, Extending GiNaC
2897 @c    node-name, next, previous, up
2898 @section Symbolic functions
2899
2900 The easiest and most instructive way to start with is probably to
2901 implement your own function.  GiNaC's functions are objects of class
2902 @code{function}.  The preprocessor is then used to convert the function
2903 names to objects with a corresponding serial number that is used
2904 internally to identify them.  You usually need not worry about this
2905 number.  New functions may be inserted into the system via a kind of
2906 `registry'.  It is your responsibility to care for some functions that
2907 are called when the user invokes certain methods.  These are usual
2908 C++-functions accepting a number of @code{ex} as arguments and returning
2909 one @code{ex}.  As an example, if we have a look at a simplified
2910 implementation of the cosine trigonometric function, we first need a
2911 function that is called when one wishes to @code{eval} it.  It could
2912 look something like this:
2913
2914 @example
2915 static ex cos_eval_method(const ex & x)
2916 @{
2917     // if (!x%(2*Pi)) return 1
2918     // if (!x%Pi) return -1
2919     // if (!x%Pi/2) return 0
2920     // care for other cases...
2921     return cos(x).hold();
2922 @}
2923 @end example
2924
2925 @cindex @code{hold()}
2926 @cindex evaluation
2927 The last line returns @code{cos(x)} if we don't know what else to do and
2928 stops a potential recursive evaluation by saying @code{.hold()}, which
2929 sets a flag to the expression signaling that it has been evaluated.  We
2930 should also implement a method for numerical evaluation and since we are
2931 lazy we sweep the problem under the rug by calling someone else's
2932 function that does so, in this case the one in class @code{numeric}:
2933
2934 @example
2935 static ex cos_evalf(const ex & x)
2936 @{
2937     return cos(ex_to_numeric(x));
2938 @}
2939 @end example
2940
2941 Differentiation will surely turn up and so we need to tell @code{cos}
2942 what the first derivative is (higher derivatives (@code{.diff(x,3)} for
2943 instance are then handled automatically by @code{basic::diff} and
2944 @code{ex::diff}):
2945
2946 @example
2947 static ex cos_deriv(const ex & x, unsigned diff_param)
2948 @{
2949     return -sin(x);
2950 @}
2951 @end example
2952
2953 @cindex product rule
2954 The second parameter is obligatory but uninteresting at this point.  It
2955 specifies which parameter to differentiate in a partial derivative in
2956 case the function has more than one parameter and its main application
2957 is for correct handling of the chain rule.  For Taylor expansion, it is
2958 enough to know how to differentiate.  But if the function you want to
2959 implement does have a pole somewhere in the complex plane, you need to
2960 write another method for Laurent expansion around that point.
2961
2962 Now that all the ingredients for @code{cos} have been set up, we need
2963 to tell the system about it.  This is done by a macro and we are not
2964 going to descibe how it expands, please consult your preprocessor if you
2965 are curious:
2966
2967 @example
2968 REGISTER_FUNCTION(cos, eval_func(cos_eval).
2969                        evalf_func(cos_evalf).
2970                        derivative_func(cos_deriv));
2971 @end example
2972
2973 The first argument is the function's name used for calling it and for
2974 output.  The second binds the corresponding methods as options to this
2975 object.  Options are separated by a dot and can be given in an arbitrary
2976 order.  GiNaC functions understand several more options which are always
2977 specified as @code{.option(params)}, for example a method for series
2978 expansion @code{.series_func(cos_series)}.  Again, if no series
2979 expansion method is given, GiNaC defaults to simple Taylor expansion,
2980 which is correct if there are no poles involved as is the case for the
2981 @code{cos} function.  The way GiNaC handles poles in case there are any
2982 is best understood by studying one of the examples, like the Gamma
2983 (@code{tgamma}) function for instance.  (In essence the function first
2984 checks if there is a pole at the evaluation point and falls back to
2985 Taylor expansion if there isn't.  Then, the pole is regularized by some
2986 suitable transformation.)  Also, the new function needs to be declared
2987 somewhere.  This may also be done by a convenient preprocessor macro:
2988
2989 @example
2990 DECLARE_FUNCTION_1P(cos)
2991 @end example
2992
2993 The suffix @code{_1P} stands for @emph{one parameter}.  Of course, this
2994 implementation of @code{cos} is very incomplete and lacks several safety
2995 mechanisms.  Please, have a look at the real implementation in GiNaC.
2996 (By the way: in case you are worrying about all the macros above we can
2997 assure you that functions are GiNaC's most macro-intense classes.  We
2998 have done our best to avoid macros where we can.)
2999
3000
3001 @node Adding classes, A Comparison With Other CAS, Symbolic functions, Extending GiNaC
3002 @c    node-name, next, previous, up
3003 @section Adding classes
3004
3005 If you are doing some very specialized things with GiNaC you may find that
3006 you have to implement your own algebraic classes to fit your needs. This
3007 section will explain how to do this by giving the example of a simple
3008 'string' class. After reading this section you will know how to properly
3009 declare a GiNaC class and what the minimum required member functions are
3010 that you have to implement. We only cover the implementation of a 'leaf'
3011 class here (i.e. one that doesn't contain subexpressions). Creating a
3012 container class like, for example, a class representing tensor products is
3013 more involved but this section should give you enough information so you can
3014 consult the source to GiNaC's predefined classes if you want to implement
3015 something more complicated.
3016
3017 @subsection GiNaC's run-time type information system
3018
3019 @cindex hierarchy of classes
3020 @cindex RTTI
3021 All algebraic classes (that is, all classes that can appear in expressions)
3022 in GiNaC are direct or indirect subclasses of the class @code{basic}. So a
3023 @code{basic *} (which is essentially what an @code{ex} is) represents a
3024 generic pointer to an algebraic class. Occasionally it is necessary to find
3025 out what the class of an object pointed to by a @code{basic *} really is.
3026 Also, for the unarchiving of expressions it must be possible to find the
3027 @code{unarchive()} function of a class given the class name (as a string). A
3028 system that provides this kind of information is called a run-time type
3029 information (RTTI) system. The C++ language provides such a thing (see the
3030 standard header file @file{<typeinfo>}) but for efficiency reasons GiNaC
3031 implements its own, simpler RTTI.
3032
3033 The RTTI in GiNaC is based on two mechanisms:
3034
3035 @itemize @bullet
3036
3037 @item
3038 The @code{basic} class declares a member variable @code{tinfo_key} which
3039 holds an unsigned integer that identifies the object's class. These numbers
3040 are defined in the @file{tinfos.h} header file for the built-in GiNaC
3041 classes. They all start with @code{TINFO_}.
3042
3043 @item
3044 By means of some clever tricks with static members, GiNaC maintains a list
3045 of information for all classes derived from @code{basic}. The information
3046 available includes the class names, the @code{tinfo_key}s, and pointers
3047 to the unarchiving functions. This class registry is defined in the
3048 @file{registrar.h} header file.
3049
3050 @end itemize
3051
3052 The disadvantage of this proprietary RTTI implementation is that there's
3053 a little more to do when implementing new classes (C++'s RTTI works more
3054 or less automatic) but don't worry, most of the work is simplified by
3055 macros.
3056
3057 @subsection A minimalistic example
3058
3059 Now we will start implementing a new class @code{mystring} that allows
3060 placing character strings in algebraic expressions (this is not very useful,
3061 but it's just an example). This class will be a direct subclass of
3062 @code{basic}. You can use this sample implementation as a starting point
3063 for your own classes.
3064
3065 The code snippets given here assume that you have included some header files
3066 as follows:
3067
3068 @example
3069 #include <iostream>
3070 #include <string>   
3071 #include <stdexcept>
3072 using namespace std;
3073
3074 #include <ginac/ginac.h>
3075 using namespace GiNaC;
3076 @end example
3077
3078 The first thing we have to do is to define a @code{tinfo_key} for our new
3079 class. This can be any arbitrary unsigned number that is not already taken
3080 by one of the existing classes but it's better to come up with something
3081 that is unlikely to clash with keys that might be added in the future. The
3082 numbers in @file{tinfos.h} are modeled somewhat after the class hierarchy
3083 which is not a requirement but we are going to stick with this scheme:
3084
3085 @example
3086 const unsigned TINFO_mystring = 0x42420001U;
3087 @end example
3088
3089 Now we can write down the class declaration. The class stores a C++
3090 @code{string} and the user shall be able to construct a @code{mystring}
3091 object from a C or C++ string:
3092
3093 @example
3094 class mystring : public basic
3095 @{
3096     GINAC_DECLARE_REGISTERED_CLASS(mystring, basic)
3097   
3098 public:
3099     mystring(const string &s);
3100     mystring(const char *s);
3101
3102 private:
3103     string str;
3104 @};
3105
3106 GIANC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
3107 @end example
3108
3109 The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
3110 macros are defined in @file{registrar.h}. They take the name of the class
3111 and its direct superclass as arguments and insert all required declarations
3112 for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
3113 the first line after the opening brace of the class definition. The
3114 @code{GINAC_IMPLEMENT_REGISTERED_CLASS} may appear anywhere else in the
3115 source (at global scope, of course, not inside a function).
3116
3117 @code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
3118 declarations of the default and copy constructor, the destructor, the
3119 assignment operator and a couple of other functions that are required. It
3120 also defines a type @code{inherited} which refers to the superclass so you
3121 don't have to modify your code every time you shuffle around the class
3122 hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} implements the copy
3123 constructor, the destructor and the assignment operator.
3124
3125 Now there are nine member functions we have to implement to get a working
3126 class:
3127
3128 @itemize
3129
3130 @item
3131 @code{mystring()}, the default constructor.
3132
3133 @item
3134 @code{void destroy(bool call_parent)}, which is used in the destructor and the
3135 assignment operator to free dynamically allocated members. The @code{call_parent}
3136 specifies whether the @code{destroy()} function of the superclass is to be
3137 called also.
3138
3139 @item
3140 @code{void copy(const mystring &other)}, which is used in the copy constructor
3141 and assignment operator to copy the member variables over from another
3142 object of the same class.
3143
3144 @item
3145 @code{void archive(archive_node &n)}, the archiving function. This stores all
3146 information needed to reconstruct an object of this class inside an
3147 @code{archive_node}.
3148
3149 @item
3150 @code{mystring(const archive_node &n, const lst &sym_lst)}, the unarchiving
3151 constructor. This constructs an instance of the class from the information
3152 found in an @code{archive_node}.
3153
3154 @item
3155 @code{ex unarchive(const archive_node &n, const lst &sym_lst)}, the static
3156 unarchiving function. It constructs a new instance by calling the unarchiving
3157 constructor.
3158
3159 @item
3160 @code{int compare_same_type(const basic &other)}, which is used internally
3161 by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
3162 -1, depending on the relative order of this object and the @code{other}
3163 object. If it returns 0, the objects are considered equal.
3164 @strong{Note:} This has nothing to do with the (numeric) ordering
3165 relationship expressed by @code{<}, @code{>=} etc (which cannot be defined
3166 for non-numeric classes). For example, @code{numeric(1).compare_same_type(numeric(2))}
3167 may return +1 even though 1 is clearly smaller than 2. Every GiNaC class
3168 must provide a @code{compare_same_type()} function, even those representing
3169 objects for which no reasonable algebraic ordering relationship can be
3170 defined.
3171
3172 @item
3173 And, of course, @code{mystring(const string &s)} and @code{mystring(const char *s)}
3174 which are the two constructors we declared.
3175
3176 @end itemize
3177
3178 Let's proceed step-by-step. The default constructor looks like this:
3179
3180 @example
3181 mystring::mystring() : inherited(TINFO_mystring)
3182 @{
3183     // dynamically allocate resources here if required
3184 @}
3185 @end example
3186
3187 The golden rule is that in all constructors you have to set the
3188 @code{tinfo_key} member to the @code{TINFO_*} value of your class. Otherwise
3189 it will be set by the constructor of the superclass and all hell will break
3190 loose in the RTTI. For your convenience, the @code{basic} class provides
3191 a constructor that takes a @code{tinfo_key} value, which we are using here
3192 (remember that in our case @code{inherited = basic}). If the superclass
3193 didn't have such a constructor, we would have to set the @code{tinfo_key}
3194 to the right value manually.
3195
3196 In the default constructor you should set all other member variables to
3197 reasonable default values (we don't need that here since our @code{str}
3198 member gets set to an empty string automatically). The constructor(s) are of
3199 course also the right place to allocate any dynamic resources you require.
3200
3201 Next, the @code{destroy()} function:
3202
3203 @example
3204 void mystring::destroy(bool call_parent)
3205 @{
3206     // free dynamically allocated resources here if required
3207     if (call_parent)
3208         inherited::destroy(call_parent);
3209 @}
3210 @end example
3211
3212 This function is where we free all dynamically allocated resources. We don't
3213 have any so we're not doing anything here, but if we had, for example, used
3214 a C-style @code{char *} to store our string, this would be the place to
3215 @code{delete[]} the string storage. If @code{call_parent} is true, we have
3216 to call the @code{destroy()} function of the superclass after we're done
3217 (to mimic C++'s automatic invocation of superclass destructors where
3218 @code{destroy()} is called from outside a destructor).
3219
3220 The @code{copy()} function just copies over the member variables from
3221 another object:
3222
3223 @example
3224 void mystring::copy(const mystring &other)
3225 @{
3226     inherited::copy(other);
3227     str = other.str;
3228 @}
3229 @end example
3230
3231 We can simply overwrite the member variables here. There's no need to worry
3232 about dynamically allocated storage. The assignment operator (which is
3233 automatically defined by @code{GINAC_IMPLEMENT_REGISTERED_CLASS}, as you
3234 recall) calls @code{destroy()} before it calls @code{copy()}. You have to
3235 explicitly call the @code{copy()} function of the superclass here so
3236 all the member variables will get copied.
3237
3238 Next are the three functions for archiving. You have to implement them even
3239 if you don't plan to use archives, but the minimum required implementation
3240 is really simple. First, the archiving function:
3241
3242 @example
3243 void mystring::archive(archive_node &n) const
3244 @{
3245     inherited::archive(n);
3246     n.add_string("string", str);
3247 @}
3248 @end example
3249
3250 The only thing that is really required is calling the @code{archive()}
3251 function of the superclass. Optionally, you can store all information you
3252 deem necessary for representing the object into the passed
3253 @code{archive_node}. We are just storing our string here. For more
3254 information on how the archiving works, consult the @file{archive.h} header
3255 file.
3256
3257 The unarchiving constructor is basically the inverse of the archiving
3258 function:
3259
3260 @example
3261 mystring::mystring(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
3262 @{
3263     n.find_string("string", str);
3264 @}
3265 @end example
3266
3267 If you don't need archiving, just leave this function empty (but you must
3268 invoke the unarchiving constructor of the superclass). Note that we don't
3269 have to set the @code{tinfo_key} here because it is done automatically
3270 by the unarchiving constructor of the @code{basic} class.
3271
3272 Finally, the unarchiving function:
3273
3274 @example
3275 ex mystring::unarchive(const archive_node &n, const lst &sym_lst)
3276 @{
3277     return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
3278 @}
3279 @end example
3280
3281 You don't have to understand how exactly this works. Just copy these four
3282 lines into your code literally (replacing the class name, of course). It
3283 calls the unarchiving constructor of the class and unless you are doing
3284 something very special (like matching @code{archive_node}s to global
3285 objects) you don't need a different implementation. For those who are
3286 interested: setting the @code{dynallocated} flag puts the object under
3287 the control of GiNaC's garbage collection. It will get deleted automatically
3288 once it is no longer referenced.
3289
3290 Our @code{compare_same_type()} function uses a provided function to compare
3291 the string members:
3292
3293 @example
3294 int mystring::compare_same_type(const basic &other) const
3295 @{
3296     const mystring &o = static_cast<const mystring &>(other);
3297     int cmpval = str.compare(o.str);
3298     if (cmpval == 0)
3299         return 0;
3300     else if (cmpval < 0)
3301         return -1;
3302     else
3303         return 1;
3304 @}
3305 @end example
3306
3307 Although this function takes a @code{basic &}, it will always be a reference
3308 to an object of exactly the same class (objects of different classes are not
3309 comparable), so the cast is safe. If this function returns 0, the two objects
3310 are considered equal (in the sense that @math{A-B=0}), so you should compare
3311 all relevant member variables.
3312
3313 Now the only thing missing is our two new constructors:
3314
3315 @example
3316 mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s)
3317 @{
3318     // dynamically allocate resources here if required
3319 @}
3320
3321 mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s)
3322 @{
3323     // dynamically allocate resources here if required
3324 @}
3325 @end example
3326
3327 No surprises here. We set the @code{str} member from the argument and
3328 remember to pass the right @code{tinfo_key} to the @code{basic} constructor.
3329
3330 That's it! We now have a minimal working GiNaC class that can store
3331 strings in algebraic expressions. Let's confirm that the RTTI works:
3332
3333 @example
3334 ex e = mystring("Hello, world!");
3335 cout << is_ex_of_type(e, mystring) << endl;
3336  // -> 1 (true)
3337
3338 cout << e.bp->class_name() << endl;
3339  // -> mystring
3340 @end example
3341
3342 Obviously it does. Let's see what the expression @code{e} looks like:
3343
3344 @example
3345 cout << e << endl;
3346  // -> [mystring object]
3347 @end example
3348
3349 Hm, not exactly what we expect, but of course the @code{mystring} class
3350 doesn't yet know how to print itself. This is done in the @code{print()}
3351 member function. Let's say that we wanted to print the string surrounded
3352 by double quotes:
3353
3354 @example
3355 class mystring : public basic
3356 @{
3357     ...
3358 public:
3359     void print(ostream &os, unsigned upper_precedence) const;
3360     ...
3361 @};
3362
3363 void mystring::print(ostream &os, unsigned upper_precedence) const
3364 @{
3365     os << '\"' << str << '\"';
3366 @}
3367 @end example
3368
3369 The @code{upper_precedence} argument is only required for container classes
3370 to correctly parenthesize the output. Let's try again to print the expression:
3371
3372 @example
3373 cout << e << endl;
3374  // -> "Hello, world!"
3375 @end example
3376
3377 Much better. The @code{mystring} class can be used in arbitrary expressions:
3378
3379 @example
3380 e += mystring("GiNaC rulez"); 
3381 cout << e << endl;
3382  // -> "GiNaC rulez"+"Hello, world!"
3383 @end example
3384
3385 (note that GiNaC's automatic term reordering is in effect here), or even
3386
3387 @example
3388 e = pow(mystring("One string"), 2*sin(Pi-mystring("Another string")));
3389 cout << e << endl;
3390  // -> "One string"^(2*sin(-"Another string"+Pi))
3391 @end example
3392
3393 Whether this makes sense is debatable but remember that this is only an
3394 example. At least it allows you to implement your own symbolic algorithms
3395 for your objects.
3396
3397 Note that GiNaC's algebraic rules remain unchanged:
3398
3399 @example
3400 e = mystring("Wow") * mystring("Wow");
3401 cout << e << endl;
3402  // -> "Wow"^2
3403
3404 e = pow(mystring("First")-mystring("Second"), 2);
3405 cout << e.expand() << endl;
3406  // -> -2*"First"*"Second"+"First"^2+"Second"^2
3407 @end example
3408
3409 There's no way to, for example, make GiNaC's @code{add} class perform string
3410 concatenation. You would have to implement this yourself.
3411
3412 @subsection Automatic evaluation
3413
3414 @cindex @code{hold()}
3415 @cindex evaluation
3416 When dealing with objects that are just a little more complicated than the
3417 simple string objects we have implemented, chances are that you will want to
3418 have some automatic simplifications or canonicalizations performed on them.
3419 This is done in the evaluation member function @code{eval()}. Let's say that
3420 we wanted all strings automatically converted to lowercase with
3421 non-alphabetic characters stripped, and empty strings removed:
3422
3423 @example
3424 class mystring : public basic
3425 @{
3426     ...
3427 public:
3428     ex eval(int level = 0) const;
3429     ...
3430 @};
3431
3432 ex mystring::eval(int level) const
3433 @{
3434     string new_str;
3435     for (int i=0; i<str.length(); i++) @{
3436         char c = str[i];
3437         if (c >= 'A' && c <= 'Z') 
3438             new_str += tolower(c);
3439         else if (c >= 'a' && c <= 'z')
3440             new_str += c;
3441     @}
3442
3443     if (new_str.length() == 0)
3444         return _ex0();
3445     else
3446         return mystring(new_str).hold();
3447 @}
3448 @end example
3449
3450 The @code{level} argument is used to limit the recursion depth of the
3451 evaluation. We don't have any subexpressions in the @code{mystring} class
3452 so we are not concerned with this. If we had, we would call the @code{eval()}
3453 functions of the subexpressions with @code{level - 1} as the argument if
3454 @code{level != 1}. The @code{hold()} member function sets a flag in the
3455 object that prevents further evaluation. Otherwise we might end up in an
3456 endless loop. When you want to return the object unmodified, use
3457 @code{return this->hold();}.
3458
3459 Let's confirm that it works:
3460
3461 @example
3462 ex e = mystring("Hello, world!") + mystring("!?#");
3463 cout << e << endl;
3464  // -> "helloworld"
3465
3466 e = mystring("Wow!") + mystring("WOW") + mystring(" W ** o ** W");  
3467 cout << e << endl;
3468  // -> 3*"wow"
3469 @end example
3470
3471 @subsection Other member functions
3472
3473 We have implemented only a small set of member functions to make the class
3474 work in the GiNaC framework. For a real algebraic class, there are probably
3475 some more functions that you will want to re-implement, such as
3476 @code{evalf()}, @code{series()} or @code{op()}. Have a look at @file{basic.h}
3477 or the header file of the class you want to make a subclass of to see
3478 what's there. You can, of course, also add your own new member functions.
3479 In this case you will probably want to define a little helper function like
3480
3481 @example
3482 inline const mystring &ex_to_mystring(const ex &e)
3483 @{
3484     return static_cast<const mystring &>(*e.bp);
3485 @}
3486 @end example
3487
3488 that let's you get at the object inside an expression (after you have verified
3489 that the type is correct) so you can call member functions that are specific
3490 to the class.
3491
3492 That's it. May the source be with you!
3493
3494
3495 @node A Comparison With Other CAS, Advantages, Adding classes, Top
3496 @c    node-name, next, previous, up
3497 @chapter A Comparison With Other CAS
3498 @cindex advocacy
3499
3500 This chapter will give you some information on how GiNaC compares to
3501 other, traditional Computer Algebra Systems, like @emph{Maple},
3502 @emph{Mathematica} or @emph{Reduce}, where it has advantages and
3503 disadvantages over these systems.
3504
3505 @menu
3506 * Advantages::                       Stengths of the GiNaC approach.
3507 * Disadvantages::                    Weaknesses of the GiNaC approach.
3508 * Why C++?::                         Attractiveness of C++.
3509 @end menu
3510
3511 @node Advantages, Disadvantages, A Comparison With Other CAS, A Comparison With Other CAS
3512 @c    node-name, next, previous, up
3513 @section Advantages
3514
3515 GiNaC has several advantages over traditional Computer
3516 Algebra Systems, like 
3517
3518 @itemize @bullet
3519
3520 @item
3521 familiar language: all common CAS implement their own proprietary
3522 grammar which you have to learn first (and maybe learn again when your
3523 vendor decides to `enhance' it).  With GiNaC you can write your program
3524 in common C++, which is standardized.
3525
3526 @cindex STL
3527 @item
3528 structured data types: you can build up structured data types using
3529 @code{struct}s or @code{class}es together with STL features instead of
3530 using unnamed lists of lists of lists.
3531
3532 @item
3533 strongly typed: in CAS, you usually have only one kind of variables
3534 which can hold contents of an arbitrary type.  This 4GL like feature is
3535 nice for novice programmers, but dangerous.
3536     
3537 @item
3538 development tools: powerful development tools exist for C++, like fancy
3539 editors (e.g. with automatic indentation and syntax highlighting),
3540 debuggers, visualization tools, documentation generators...
3541
3542 @item
3543 modularization: C++ programs can easily be split into modules by
3544 separating interface and implementation.
3545
3546 @item
3547 price: GiNaC is distributed under the GNU Public License which means
3548 that it is free and available with source code.  And there are excellent
3549 C++-compilers for free, too.
3550     
3551 @item
3552 extendable: you can add your own classes to GiNaC, thus extending it on
3553 a very low level.  Compare this to a traditional CAS that you can
3554 usually only extend on a high level by writing in the language defined
3555 by the parser.  In particular, it turns out to be almost impossible to
3556 fix bugs in a traditional system.
3557
3558 @item
3559 multiple interfaces: Though real GiNaC programs have to be written in
3560 some editor, then be compiled, linked and executed, there are more ways
3561 to work with the GiNaC engine.  Many people want to play with
3562 expressions interactively, as in traditional CASs.  Currently, two such
3563 windows into GiNaC have been implemented and many more are possible: the
3564 tiny @command{ginsh} that is part of the distribution exposes GiNaC's
3565 types to a command line and second, as a more consistent approach, an
3566 interactive interface to the @acronym{Cint} C++ interpreter has been put
3567 together (called @acronym{GiNaC-cint}) that allows an interactive
3568 scripting interface consistent with the C++ language.
3569
3570 @item
3571 seemless integration: it is somewhere between difficult and impossible
3572 to call CAS functions from within a program written in C++ or any other
3573 programming language and vice versa.  With GiNaC, your symbolic routines
3574 are part of your program.  You can easily call third party libraries,
3575 e.g. for numerical evaluation or graphical interaction.  All other
3576 approaches are much more cumbersome: they range from simply ignoring the
3577 problem (i.e. @emph{Maple}) to providing a method for `embedding' the
3578 system (i.e. @emph{Yacas}).
3579
3580 @item
3581 efficiency: often large parts of a program do not need symbolic
3582 calculations at all.  Why use large integers for loop variables or
3583 arbitrary precision arithmetics where @code{int} and @code{double} are
3584 sufficient?  For pure symbolic applications, GiNaC is comparable in
3585 speed with other CAS.
3586
3587 @end itemize
3588
3589
3590 @node Disadvantages, Why C++?, Advantages, A Comparison With Other CAS
3591 @c    node-name, next, previous, up
3592 @section Disadvantages
3593
3594 Of course it also has some disadvantages:
3595
3596 @itemize @bullet
3597
3598 @item
3599 advanced features: GiNaC cannot compete with a program like
3600 @emph{Reduce} which exists for more than 30 years now or @emph{Maple}
3601 which grows since 1981 by the work of dozens of programmers, with
3602 respect to mathematical features.  Integration, factorization,
3603 non-trivial simplifications, limits etc. are missing in GiNaC (and are
3604 not planned for the near future).
3605
3606 @item
3607 portability: While the GiNaC library itself is designed to avoid any
3608 platform dependent features (it should compile on any ANSI compliant C++
3609 compiler), the currently used version of the CLN library (fast large
3610 integer and arbitrary precision arithmetics) can be compiled only on
3611 systems with a recently new C++ compiler from the GNU Compiler
3612 Collection (@acronym{GCC}).@footnote{This is because CLN uses
3613 PROVIDE/REQUIRE like macros to let the compiler gather all static
3614 initializations, which works for GNU C++ only.}  GiNaC uses recent
3615 language features like explicit constructors, mutable members, RTTI,
3616 @code{dynamic_cast}s and STL, so ANSI compliance is meant literally.
3617 Recent @acronym{GCC} versions starting at 2.95, although itself not yet
3618 ANSI compliant, support all needed features.
3619     
3620 @end itemize
3621
3622
3623 @node Why C++?, Internal Structures, Disadvantages, A Comparison With Other CAS
3624 @c    node-name, next, previous, up
3625 @section Why C++?
3626
3627 Why did we choose to implement GiNaC in C++ instead of Java or any other
3628 language?  C++ is not perfect: type checking is not strict (casting is
3629 possible), separation between interface and implementation is not
3630 complete, object oriented design is not enforced.  The main reason is
3631 the often scolded feature of operator overloading in C++.  While it may
3632 be true that operating on classes with a @code{+} operator is rarely
3633 meaningful, it is perfectly suited for algebraic expressions.  Writing
3634 @math{3x+5y} as @code{3*x+5*y} instead of
3635 @code{x.times(3).plus(y.times(5))} looks much more natural.
3636 Furthermore, the main developers are more familiar with C++ than with
3637 any other programming language.
3638
3639
3640 @node Internal Structures, Expressions are reference counted, Why C++? , Top
3641 @c    node-name, next, previous, up
3642 @appendix Internal Structures
3643
3644 @menu
3645 * Expressions are reference counted::
3646 * Internal representation of products and sums::
3647 @end menu
3648
3649 @node Expressions are reference counted, Internal representation of products and sums, Internal Structures, Internal Structures
3650 @c    node-name, next, previous, up
3651 @appendixsection Expressions are reference counted
3652
3653 @cindex reference counting
3654 @cindex copy-on-write
3655 @cindex garbage collection
3656 An expression is extremely light-weight since internally it works like a
3657 handle to the actual representation and really holds nothing more than a
3658 pointer to some other object. What this means in practice is that
3659 whenever you create two @code{ex} and set the second equal to the first
3660 no copying process is involved. Instead, the copying takes place as soon
3661 as you try to change the second.  Consider the simple sequence of code:
3662
3663 @example
3664 #include <ginac/ginac.h>
3665 using namespace GiNaC;
3666
3667 int main()
3668 @{
3669     symbol x("x"), y("y"), z("z");
3670     ex e1, e2;
3671
3672     e1 = sin(x + 2*y) + 3*z + 41;
3673     e2 = e1;                // e2 points to same object as e1
3674     cout << e2 << endl;     // prints sin(x+2*y)+3*z+41
3675     e2 += 1;                // e2 is copied into a new object
3676     cout << e2 << endl;     // prints sin(x+2*y)+3*z+42
3677 @}
3678 @end example
3679
3680 The line @code{e2 = e1;} creates a second expression pointing to the
3681 object held already by @code{e1}.  The time involved for this operation
3682 is therefore constant, no matter how large @code{e1} was.  Actual
3683 copying, however, must take place in the line @code{e2 += 1;} because
3684 @code{e1} and @code{e2} are not handles for the same object any more.
3685 This concept is called @dfn{copy-on-write semantics}.  It increases
3686 performance considerably whenever one object occurs multiple times and
3687 represents a simple garbage collection scheme because when an @code{ex}
3688 runs out of scope its destructor checks whether other expressions handle
3689 the object it points to too and deletes the object from memory if that
3690 turns out not to be the case.  A slightly less trivial example of
3691 differentiation using the chain-rule should make clear how powerful this
3692 can be:
3693
3694 @example
3695 #include <ginac/ginac.h>
3696 using namespace GiNaC;
3697
3698 int main()
3699 @{
3700     symbol x("x"), y("y");
3701
3702     ex e1 = x + 3*y;
3703     ex e2 = pow(e1, 3);
3704     ex e3 = diff(sin(e2), x);   // first derivative of sin(e2) by x
3705     cout << e1 << endl          // prints x+3*y
3706          << e2 << endl          // prints (x+3*y)^3
3707          << e3 << endl;         // prints 3*(x+3*y)^2*cos((x+3*y)^3)
3708 @}
3709 @end example
3710
3711 Here, @code{e1} will actually be referenced three times while @code{e2}
3712 will be referenced two times.  When the power of an expression is built,
3713 that expression needs not be copied.  Likewise, since the derivative of
3714 a power of an expression can be easily expressed in terms of that
3715 expression, no copying of @code{e1} is involved when @code{e3} is
3716 constructed.  So, when @code{e3} is constructed it will print as
3717 @code{3*(x+3*y)^2*cos((x+3*y)^3)} but the argument of @code{cos()} only
3718 holds a reference to @code{e2} and the factor in front is just
3719 @code{3*e1^2}.
3720
3721 As a user of GiNaC, you cannot see this mechanism of copy-on-write
3722 semantics.  When you insert an expression into a second expression, the