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