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