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