1 \input texinfo @c -*-texinfo-*-
3 @setfilename ginac.info
4 @settitle GiNaC, an open framework for symbolic computation within the C++ programming language
11 @c I hate putting "@noindent" in front of every paragraph.
19 * ginac: (ginac). C++ library for symbolic computation.
23 This is a tutorial that documents GiNaC @value{VERSION}, an open
24 framework for symbolic computation within the C++ programming language.
26 Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany
28 Permission is granted to make and distribute verbatim copies of
29 this manual provided the copyright notice and this permission notice
30 are preserved on all copies.
33 Permission is granted to process this file through TeX and print the
34 results, provided the printed document carries copying permission
35 notice identical to this one except for the removal of this paragraph
38 Permission is granted to copy and distribute modified versions of this
39 manual under the conditions for verbatim copying, provided that the entire
40 resulting derived work is distributed under the terms of a permission
41 notice identical to this one.
45 @c finalout prevents ugly black rectangles on overfull hbox lines
47 @title GiNaC @value{VERSION}
48 @subtitle An open framework for symbolic computation within the C++ programming language
49 @subtitle @value{UPDATED}
50 @author @uref{http://www.ginac.de}
53 @vskip 0pt plus 1filll
54 Copyright @copyright{} 1999-2007 Johannes Gutenberg University Mainz, Germany
56 Permission is granted to make and distribute verbatim copies of
57 this manual provided the copyright notice and this permission notice
58 are preserved on all copies.
60 Permission is granted to copy and distribute modified versions of this
61 manual under the conditions for verbatim copying, provided that the entire
62 resulting derived work is distributed under the terms of a permission
63 notice identical to this one.
72 @node Top, Introduction, (dir), (dir)
73 @c node-name, next, previous, up
76 This is a tutorial that documents GiNaC @value{VERSION}, an open
77 framework for symbolic computation within the C++ programming language.
80 * Introduction:: GiNaC's purpose.
81 * A tour of GiNaC:: A quick tour of the library.
82 * Installation:: How to install the package.
83 * Basic concepts:: Description of fundamental classes.
84 * Methods and functions:: Algorithms for symbolic manipulations.
85 * Extending GiNaC:: How to extend the library.
86 * A comparison with other CAS:: Compares GiNaC to traditional CAS.
87 * Internal structures:: Description of some internal structures.
88 * Package tools:: Configuring packages to work with GiNaC.
94 @node Introduction, A tour of GiNaC, Top, Top
95 @c node-name, next, previous, up
97 @cindex history of GiNaC
99 The motivation behind GiNaC derives from the observation that most
100 present day computer algebra systems (CAS) are linguistically and
101 semantically impoverished. Although they are quite powerful tools for
102 learning math and solving particular problems they lack modern
103 linguistic structures that allow for the creation of large-scale
104 projects. GiNaC is an attempt to overcome this situation by extending a
105 well established and standardized computer language (C++) by some
106 fundamental symbolic capabilities, thus allowing for integrated systems
107 that embed symbolic manipulations together with more established areas
108 of computer science (like computation-intense numeric applications,
109 graphical interfaces, etc.) under one roof.
111 The particular problem that led to the writing of the GiNaC framework is
112 still a very active field of research, namely the calculation of higher
113 order corrections to elementary particle interactions. There,
114 theoretical physicists are interested in matching present day theories
115 against experiments taking place at particle accelerators. The
116 computations involved are so complex they call for a combined symbolical
117 and numerical approach. This turned out to be quite difficult to
118 accomplish with the present day CAS we have worked with so far and so we
119 tried to fill the gap by writing GiNaC. But of course its applications
120 are in no way restricted to theoretical physics.
122 This tutorial is intended for the novice user who is new to GiNaC but
123 already has some background in C++ programming. However, since a
124 hand-made documentation like this one is difficult to keep in sync with
125 the development, the actual documentation is inside the sources in the
126 form of comments. That documentation may be parsed by one of the many
127 Javadoc-like documentation systems. If you fail at generating it you
128 may access it from @uref{http://www.ginac.de/reference/, the GiNaC home
129 page}. It is an invaluable resource not only for the advanced user who
130 wishes to extend the system (or chase bugs) but for everybody who wants
131 to comprehend the inner workings of GiNaC. This little tutorial on the
132 other hand only covers the basic things that are unlikely to change in
136 The GiNaC framework for symbolic computation within the C++ programming
137 language is Copyright @copyright{} 1999-2007 Johannes Gutenberg
138 University Mainz, Germany.
140 This program is free software; you can redistribute it and/or
141 modify it under the terms of the GNU General Public License as
142 published by the Free Software Foundation; either version 2 of the
143 License, or (at your option) any later version.
145 This program is distributed in the hope that it will be useful, but
146 WITHOUT ANY WARRANTY; without even the implied warranty of
147 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
148 General Public License for more details.
150 You should have received a copy of the GNU General Public License
151 along with this program; see the file COPYING. If not, write to the
152 Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
156 @node A tour of GiNaC, How to use it from within C++, Introduction, Top
157 @c node-name, next, previous, up
158 @chapter A Tour of GiNaC
160 This quick tour of GiNaC wants to arise your interest in the
161 subsequent chapters by showing off a bit. Please excuse us if it
162 leaves many open questions.
165 * How to use it from within C++:: Two simple examples.
166 * What it can do for you:: A Tour of GiNaC's features.
170 @node How to use it from within C++, What it can do for you, A tour of GiNaC, A tour of GiNaC
171 @c node-name, next, previous, up
172 @section How to use it from within C++
174 The GiNaC open framework for symbolic computation within the C++ programming
175 language does not try to define a language of its own as conventional
176 CAS do. Instead, it extends the capabilities of C++ by symbolic
177 manipulations. Here is how to generate and print a simple (and rather
178 pointless) bivariate polynomial with some large coefficients:
182 #include <ginac/ginac.h>
184 using namespace GiNaC;
188 symbol x("x"), y("y");
191 for (int i=0; i<3; ++i)
192 poly += factorial(i+16)*pow(x,i)*pow(y,2-i);
194 cout << poly << endl;
199 Assuming the file is called @file{hello.cc}, on our system we can compile
200 and run it like this:
203 $ c++ hello.cc -o hello -lcln -lginac
205 355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
208 (@xref{Package tools}, for tools that help you when creating a software
209 package that uses GiNaC.)
211 @cindex Hermite polynomial
212 Next, there is a more meaningful C++ program that calls a function which
213 generates Hermite polynomials in a specified free variable.
217 #include <ginac/ginac.h>
219 using namespace GiNaC;
221 ex HermitePoly(const symbol & x, int n)
223 ex HKer=exp(-pow(x, 2));
224 // uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2)
225 return normal(pow(-1, n) * diff(HKer, x, n) / HKer);
232 for (int i=0; i<6; ++i)
233 cout << "H_" << i << "(z) == " << HermitePoly(z,i) << endl;
239 When run, this will type out
245 H_3(z) == -12*z+8*z^3
246 H_4(z) == -48*z^2+16*z^4+12
247 H_5(z) == 120*z-160*z^3+32*z^5
250 This method of generating the coefficients is of course far from optimal
251 for production purposes.
253 In order to show some more examples of what GiNaC can do we will now use
254 the @command{ginsh}, a simple GiNaC interactive shell that provides a
255 convenient window into GiNaC's capabilities.
258 @node What it can do for you, Installation, How to use it from within C++, A tour of GiNaC
259 @c node-name, next, previous, up
260 @section What it can do for you
262 @cindex @command{ginsh}
263 After invoking @command{ginsh} one can test and experiment with GiNaC's
264 features much like in other Computer Algebra Systems except that it does
265 not provide programming constructs like loops or conditionals. For a
266 concise description of the @command{ginsh} syntax we refer to its
267 accompanied man page. Suffice to say that assignments and comparisons in
268 @command{ginsh} are written as they are in C, i.e. @code{=} assigns and
271 It can manipulate arbitrary precision integers in a very fast way.
272 Rational numbers are automatically converted to fractions of coprime
277 369988485035126972924700782451696644186473100389722973815184405301748249
279 123329495011708990974900260817232214728824366796574324605061468433916083
286 Exact numbers are always retained as exact numbers and only evaluated as
287 floating point numbers if requested. For instance, with numeric
288 radicals is dealt pretty much as with symbols. Products of sums of them
292 > expand((1+a^(1/5)-a^(2/5))^3);
293 1+3*a+3*a^(1/5)-5*a^(3/5)-a^(6/5)
294 > expand((1+3^(1/5)-3^(2/5))^3);
296 > evalf((1+3^(1/5)-3^(2/5))^3);
297 0.33408977534118624228
300 The function @code{evalf} that was used above converts any number in
301 GiNaC's expressions into floating point numbers. This can be done to
302 arbitrary predefined accuracy:
306 0.14285714285714285714
310 0.1428571428571428571428571428571428571428571428571428571428571428571428
311 5714285714285714285714285714285714285
314 Exact numbers other than rationals that can be manipulated in GiNaC
315 include predefined constants like Archimedes' @code{Pi}. They can both
316 be used in symbolic manipulations (as an exact number) as well as in
317 numeric expressions (as an inexact number):
323 9.869604401089358619+x
327 11.869604401089358619
330 Built-in functions evaluate immediately to exact numbers if
331 this is possible. Conversions that can be safely performed are done
332 immediately; conversions that are not generally valid are not done:
343 (Note that converting the last input to @code{x} would allow one to
344 conclude that @code{42*Pi} is equal to @code{0}.)
346 Linear equation systems can be solved along with basic linear
347 algebra manipulations over symbolic expressions. In C++ GiNaC offers
348 a matrix class for this purpose but we can see what it can do using
349 @command{ginsh}'s bracket notation to type them in:
352 > lsolve(a+x*y==z,x);
354 > lsolve(@{3*x+5*y == 7, -2*x+10*y == -5@}, @{x, y@});
356 > M = [ [1, 3], [-3, 2] ];
360 > charpoly(M,lambda);
362 > A = [ [1, 1], [2, -1] ];
365 [[1,1],[2,-1]]+2*[[1,3],[-3,2]]
368 > B = [ [0, 0, a], [b, 1, -b], [-1/a, 0, 0] ];
369 > evalm(B^(2^12345));
370 [[1,0,0],[0,1,0],[0,0,1]]
373 Multivariate polynomials and rational functions may be expanded,
374 collected and normalized (i.e. converted to a ratio of two coprime
378 > a = x^4 + 2*x^2*y^2 + 4*x^3*y + 12*x*y^3 - 3*y^4;
379 12*x*y^3+2*x^2*y^2+4*x^3*y-3*y^4+x^4
380 > b = x^2 + 4*x*y - y^2;
383 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 4*x^3*y-y^2-3*y^4+(12*y^3+4*y)*x+x^4+x^2*(1+2*y^2)
387 12*x*y^3-3*y^4+(-1+2*x^2)*y^2+(4*x+4*x^3)*y+x^2+x^4
392 You can differentiate functions and expand them as Taylor or Laurent
393 series in a very natural syntax (the second argument of @code{series} is
394 a relation defining the evaluation point, the third specifies the
397 @cindex Zeta function
401 > series(sin(x),x==0,4);
403 > series(1/tan(x),x==0,4);
404 x^(-1)-1/3*x+Order(x^2)
405 > series(tgamma(x),x==0,3);
406 x^(-1)-Euler+(1/12*Pi^2+1/2*Euler^2)*x+
407 (-1/3*zeta(3)-1/12*Pi^2*Euler-1/6*Euler^3)*x^2+Order(x^3)
409 x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
410 -(0.90747907608088628905)*x^2+Order(x^3)
411 > series(tgamma(2*sin(x)-2),x==Pi/2,6);
412 -(x-1/2*Pi)^(-2)+(-1/12*Pi^2-1/2*Euler^2-1/240)*(x-1/2*Pi)^2
413 -Euler-1/12+Order((x-1/2*Pi)^3)
416 Here we have made use of the @command{ginsh}-command @code{%} to pop the
417 previously evaluated element from @command{ginsh}'s internal stack.
419 Often, functions don't have roots in closed form. Nevertheless, it's
420 quite easy to compute a solution numerically, to arbitrary precision:
425 > fsolve(cos(x)==x,x,0,2);
426 0.7390851332151606416553120876738734040134117589007574649658
428 > X=fsolve(f,x,-10,10);
429 2.2191071489137460325957851882042901681753665565320678854155
431 -6.372367644529809108115521591070847222364418220770475144296E-58
434 Notice how the final result above differs slightly from zero by about
435 @math{6*10^(-58)}. This is because with 50 decimal digits precision the
436 root cannot be represented more accurately than @code{X}. Such
437 inaccuracies are to be expected when computing with finite floating
440 If you ever wanted to convert units in C or C++ and found this is
441 cumbersome, here is the solution. Symbolic types can always be used as
442 tags for different types of objects. Converting from wrong units to the
443 metric system is now easy:
451 140613.91592783185568*kg*m^(-2)
455 @node Installation, Prerequisites, What it can do for you, Top
456 @c node-name, next, previous, up
457 @chapter Installation
460 GiNaC's installation follows the spirit of most GNU software. It is
461 easily installed on your system by three steps: configuration, build,
465 * Prerequisites:: Packages upon which GiNaC depends.
466 * Configuration:: How to configure GiNaC.
467 * Building GiNaC:: How to compile GiNaC.
468 * Installing GiNaC:: How to install GiNaC on your system.
472 @node Prerequisites, Configuration, Installation, Installation
473 @c node-name, next, previous, up
474 @section Prerequisites
476 In order to install GiNaC on your system, some prerequisites need to be
477 met. First of all, you need to have a C++-compiler adhering to the
478 ANSI-standard @cite{ISO/IEC 14882:1998(E)}. We used GCC for development
479 so if you have a different compiler you are on your own. For the
480 configuration to succeed you need a Posix compliant shell installed in
481 @file{/bin/sh}, GNU @command{bash} is fine. Perl is needed by the built
482 process as well, since some of the source files are automatically
483 generated by Perl scripts. Last but not least, the CLN library
484 is used extensively and needs to be installed on your system.
485 Please get it from @uref{ftp://ftpthep.physik.uni-mainz.de/pub/gnu/}
486 (it is covered by GPL) and install it prior to trying to install
487 GiNaC. The configure script checks if it can find it and if it cannot
488 it will refuse to continue.
491 @node Configuration, Building GiNaC, Prerequisites, Installation
492 @c node-name, next, previous, up
493 @section Configuration
494 @cindex configuration
497 To configure GiNaC means to prepare the source distribution for
498 building. It is done via a shell script called @command{configure} that
499 is shipped with the sources and was originally generated by GNU
500 Autoconf. Since a configure script generated by GNU Autoconf never
501 prompts, all customization must be done either via command line
502 parameters or environment variables. It accepts a list of parameters,
503 the complete set of which can be listed by calling it with the
504 @option{--help} option. The most important ones will be shortly
505 described in what follows:
510 @option{--disable-shared}: When given, this option switches off the
511 build of a shared library, i.e. a @file{.so} file. This may be convenient
512 when developing because it considerably speeds up compilation.
515 @option{--prefix=@var{PREFIX}}: The directory where the compiled library
516 and headers are installed. It defaults to @file{/usr/local} which means
517 that the library is installed in the directory @file{/usr/local/lib},
518 the header files in @file{/usr/local/include/ginac} and the documentation
519 (like this one) into @file{/usr/local/share/doc/GiNaC}.
522 @option{--libdir=@var{LIBDIR}}: Use this option in case you want to have
523 the library installed in some other directory than
524 @file{@var{PREFIX}/lib/}.
527 @option{--includedir=@var{INCLUDEDIR}}: Use this option in case you want
528 to have the header files installed in some other directory than
529 @file{@var{PREFIX}/include/ginac/}. For instance, if you specify
530 @option{--includedir=/usr/include} you will end up with the header files
531 sitting in the directory @file{/usr/include/ginac/}. Note that the
532 subdirectory @file{ginac} is enforced by this process in order to
533 keep the header files separated from others. This avoids some
534 clashes and allows for an easier deinstallation of GiNaC. This ought
535 to be considered A Good Thing (tm).
538 @option{--datadir=@var{DATADIR}}: This option may be given in case you
539 want to have the documentation installed in some other directory than
540 @file{@var{PREFIX}/share/doc/GiNaC/}.
544 In addition, you may specify some environment variables. @env{CXX}
545 holds the path and the name of the C++ compiler in case you want to
546 override the default in your path. (The @command{configure} script
547 searches your path for @command{c++}, @command{g++}, @command{gcc},
548 @command{CC}, @command{cxx} and @command{cc++} in that order.) It may
549 be very useful to define some compiler flags with the @env{CXXFLAGS}
550 environment variable, like optimization, debugging information and
551 warning levels. If omitted, it defaults to @option{-g
552 -O2}.@footnote{The @command{configure} script is itself generated from
553 the file @file{configure.ac}. It is only distributed in packaged
554 releases of GiNaC. If you got the naked sources, e.g. from CVS, you
555 must generate @command{configure} along with the various
556 @file{Makefile.in} by using the @command{autogen.sh} script. This will
557 require a fair amount of support from your local toolchain, though.}
559 The whole process is illustrated in the following two
560 examples. (Substitute @command{setenv @var{VARIABLE} @var{value}} for
561 @command{export @var{VARIABLE}=@var{value}} if the Berkeley C shell is
564 Here is a simple configuration for a site-wide GiNaC library assuming
565 everything is in default paths:
568 $ export CXXFLAGS="-Wall -O2"
572 And here is a configuration for a private static GiNaC library with
573 several components sitting in custom places (site-wide GCC and private
574 CLN). The compiler is persuaded to be picky and full assertions and
575 debugging information are switched on:
578 $ export CXX=/usr/local/gnu/bin/c++
579 $ export CPPFLAGS="$(CPPFLAGS) -I$(HOME)/include"
580 $ export CXXFLAGS="$(CXXFLAGS) -DDO_GINAC_ASSERT -ggdb -Wall -pedantic"
581 $ export LDFLAGS="$(LDFLAGS) -L$(HOME)/lib"
582 $ ./configure --disable-shared --prefix=$(HOME)
586 @node Building GiNaC, Installing GiNaC, Configuration, Installation
587 @c node-name, next, previous, up
588 @section Building GiNaC
589 @cindex building GiNaC
591 After proper configuration you should just build the whole
596 at the command prompt and go for a cup of coffee. The exact time it
597 takes to compile GiNaC depends not only on the speed of your machines
598 but also on other parameters, for instance what value for @env{CXXFLAGS}
599 you entered. Optimization may be very time-consuming.
601 Just to make sure GiNaC works properly you may run a collection of
602 regression tests by typing
608 This will compile some sample programs, run them and check the output
609 for correctness. The regression tests fall in three categories. First,
610 the so called @emph{exams} are performed, simple tests where some
611 predefined input is evaluated (like a pupils' exam). Second, the
612 @emph{checks} test the coherence of results among each other with
613 possible random input. Third, some @emph{timings} are performed, which
614 benchmark some predefined problems with different sizes and display the
615 CPU time used in seconds. Each individual test should return a message
616 @samp{passed}. This is mostly intended to be a QA-check if something
617 was broken during development, not a sanity check of your system. Some
618 of the tests in sections @emph{checks} and @emph{timings} may require
619 insane amounts of memory and CPU time. Feel free to kill them if your
620 machine catches fire. Another quite important intent is to allow people
621 to fiddle around with optimization.
623 By default, the only documentation that will be built is this tutorial
624 in @file{.info} format. To build the GiNaC tutorial and reference manual
625 in HTML, DVI, PostScript, or PDF formats, use one of
634 Generally, the top-level Makefile runs recursively to the
635 subdirectories. It is therefore safe to go into any subdirectory
636 (@code{doc/}, @code{ginsh/}, @dots{}) and simply type @code{make}
637 @var{target} there in case something went wrong.
640 @node Installing GiNaC, Basic concepts, Building GiNaC, Installation
641 @c node-name, next, previous, up
642 @section Installing GiNaC
645 To install GiNaC on your system, simply type
651 As described in the section about configuration the files will be
652 installed in the following directories (the directories will be created
653 if they don't already exist):
658 @file{libginac.a} will go into @file{@var{PREFIX}/lib/} (or
659 @file{@var{LIBDIR}}) which defaults to @file{/usr/local/lib/}.
660 So will @file{libginac.so} unless the configure script was
661 given the option @option{--disable-shared}. The proper symlinks
662 will be established as well.
665 All the header files will be installed into @file{@var{PREFIX}/include/ginac/}
666 (or @file{@var{INCLUDEDIR}/ginac/}, if specified).
669 All documentation (info) will be stuffed into
670 @file{@var{PREFIX}/share/doc/GiNaC/} (or
671 @file{@var{DATADIR}/doc/GiNaC/}, if @var{DATADIR} was specified).
675 For the sake of completeness we will list some other useful make
676 targets: @command{make clean} deletes all files generated by
677 @command{make}, i.e. all the object files. In addition @command{make
678 distclean} removes all files generated by the configuration and
679 @command{make maintainer-clean} goes one step further and deletes files
680 that may require special tools to rebuild (like the @command{libtool}
681 for instance). Finally @command{make uninstall} removes the installed
682 library, header files and documentation@footnote{Uninstallation does not
683 work after you have called @command{make distclean} since the
684 @file{Makefile} is itself generated by the configuration from
685 @file{Makefile.in} and hence deleted by @command{make distclean}. There
686 are two obvious ways out of this dilemma. First, you can run the
687 configuration again with the same @var{PREFIX} thus creating a
688 @file{Makefile} with a working @samp{uninstall} target. Second, you can
689 do it by hand since you now know where all the files went during
693 @node Basic concepts, Expressions, Installing GiNaC, Top
694 @c node-name, next, previous, up
695 @chapter Basic concepts
697 This chapter will describe the different fundamental objects that can be
698 handled by GiNaC. But before doing so, it is worthwhile introducing you
699 to the more commonly used class of expressions, representing a flexible
700 meta-class for storing all mathematical objects.
703 * Expressions:: The fundamental GiNaC class.
704 * Automatic evaluation:: Evaluation and canonicalization.
705 * Error handling:: How the library reports errors.
706 * The class hierarchy:: Overview of GiNaC's classes.
707 * Symbols:: Symbolic objects.
708 * Numbers:: Numerical objects.
709 * Constants:: Pre-defined constants.
710 * Fundamental containers:: Sums, products and powers.
711 * Lists:: Lists of expressions.
712 * Mathematical functions:: Mathematical functions.
713 * Relations:: Equality, Inequality and all that.
714 * Integrals:: Symbolic integrals.
715 * Matrices:: Matrices.
716 * Indexed objects:: Handling indexed quantities.
717 * Non-commutative objects:: Algebras with non-commutative products.
718 * Hash maps:: A faster alternative to std::map<>.
722 @node Expressions, Automatic evaluation, Basic concepts, Basic concepts
723 @c node-name, next, previous, up
725 @cindex expression (class @code{ex})
728 The most common class of objects a user deals with is the expression
729 @code{ex}, representing a mathematical object like a variable, number,
730 function, sum, product, etc@dots{} Expressions may be put together to form
731 new expressions, passed as arguments to functions, and so on. Here is a
732 little collection of valid expressions:
735 ex MyEx1 = 5; // simple number
736 ex MyEx2 = x + 2*y; // polynomial in x and y
737 ex MyEx3 = (x + 1)/(x - 1); // rational expression
738 ex MyEx4 = sin(x + 2*y) + 3*z + 41; // containing a function
739 ex MyEx5 = MyEx4 + 1; // similar to above
742 Expressions are handles to other more fundamental objects, that often
743 contain other expressions thus creating a tree of expressions
744 (@xref{Internal structures}, for particular examples). Most methods on
745 @code{ex} therefore run top-down through such an expression tree. For
746 example, the method @code{has()} scans recursively for occurrences of
747 something inside an expression. Thus, if you have declared @code{MyEx4}
748 as in the example above @code{MyEx4.has(y)} will find @code{y} inside
749 the argument of @code{sin} and hence return @code{true}.
751 The next sections will outline the general picture of GiNaC's class
752 hierarchy and describe the classes of objects that are handled by
755 @subsection Note: Expressions and STL containers
757 GiNaC expressions (@code{ex} objects) have value semantics (they can be
758 assigned, reassigned and copied like integral types) but the operator
759 @code{<} doesn't provide a well-defined ordering on them. In STL-speak,
760 expressions are @samp{Assignable} but not @samp{LessThanComparable}.
762 This implies that in order to use expressions in sorted containers such as
763 @code{std::map<>} and @code{std::set<>} you have to supply a suitable
764 comparison predicate. GiNaC provides such a predicate, called
765 @code{ex_is_less}. For example, a set of expressions should be defined
766 as @code{std::set<ex, ex_is_less>}.
768 Unsorted containers such as @code{std::vector<>} and @code{std::list<>}
769 don't pose a problem. A @code{std::vector<ex>} works as expected.
771 @xref{Information about expressions}, for more about comparing and ordering
775 @node Automatic evaluation, Error handling, Expressions, Basic concepts
776 @c node-name, next, previous, up
777 @section Automatic evaluation and canonicalization of expressions
780 GiNaC performs some automatic transformations on expressions, to simplify
781 them and put them into a canonical form. Some examples:
784 ex MyEx1 = 2*x - 1 + x; // 3*x-1
785 ex MyEx2 = x - x; // 0
786 ex MyEx3 = cos(2*Pi); // 1
787 ex MyEx4 = x*y/x; // y
790 This behavior is usually referred to as @dfn{automatic} or @dfn{anonymous
791 evaluation}. GiNaC only performs transformations that are
795 at most of complexity
803 algebraically correct, possibly except for a set of measure zero (e.g.
804 @math{x/x} is transformed to @math{1} although this is incorrect for @math{x=0})
807 There are two types of automatic transformations in GiNaC that may not
808 behave in an entirely obvious way at first glance:
812 The terms of sums and products (and some other things like the arguments of
813 symmetric functions, the indices of symmetric tensors etc.) are re-ordered
814 into a canonical form that is deterministic, but not lexicographical or in
815 any other way easy to guess (it almost always depends on the number and
816 order of the symbols you define). However, constructing the same expression
817 twice, either implicitly or explicitly, will always result in the same
820 Expressions of the form 'number times sum' are automatically expanded (this
821 has to do with GiNaC's internal representation of sums and products). For
824 ex MyEx5 = 2*(x + y); // 2*x+2*y
825 ex MyEx6 = z*(x + y); // z*(x+y)
829 The general rule is that when you construct expressions, GiNaC automatically
830 creates them in canonical form, which might differ from the form you typed in
831 your program. This may create some awkward looking output (@samp{-y+x} instead
832 of @samp{x-y}) but allows for more efficient operation and usually yields
833 some immediate simplifications.
835 @cindex @code{eval()}
836 Internally, the anonymous evaluator in GiNaC is implemented by the methods
839 ex ex::eval(int level = 0) const;
840 ex basic::eval(int level = 0) const;
843 but unless you are extending GiNaC with your own classes or functions, there
844 should never be any reason to call them explicitly. All GiNaC methods that
845 transform expressions, like @code{subs()} or @code{normal()}, automatically
846 re-evaluate their results.
849 @node Error handling, The class hierarchy, Automatic evaluation, Basic concepts
850 @c node-name, next, previous, up
851 @section Error handling
853 @cindex @code{pole_error} (class)
855 GiNaC reports run-time errors by throwing C++ exceptions. All exceptions
856 generated by GiNaC are subclassed from the standard @code{exception} class
857 defined in the @file{<stdexcept>} header. In addition to the predefined
858 @code{logic_error}, @code{domain_error}, @code{out_of_range},
859 @code{invalid_argument}, @code{runtime_error}, @code{range_error} and
860 @code{overflow_error} types, GiNaC also defines a @code{pole_error}
861 exception that gets thrown when trying to evaluate a mathematical function
864 The @code{pole_error} class has a member function
867 int pole_error::degree() const;
870 that returns the order of the singularity (or 0 when the pole is
871 logarithmic or the order is undefined).
873 When using GiNaC it is useful to arrange for exceptions to be caught in
874 the main program even if you don't want to do any special error handling.
875 Otherwise whenever an error occurs in GiNaC, it will be delegated to the
876 default exception handler of your C++ compiler's run-time system which
877 usually only aborts the program without giving any information what went
880 Here is an example for a @code{main()} function that catches and prints
881 exceptions generated by GiNaC:
886 #include <ginac/ginac.h>
888 using namespace GiNaC;
896 @} catch (exception &p) @{
897 cerr << p.what() << endl;
905 @node The class hierarchy, Symbols, Error handling, Basic concepts
906 @c node-name, next, previous, up
907 @section The class hierarchy
909 GiNaC's class hierarchy consists of several classes representing
910 mathematical objects, all of which (except for @code{ex} and some
911 helpers) are internally derived from one abstract base class called
912 @code{basic}. You do not have to deal with objects of class
913 @code{basic}, instead you'll be dealing with symbols, numbers,
914 containers of expressions and so on.
918 To get an idea about what kinds of symbolic composites may be built we
919 have a look at the most important classes in the class hierarchy and
920 some of the relations among the classes:
922 @image{classhierarchy}
924 The abstract classes shown here (the ones without drop-shadow) are of no
925 interest for the user. They are used internally in order to avoid code
926 duplication if two or more classes derived from them share certain
927 features. An example is @code{expairseq}, a container for a sequence of
928 pairs each consisting of one expression and a number (@code{numeric}).
929 What @emph{is} visible to the user are the derived classes @code{add}
930 and @code{mul}, representing sums and products. @xref{Internal
931 structures}, where these two classes are described in more detail. The
932 following table shortly summarizes what kinds of mathematical objects
933 are stored in the different classes:
936 @multitable @columnfractions .22 .78
937 @item @code{symbol} @tab Algebraic symbols @math{a}, @math{x}, @math{y}@dots{}
938 @item @code{constant} @tab Constants like
945 @item @code{numeric} @tab All kinds of numbers, @math{42}, @math{7/3*I}, @math{3.14159}@dots{}
946 @item @code{add} @tab Sums like @math{x+y} or @math{a-(2*b)+3}
947 @item @code{mul} @tab Products like @math{x*y} or @math{2*a^2*(x+y+z)/b}
948 @item @code{ncmul} @tab Products of non-commutative objects
949 @item @code{power} @tab Exponentials such as @math{x^2}, @math{a^b},
954 @code{sqrt(}@math{2}@code{)}
957 @item @code{pseries} @tab Power Series, e.g. @math{x-1/6*x^3+1/120*x^5+O(x^7)}
958 @item @code{function} @tab A symbolic function like
965 @item @code{lst} @tab Lists of expressions @{@math{x}, @math{2*y}, @math{3+z}@}
966 @item @code{matrix} @tab @math{m}x@math{n} matrices of expressions
967 @item @code{relational} @tab A relation like the identity @math{x}@code{==}@math{y}
968 @item @code{indexed} @tab Indexed object like @math{A_ij}
969 @item @code{tensor} @tab Special tensor like the delta and metric tensors
970 @item @code{idx} @tab Index of an indexed object
971 @item @code{varidx} @tab Index with variance
972 @item @code{spinidx} @tab Index with variance and dot (used in Weyl-van-der-Waerden spinor formalism)
973 @item @code{wildcard} @tab Wildcard for pattern matching
974 @item @code{structure} @tab Template for user-defined classes
979 @node Symbols, Numbers, The class hierarchy, Basic concepts
980 @c node-name, next, previous, up
982 @cindex @code{symbol} (class)
983 @cindex hierarchy of classes
986 Symbolic indeterminates, or @dfn{symbols} for short, are for symbolic
987 manipulation what atoms are for chemistry.
989 A typical symbol definition looks like this:
994 This definition actually contains three very different things:
996 @item a C++ variable named @code{x}
997 @item a @code{symbol} object stored in this C++ variable; this object
998 represents the symbol in a GiNaC expression
999 @item the string @code{"x"} which is the name of the symbol, used (almost)
1000 exclusively for printing expressions holding the symbol
1003 Symbols have an explicit name, supplied as a string during construction,
1004 because in C++, variable names can't be used as values, and the C++ compiler
1005 throws them away during compilation.
1007 It is possible to omit the symbol name in the definition:
1012 In this case, GiNaC will assign the symbol an internal, unique name of the
1013 form @code{symbolNNN}. This won't affect the usability of the symbol but
1014 the output of your calculations will become more readable if you give your
1015 symbols sensible names (for intermediate expressions that are only used
1016 internally such anonymous symbols can be quite useful, however).
1018 Now, here is one important property of GiNaC that differentiates it from
1019 other computer algebra programs you may have used: GiNaC does @emph{not} use
1020 the names of symbols to tell them apart, but a (hidden) serial number that
1021 is unique for each newly created @code{symbol} object. In you want to use
1022 one and the same symbol in different places in your program, you must only
1023 create one @code{symbol} object and pass that around. If you create another
1024 symbol, even if it has the same name, GiNaC will treat it as a different
1041 // prints "x^6" which looks right, but...
1043 cout << e.degree(x) << endl;
1044 // ...this doesn't work. The symbol "x" here is different from the one
1045 // in f() and in the expression returned by f(). Consequently, it
1050 One possibility to ensure that @code{f()} and @code{main()} use the same
1051 symbol is to pass the symbol as an argument to @code{f()}:
1053 ex f(int n, const ex & x)
1062 // Now, f() uses the same symbol.
1065 cout << e.degree(x) << endl;
1066 // prints "6", as expected
1070 Another possibility would be to define a global symbol @code{x} that is used
1071 by both @code{f()} and @code{main()}. If you are using global symbols and
1072 multiple compilation units you must take special care, however. Suppose
1073 that you have a header file @file{globals.h} in your program that defines
1074 a @code{symbol x("x");}. In this case, every unit that includes
1075 @file{globals.h} would also get its own definition of @code{x} (because
1076 header files are just inlined into the source code by the C++ preprocessor),
1077 and hence you would again end up with multiple equally-named, but different,
1078 symbols. Instead, the @file{globals.h} header should only contain a
1079 @emph{declaration} like @code{extern symbol x;}, with the definition of
1080 @code{x} moved into a C++ source file such as @file{globals.cpp}.
1082 A different approach to ensuring that symbols used in different parts of
1083 your program are identical is to create them with a @emph{factory} function
1086 const symbol & get_symbol(const string & s)
1088 static map<string, symbol> directory;
1089 map<string, symbol>::iterator i = directory.find(s);
1090 if (i != directory.end())
1093 return directory.insert(make_pair(s, symbol(s))).first->second;
1097 This function returns one newly constructed symbol for each name that is
1098 passed in, and it returns the same symbol when called multiple times with
1099 the same name. Using this symbol factory, we can rewrite our example like
1104 return pow(get_symbol("x"), n);
1111 // Both calls of get_symbol("x") yield the same symbol.
1112 cout << e.degree(get_symbol("x")) << endl;
1117 Instead of creating symbols from strings we could also have
1118 @code{get_symbol()} take, for example, an integer number as its argument.
1119 In this case, we would probably want to give the generated symbols names
1120 that include this number, which can be accomplished with the help of an
1121 @code{ostringstream}.
1123 In general, if you're getting weird results from GiNaC such as an expression
1124 @samp{x-x} that is not simplified to zero, you should check your symbol
1127 As we said, the names of symbols primarily serve for purposes of expression
1128 output. But there are actually two instances where GiNaC uses the names for
1129 identifying symbols: When constructing an expression from a string, and when
1130 recreating an expression from an archive (@pxref{Input/output}).
1132 In addition to its name, a symbol may contain a special string that is used
1135 symbol x("x", "\\Box");
1138 This creates a symbol that is printed as "@code{x}" in normal output, but
1139 as "@code{\Box}" in LaTeX code (@xref{Input/output}, for more
1140 information about the different output formats of expressions in GiNaC).
1141 GiNaC automatically creates proper LaTeX code for symbols having names of
1142 greek letters (@samp{alpha}, @samp{mu}, etc.).
1144 @cindex @code{subs()}
1145 Symbols in GiNaC can't be assigned values. If you need to store results of
1146 calculations and give them a name, use C++ variables of type @code{ex}.
1147 If you want to replace a symbol in an expression with something else, you
1148 can invoke the expression's @code{.subs()} method
1149 (@pxref{Substituting expressions}).
1151 @cindex @code{realsymbol()}
1152 By default, symbols are expected to stand in for complex values, i.e. they live
1153 in the complex domain. As a consequence, operations like complex conjugation,
1154 for example (@pxref{Complex expressions}), do @emph{not} evaluate if applied
1155 to such symbols. Likewise @code{log(exp(x))} does not evaluate to @code{x},
1156 because of the unknown imaginary part of @code{x}.
1157 On the other hand, if you are sure that your symbols will hold only real
1158 values, you would like to have such functions evaluated. Therefore GiNaC
1159 allows you to specify
1160 the domain of the symbol. Instead of @code{symbol x("x");} you can write
1161 @code{realsymbol x("x");} to tell GiNaC that @code{x} stands in for real values.
1163 @cindex @code{possymbol()}
1164 Furthermore, it is also possible to declare a symbol as positive. This will,
1165 for instance, enable the automatic simplification of @code{abs(x)} into
1166 @code{x}. This is done by declaying the symbol as @code{possymbol x("x");}.
1169 @node Numbers, Constants, Symbols, Basic concepts
1170 @c node-name, next, previous, up
1172 @cindex @code{numeric} (class)
1178 For storing numerical things, GiNaC uses Bruno Haible's library CLN.
1179 The classes therein serve as foundation classes for GiNaC. CLN stands
1180 for Class Library for Numbers or alternatively for Common Lisp Numbers.
1181 In order to find out more about CLN's internals, the reader is referred to
1182 the documentation of that library. @inforef{Introduction, , cln}, for
1183 more information. Suffice to say that it is by itself build on top of
1184 another library, the GNU Multiple Precision library GMP, which is an
1185 extremely fast library for arbitrary long integers and rationals as well
1186 as arbitrary precision floating point numbers. It is very commonly used
1187 by several popular cryptographic applications. CLN extends GMP by
1188 several useful things: First, it introduces the complex number field
1189 over either reals (i.e. floating point numbers with arbitrary precision)
1190 or rationals. Second, it automatically converts rationals to integers
1191 if the denominator is unity and complex numbers to real numbers if the
1192 imaginary part vanishes and also correctly treats algebraic functions.
1193 Third it provides good implementations of state-of-the-art algorithms
1194 for all trigonometric and hyperbolic functions as well as for
1195 calculation of some useful constants.
1197 The user can construct an object of class @code{numeric} in several
1198 ways. The following example shows the four most important constructors.
1199 It uses construction from C-integer, construction of fractions from two
1200 integers, construction from C-float and construction from a string:
1204 #include <ginac/ginac.h>
1205 using namespace GiNaC;
1209 numeric two = 2; // exact integer 2
1210 numeric r(2,3); // exact fraction 2/3
1211 numeric e(2.71828); // floating point number
1212 numeric p = "3.14159265358979323846"; // constructor from string
1213 // Trott's constant in scientific notation:
1214 numeric trott("1.0841015122311136151E-2");
1216 std::cout << two*p << std::endl; // floating point 6.283...
1221 @cindex complex numbers
1222 The imaginary unit in GiNaC is a predefined @code{numeric} object with the
1227 numeric z1 = 2-3*I; // exact complex number 2-3i
1228 numeric z2 = 5.9+1.6*I; // complex floating point number
1232 It may be tempting to construct fractions by writing @code{numeric r(3/2)}.
1233 This would, however, call C's built-in operator @code{/} for integers
1234 first and result in a numeric holding a plain integer 1. @strong{Never
1235 use the operator @code{/} on integers} unless you know exactly what you
1236 are doing! Use the constructor from two integers instead, as shown in
1237 the example above. Writing @code{numeric(1)/2} may look funny but works
1240 @cindex @code{Digits}
1242 We have seen now the distinction between exact numbers and floating
1243 point numbers. Clearly, the user should never have to worry about
1244 dynamically created exact numbers, since their `exactness' always
1245 determines how they ought to be handled, i.e. how `long' they are. The
1246 situation is different for floating point numbers. Their accuracy is
1247 controlled by one @emph{global} variable, called @code{Digits}. (For
1248 those readers who know about Maple: it behaves very much like Maple's
1249 @code{Digits}). All objects of class numeric that are constructed from
1250 then on will be stored with a precision matching that number of decimal
1255 #include <ginac/ginac.h>
1256 using namespace std;
1257 using namespace GiNaC;
1261 numeric three(3.0), one(1.0);
1262 numeric x = one/three;
1264 cout << "in " << Digits << " digits:" << endl;
1266 cout << Pi.evalf() << endl;
1278 The above example prints the following output to screen:
1282 0.33333333333333333334
1283 3.1415926535897932385
1285 0.33333333333333333333333333333333333333333333333333333333333333333334
1286 3.1415926535897932384626433832795028841971693993751058209749445923078
1290 Note that the last number is not necessarily rounded as you would
1291 naively expect it to be rounded in the decimal system. But note also,
1292 that in both cases you got a couple of extra digits. This is because
1293 numbers are internally stored by CLN as chunks of binary digits in order
1294 to match your machine's word size and to not waste precision. Thus, on
1295 architectures with different word size, the above output might even
1296 differ with regard to actually computed digits.
1298 It should be clear that objects of class @code{numeric} should be used
1299 for constructing numbers or for doing arithmetic with them. The objects
1300 one deals with most of the time are the polymorphic expressions @code{ex}.
1302 @subsection Tests on numbers
1304 Once you have declared some numbers, assigned them to expressions and
1305 done some arithmetic with them it is frequently desired to retrieve some
1306 kind of information from them like asking whether that number is
1307 integer, rational, real or complex. For those cases GiNaC provides
1308 several useful methods. (Internally, they fall back to invocations of
1309 certain CLN functions.)
1311 As an example, let's construct some rational number, multiply it with
1312 some multiple of its denominator and test what comes out:
1316 #include <ginac/ginac.h>
1317 using namespace std;
1318 using namespace GiNaC;
1320 // some very important constants:
1321 const numeric twentyone(21);
1322 const numeric ten(10);
1323 const numeric five(5);
1327 numeric answer = twentyone;
1330 cout << answer.is_integer() << endl; // false, it's 21/5
1332 cout << answer.is_integer() << endl; // true, it's 42 now!
1336 Note that the variable @code{answer} is constructed here as an integer
1337 by @code{numeric}'s copy constructor but in an intermediate step it
1338 holds a rational number represented as integer numerator and integer
1339 denominator. When multiplied by 10, the denominator becomes unity and
1340 the result is automatically converted to a pure integer again.
1341 Internally, the underlying CLN is responsible for this behavior and we
1342 refer the reader to CLN's documentation. Suffice to say that
1343 the same behavior applies to complex numbers as well as return values of
1344 certain functions. Complex numbers are automatically converted to real
1345 numbers if the imaginary part becomes zero. The full set of tests that
1346 can be applied is listed in the following table.
1349 @multitable @columnfractions .30 .70
1350 @item @strong{Method} @tab @strong{Returns true if the object is@dots{}}
1351 @item @code{.is_zero()}
1352 @tab @dots{}equal to zero
1353 @item @code{.is_positive()}
1354 @tab @dots{}not complex and greater than 0
1355 @item @code{.is_integer()}
1356 @tab @dots{}a (non-complex) integer
1357 @item @code{.is_pos_integer()}
1358 @tab @dots{}an integer and greater than 0
1359 @item @code{.is_nonneg_integer()}
1360 @tab @dots{}an integer and greater equal 0
1361 @item @code{.is_even()}
1362 @tab @dots{}an even integer
1363 @item @code{.is_odd()}
1364 @tab @dots{}an odd integer
1365 @item @code{.is_prime()}
1366 @tab @dots{}a prime integer (probabilistic primality test)
1367 @item @code{.is_rational()}
1368 @tab @dots{}an exact rational number (integers are rational, too)
1369 @item @code{.is_real()}
1370 @tab @dots{}a real integer, rational or float (i.e. is not complex)
1371 @item @code{.is_cinteger()}
1372 @tab @dots{}a (complex) integer (such as @math{2-3*I})
1373 @item @code{.is_crational()}
1374 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
1378 @subsection Numeric functions
1380 The following functions can be applied to @code{numeric} objects and will be
1381 evaluated immediately:
1384 @multitable @columnfractions .30 .70
1385 @item @strong{Name} @tab @strong{Function}
1386 @item @code{inverse(z)}
1387 @tab returns @math{1/z}
1388 @cindex @code{inverse()} (numeric)
1389 @item @code{pow(a, b)}
1390 @tab exponentiation @math{a^b}
1393 @item @code{real(z)}
1395 @cindex @code{real()}
1396 @item @code{imag(z)}
1398 @cindex @code{imag()}
1399 @item @code{csgn(z)}
1400 @tab complex sign (returns an @code{int})
1401 @item @code{step(x)}
1402 @tab step function (returns an @code{numeric})
1403 @item @code{numer(z)}
1404 @tab numerator of rational or complex rational number
1405 @item @code{denom(z)}
1406 @tab denominator of rational or complex rational number
1407 @item @code{sqrt(z)}
1409 @item @code{isqrt(n)}
1410 @tab integer square root
1411 @cindex @code{isqrt()}
1418 @item @code{asin(z)}
1420 @item @code{acos(z)}
1422 @item @code{atan(z)}
1423 @tab inverse tangent
1424 @item @code{atan(y, x)}
1425 @tab inverse tangent with two arguments
1426 @item @code{sinh(z)}
1427 @tab hyperbolic sine
1428 @item @code{cosh(z)}
1429 @tab hyperbolic cosine
1430 @item @code{tanh(z)}
1431 @tab hyperbolic tangent
1432 @item @code{asinh(z)}
1433 @tab inverse hyperbolic sine
1434 @item @code{acosh(z)}
1435 @tab inverse hyperbolic cosine
1436 @item @code{atanh(z)}
1437 @tab inverse hyperbolic tangent
1439 @tab exponential function
1441 @tab natural logarithm
1444 @item @code{zeta(z)}
1445 @tab Riemann's zeta function
1446 @item @code{tgamma(z)}
1448 @item @code{lgamma(z)}
1449 @tab logarithm of gamma function
1451 @tab psi (digamma) function
1452 @item @code{psi(n, z)}
1453 @tab derivatives of psi function (polygamma functions)
1454 @item @code{factorial(n)}
1455 @tab factorial function @math{n!}
1456 @item @code{doublefactorial(n)}
1457 @tab double factorial function @math{n!!}
1458 @cindex @code{doublefactorial()}
1459 @item @code{binomial(n, k)}
1460 @tab binomial coefficients
1461 @item @code{bernoulli(n)}
1462 @tab Bernoulli numbers
1463 @cindex @code{bernoulli()}
1464 @item @code{fibonacci(n)}
1465 @tab Fibonacci numbers
1466 @cindex @code{fibonacci()}
1467 @item @code{mod(a, b)}
1468 @tab modulus in positive representation (in the range @code{[0, abs(b)-1]} with the sign of b, or zero)
1469 @cindex @code{mod()}
1470 @item @code{smod(a, b)}
1471 @tab modulus in symmetric representation (in the range @code{[-iquo(abs(b)-1, 2), iquo(abs(b), 2)]})
1472 @cindex @code{smod()}
1473 @item @code{irem(a, b)}
1474 @tab integer remainder (has the sign of @math{a}, or is zero)
1475 @cindex @code{irem()}
1476 @item @code{irem(a, b, q)}
1477 @tab integer remainder and quotient, @code{irem(a, b, q) == a-q*b}
1478 @item @code{iquo(a, b)}
1479 @tab integer quotient
1480 @cindex @code{iquo()}
1481 @item @code{iquo(a, b, r)}
1482 @tab integer quotient and remainder, @code{r == a-iquo(a, b)*b}
1483 @item @code{gcd(a, b)}
1484 @tab greatest common divisor
1485 @item @code{lcm(a, b)}
1486 @tab least common multiple
1490 Most of these functions are also available as symbolic functions that can be
1491 used in expressions (@pxref{Mathematical functions}) or, like @code{gcd()},
1492 as polynomial algorithms.
1494 @subsection Converting numbers
1496 Sometimes it is desirable to convert a @code{numeric} object back to a
1497 built-in arithmetic type (@code{int}, @code{double}, etc.). The @code{numeric}
1498 class provides a couple of methods for this purpose:
1500 @cindex @code{to_int()}
1501 @cindex @code{to_long()}
1502 @cindex @code{to_double()}
1503 @cindex @code{to_cl_N()}
1505 int numeric::to_int() const;
1506 long numeric::to_long() const;
1507 double numeric::to_double() const;
1508 cln::cl_N numeric::to_cl_N() const;
1511 @code{to_int()} and @code{to_long()} only work when the number they are
1512 applied on is an exact integer. Otherwise the program will halt with a
1513 message like @samp{Not a 32-bit integer}. @code{to_double()} applied on a
1514 rational number will return a floating-point approximation. Both
1515 @code{to_int()/to_long()} and @code{to_double()} discard the imaginary
1516 part of complex numbers.
1519 @node Constants, Fundamental containers, Numbers, Basic concepts
1520 @c node-name, next, previous, up
1522 @cindex @code{constant} (class)
1525 @cindex @code{Catalan}
1526 @cindex @code{Euler}
1527 @cindex @code{evalf()}
1528 Constants behave pretty much like symbols except that they return some
1529 specific number when the method @code{.evalf()} is called.
1531 The predefined known constants are:
1534 @multitable @columnfractions .14 .30 .56
1535 @item @strong{Name} @tab @strong{Common Name} @tab @strong{Numerical Value (to 35 digits)}
1537 @tab Archimedes' constant
1538 @tab 3.14159265358979323846264338327950288
1539 @item @code{Catalan}
1540 @tab Catalan's constant
1541 @tab 0.91596559417721901505460351493238411
1543 @tab Euler's (or Euler-Mascheroni) constant
1544 @tab 0.57721566490153286060651209008240243
1549 @node Fundamental containers, Lists, Constants, Basic concepts
1550 @c node-name, next, previous, up
1551 @section Sums, products and powers
1555 @cindex @code{power}
1557 Simple rational expressions are written down in GiNaC pretty much like
1558 in other CAS or like expressions involving numerical variables in C.
1559 The necessary operators @code{+}, @code{-}, @code{*} and @code{/} have
1560 been overloaded to achieve this goal. When you run the following
1561 code snippet, the constructor for an object of type @code{mul} is
1562 automatically called to hold the product of @code{a} and @code{b} and
1563 then the constructor for an object of type @code{add} is called to hold
1564 the sum of that @code{mul} object and the number one:
1568 symbol a("a"), b("b");
1573 @cindex @code{pow()}
1574 For exponentiation, you have already seen the somewhat clumsy (though C-ish)
1575 statement @code{pow(x,2);} to represent @code{x} squared. This direct
1576 construction is necessary since we cannot safely overload the constructor
1577 @code{^} in C++ to construct a @code{power} object. If we did, it would
1578 have several counterintuitive and undesired effects:
1582 Due to C's operator precedence, @code{2*x^2} would be parsed as @code{(2*x)^2}.
1584 Due to the binding of the operator @code{^}, @code{x^a^b} would result in
1585 @code{(x^a)^b}. This would be confusing since most (though not all) other CAS
1586 interpret this as @code{x^(a^b)}.
1588 Also, expressions involving integer exponents are very frequently used,
1589 which makes it even more dangerous to overload @code{^} since it is then
1590 hard to distinguish between the semantics as exponentiation and the one
1591 for exclusive or. (It would be embarrassing to return @code{1} where one
1592 has requested @code{2^3}.)
1595 @cindex @command{ginsh}
1596 All effects are contrary to mathematical notation and differ from the
1597 way most other CAS handle exponentiation, therefore overloading @code{^}
1598 is ruled out for GiNaC's C++ part. The situation is different in
1599 @command{ginsh}, there the exponentiation-@code{^} exists. (Also note
1600 that the other frequently used exponentiation operator @code{**} does
1601 not exist at all in C++).
1603 To be somewhat more precise, objects of the three classes described
1604 here, are all containers for other expressions. An object of class
1605 @code{power} is best viewed as a container with two slots, one for the
1606 basis, one for the exponent. All valid GiNaC expressions can be
1607 inserted. However, basic transformations like simplifying
1608 @code{pow(pow(x,2),3)} to @code{x^6} automatically are only performed
1609 when this is mathematically possible. If we replace the outer exponent
1610 three in the example by some symbols @code{a}, the simplification is not
1611 safe and will not be performed, since @code{a} might be @code{1/2} and
1614 Objects of type @code{add} and @code{mul} are containers with an
1615 arbitrary number of slots for expressions to be inserted. Again, simple
1616 and safe simplifications are carried out like transforming
1617 @code{3*x+4-x} to @code{2*x+4}.
1620 @node Lists, Mathematical functions, Fundamental containers, Basic concepts
1621 @c node-name, next, previous, up
1622 @section Lists of expressions
1623 @cindex @code{lst} (class)
1625 @cindex @code{nops()}
1627 @cindex @code{append()}
1628 @cindex @code{prepend()}
1629 @cindex @code{remove_first()}
1630 @cindex @code{remove_last()}
1631 @cindex @code{remove_all()}
1633 The GiNaC class @code{lst} serves for holding a @dfn{list} of arbitrary
1634 expressions. They are not as ubiquitous as in many other computer algebra
1635 packages, but are sometimes used to supply a variable number of arguments of
1636 the same type to GiNaC methods such as @code{subs()} and some @code{matrix}
1637 constructors, so you should have a basic understanding of them.
1639 Lists can be constructed by assigning a comma-separated sequence of
1644 symbol x("x"), y("y");
1647 // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y',
1652 There are also constructors that allow direct creation of lists of up to
1653 16 expressions, which is often more convenient but slightly less efficient:
1657 // This produces the same list 'l' as above:
1658 // lst l(x, 2, y, x+y);
1659 // lst l = lst(x, 2, y, x+y);
1663 Use the @code{nops()} method to determine the size (number of expressions) of
1664 a list and the @code{op()} method or the @code{[]} operator to access
1665 individual elements:
1669 cout << l.nops() << endl; // prints '4'
1670 cout << l.op(2) << " " << l[0] << endl; // prints 'y x'
1674 As with the standard @code{list<T>} container, accessing random elements of a
1675 @code{lst} is generally an operation of order @math{O(N)}. Faster read-only
1676 sequential access to the elements of a list is possible with the
1677 iterator types provided by the @code{lst} class:
1680 typedef ... lst::const_iterator;
1681 typedef ... lst::const_reverse_iterator;
1682 lst::const_iterator lst::begin() const;
1683 lst::const_iterator lst::end() const;
1684 lst::const_reverse_iterator lst::rbegin() const;
1685 lst::const_reverse_iterator lst::rend() const;
1688 For example, to print the elements of a list individually you can use:
1693 for (lst::const_iterator i = l.begin(); i != l.end(); ++i)
1698 which is one order faster than
1703 for (size_t i = 0; i < l.nops(); ++i)
1704 cout << l.op(i) << endl;
1708 These iterators also allow you to use some of the algorithms provided by
1709 the C++ standard library:
1713 // print the elements of the list (requires #include <iterator>)
1714 std::copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));
1716 // sum up the elements of the list (requires #include <numeric>)
1717 ex sum = std::accumulate(l.begin(), l.end(), ex(0));
1718 cout << sum << endl; // prints '2+2*x+2*y'
1722 @code{lst} is one of the few GiNaC classes that allow in-place modifications
1723 (the only other one is @code{matrix}). You can modify single elements:
1727 l[1] = 42; // l is now @{x, 42, y, x+y@}
1728 l.let_op(1) = 7; // l is now @{x, 7, y, x+y@}
1732 You can append or prepend an expression to a list with the @code{append()}
1733 and @code{prepend()} methods:
1737 l.append(4*x); // l is now @{x, 7, y, x+y, 4*x@}
1738 l.prepend(0); // l is now @{0, x, 7, y, x+y, 4*x@}
1742 You can remove the first or last element of a list with @code{remove_first()}
1743 and @code{remove_last()}:
1747 l.remove_first(); // l is now @{x, 7, y, x+y, 4*x@}
1748 l.remove_last(); // l is now @{x, 7, y, x+y@}
1752 You can remove all the elements of a list with @code{remove_all()}:
1756 l.remove_all(); // l is now empty
1760 You can bring the elements of a list into a canonical order with @code{sort()}:
1769 // l1 and l2 are now equal
1773 Finally, you can remove all but the first element of consecutive groups of
1774 elements with @code{unique()}:
1779 l3 = x, 2, 2, 2, y, x+y, y+x;
1780 l3.unique(); // l3 is now @{x, 2, y, x+y@}
1785 @node Mathematical functions, Relations, Lists, Basic concepts
1786 @c node-name, next, previous, up
1787 @section Mathematical functions
1788 @cindex @code{function} (class)
1789 @cindex trigonometric function
1790 @cindex hyperbolic function
1792 There are quite a number of useful functions hard-wired into GiNaC. For
1793 instance, all trigonometric and hyperbolic functions are implemented
1794 (@xref{Built-in functions}, for a complete list).
1796 These functions (better called @emph{pseudofunctions}) are all objects
1797 of class @code{function}. They accept one or more expressions as
1798 arguments and return one expression. If the arguments are not
1799 numerical, the evaluation of the function may be halted, as it does in
1800 the next example, showing how a function returns itself twice and
1801 finally an expression that may be really useful:
1803 @cindex Gamma function
1804 @cindex @code{subs()}
1807 symbol x("x"), y("y");
1809 cout << tgamma(foo) << endl;
1810 // -> tgamma(x+(1/2)*y)
1811 ex bar = foo.subs(y==1);
1812 cout << tgamma(bar) << endl;
1814 ex foobar = bar.subs(x==7);
1815 cout << tgamma(foobar) << endl;
1816 // -> (135135/128)*Pi^(1/2)
1820 Besides evaluation most of these functions allow differentiation, series
1821 expansion and so on. Read the next chapter in order to learn more about
1824 It must be noted that these pseudofunctions are created by inline
1825 functions, where the argument list is templated. This means that
1826 whenever you call @code{GiNaC::sin(1)} it is equivalent to
1827 @code{sin(ex(1))} and will therefore not result in a floating point
1828 number. Unless of course the function prototype is explicitly
1829 overridden -- which is the case for arguments of type @code{numeric}
1830 (not wrapped inside an @code{ex}). Hence, in order to obtain a floating
1831 point number of class @code{numeric} you should call
1832 @code{sin(numeric(1))}. This is almost the same as calling
1833 @code{sin(1).evalf()} except that the latter will return a numeric
1834 wrapped inside an @code{ex}.
1837 @node Relations, Integrals, Mathematical functions, Basic concepts
1838 @c node-name, next, previous, up
1840 @cindex @code{relational} (class)
1842 Sometimes, a relation holding between two expressions must be stored
1843 somehow. The class @code{relational} is a convenient container for such
1844 purposes. A relation is by definition a container for two @code{ex} and
1845 a relation between them that signals equality, inequality and so on.
1846 They are created by simply using the C++ operators @code{==}, @code{!=},
1847 @code{<}, @code{<=}, @code{>} and @code{>=} between two expressions.
1849 @xref{Mathematical functions}, for examples where various applications
1850 of the @code{.subs()} method show how objects of class relational are
1851 used as arguments. There they provide an intuitive syntax for
1852 substitutions. They are also used as arguments to the @code{ex::series}
1853 method, where the left hand side of the relation specifies the variable
1854 to expand in and the right hand side the expansion point. They can also
1855 be used for creating systems of equations that are to be solved for
1856 unknown variables. But the most common usage of objects of this class
1857 is rather inconspicuous in statements of the form @code{if
1858 (expand(pow(a+b,2))==a*a+2*a*b+b*b) @{...@}}. Here, an implicit
1859 conversion from @code{relational} to @code{bool} takes place. Note,
1860 however, that @code{==} here does not perform any simplifications, hence
1861 @code{expand()} must be called explicitly.
1863 @node Integrals, Matrices, Relations, Basic concepts
1864 @c node-name, next, previous, up
1866 @cindex @code{integral} (class)
1868 An object of class @dfn{integral} can be used to hold a symbolic integral.
1869 If you want to symbolically represent the integral of @code{x*x} from 0 to
1870 1, you would write this as
1872 integral(x, 0, 1, x*x)
1874 The first argument is the integration variable. It should be noted that
1875 GiNaC is not very good (yet?) at symbolically evaluating integrals. In
1876 fact, it can only integrate polynomials. An expression containing integrals
1877 can be evaluated symbolically by calling the
1881 method on it. Numerical evaluation is available by calling the
1885 method on an expression containing the integral. This will only evaluate
1886 integrals into a number if @code{subs}ing the integration variable by a
1887 number in the fourth argument of an integral and then @code{evalf}ing the
1888 result always results in a number. Of course, also the boundaries of the
1889 integration domain must @code{evalf} into numbers. It should be noted that
1890 trying to @code{evalf} a function with discontinuities in the integration
1891 domain is not recommended. The accuracy of the numeric evaluation of
1892 integrals is determined by the static member variable
1894 ex integral::relative_integration_error
1896 of the class @code{integral}. The default value of this is 10^-8.
1897 The integration works by halving the interval of integration, until numeric
1898 stability of the answer indicates that the requested accuracy has been
1899 reached. The maximum depth of the halving can be set via the static member
1902 int integral::max_integration_level
1904 The default value is 15. If this depth is exceeded, @code{evalf} will simply
1905 return the integral unevaluated. The function that performs the numerical
1906 evaluation, is also available as
1908 ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f,
1911 This function will throw an exception if the maximum depth is exceeded. The
1912 last parameter of the function is optional and defaults to the
1913 @code{relative_integration_error}. To make sure that we do not do too
1914 much work if an expression contains the same integral multiple times,
1915 a lookup table is used.
1917 If you know that an expression holds an integral, you can get the
1918 integration variable, the left boundary, right boundary and integrand by
1919 respectively calling @code{.op(0)}, @code{.op(1)}, @code{.op(2)}, and
1920 @code{.op(3)}. Differentiating integrals with respect to variables works
1921 as expected. Note that it makes no sense to differentiate an integral
1922 with respect to the integration variable.
1924 @node Matrices, Indexed objects, Integrals, Basic concepts
1925 @c node-name, next, previous, up
1927 @cindex @code{matrix} (class)
1929 A @dfn{matrix} is a two-dimensional array of expressions. The elements of a
1930 matrix with @math{m} rows and @math{n} columns are accessed with two
1931 @code{unsigned} indices, the first one in the range 0@dots{}@math{m-1}, the
1932 second one in the range 0@dots{}@math{n-1}.
1934 There are a couple of ways to construct matrices, with or without preset
1935 elements. The constructor
1938 matrix::matrix(unsigned r, unsigned c);
1941 creates a matrix with @samp{r} rows and @samp{c} columns with all elements
1944 The fastest way to create a matrix with preinitialized elements is to assign
1945 a list of comma-separated expressions to an empty matrix (see below for an
1946 example). But you can also specify the elements as a (flat) list with
1949 matrix::matrix(unsigned r, unsigned c, const lst & l);
1954 @cindex @code{lst_to_matrix()}
1956 ex lst_to_matrix(const lst & l);
1959 constructs a matrix from a list of lists, each list representing a matrix row.
1961 There is also a set of functions for creating some special types of
1964 @cindex @code{diag_matrix()}
1965 @cindex @code{unit_matrix()}
1966 @cindex @code{symbolic_matrix()}
1968 ex diag_matrix(const lst & l);
1969 ex unit_matrix(unsigned x);
1970 ex unit_matrix(unsigned r, unsigned c);
1971 ex symbolic_matrix(unsigned r, unsigned c, const string & base_name);
1972 ex symbolic_matrix(unsigned r, unsigned c, const string & base_name,
1973 const string & tex_base_name);
1976 @code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
1977 elements. @code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r}
1978 by @samp{c}) unit matrix. And finally, @code{symbolic_matrix} constructs a
1979 matrix filled with newly generated symbols made of the specified base name
1980 and the position of each element in the matrix.
1982 Matrices often arise by omitting elements of another matrix. For
1983 instance, the submatrix @code{S} of a matrix @code{M} takes a
1984 rectangular block from @code{M}. The reduced matrix @code{R} is defined
1985 by removing one row and one column from a matrix @code{M}. (The
1986 determinant of a reduced matrix is called a @emph{Minor} of @code{M} and
1987 can be used for computing the inverse using Cramer's rule.)
1989 @cindex @code{sub_matrix()}
1990 @cindex @code{reduced_matrix()}
1992 ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc);
1993 ex reduced_matrix(const matrix& m, unsigned r, unsigned c);
1996 The function @code{sub_matrix()} takes a row offset @code{r} and a
1997 column offset @code{c} and takes a block of @code{nr} rows and @code{nc}
1998 columns. The function @code{reduced_matrix()} has two integer arguments
1999 that specify which row and column to remove:
2007 cout << reduced_matrix(m, 1, 1) << endl;
2008 // -> [[11,13],[31,33]]
2009 cout << sub_matrix(m, 1, 2, 1, 2) << endl;
2010 // -> [[22,23],[32,33]]
2014 Matrix elements can be accessed and set using the parenthesis (function call)
2018 const ex & matrix::operator()(unsigned r, unsigned c) const;
2019 ex & matrix::operator()(unsigned r, unsigned c);
2022 It is also possible to access the matrix elements in a linear fashion with
2023 the @code{op()} method. But C++-style subscripting with square brackets
2024 @samp{[]} is not available.
2026 Here are a couple of examples for constructing matrices:
2030 symbol a("a"), b("b");
2044 cout << matrix(2, 2, lst(a, 0, 0, b)) << endl;
2047 cout << lst_to_matrix(lst(lst(a, 0), lst(0, b))) << endl;
2050 cout << diag_matrix(lst(a, b)) << endl;
2053 cout << unit_matrix(3) << endl;
2054 // -> [[1,0,0],[0,1,0],[0,0,1]]
2056 cout << symbolic_matrix(2, 3, "x") << endl;
2057 // -> [[x00,x01,x02],[x10,x11,x12]]
2061 @cindex @code{is_zero_matrix()}
2062 The method @code{matrix::is_zero_matrix()} returns @code{true} only if
2063 all entries of the matrix are zeros. There is also method
2064 @code{ex::is_zero_matrix()} which returns @code{true} only if the
2065 expression is zero or a zero matrix.
2067 @cindex @code{transpose()}
2068 There are three ways to do arithmetic with matrices. The first (and most
2069 direct one) is to use the methods provided by the @code{matrix} class:
2072 matrix matrix::add(const matrix & other) const;
2073 matrix matrix::sub(const matrix & other) const;
2074 matrix matrix::mul(const matrix & other) const;
2075 matrix matrix::mul_scalar(const ex & other) const;
2076 matrix matrix::pow(const ex & expn) const;
2077 matrix matrix::transpose() const;
2080 All of these methods return the result as a new matrix object. Here is an
2081 example that calculates @math{A*B-2*C} for three matrices @math{A}, @math{B}
2086 matrix A(2, 2), B(2, 2), C(2, 2);
2094 matrix result = A.mul(B).sub(C.mul_scalar(2));
2095 cout << result << endl;
2096 // -> [[-13,-6],[1,2]]
2101 @cindex @code{evalm()}
2102 The second (and probably the most natural) way is to construct an expression
2103 containing matrices with the usual arithmetic operators and @code{pow()}.
2104 For efficiency reasons, expressions with sums, products and powers of
2105 matrices are not automatically evaluated in GiNaC. You have to call the
2109 ex ex::evalm() const;
2112 to obtain the result:
2119 // -> [[1,2],[3,4]]*[[-1,0],[2,1]]-2*[[8,4],[2,1]]
2120 cout << e.evalm() << endl;
2121 // -> [[-13,-6],[1,2]]
2126 The non-commutativity of the product @code{A*B} in this example is
2127 automatically recognized by GiNaC. There is no need to use a special
2128 operator here. @xref{Non-commutative objects}, for more information about
2129 dealing with non-commutative expressions.
2131 Finally, you can work with indexed matrices and call @code{simplify_indexed()}
2132 to perform the arithmetic:
2137 idx i(symbol("i"), 2), j(symbol("j"), 2), k(symbol("k"), 2);
2138 e = indexed(A, i, k) * indexed(B, k, j) - 2 * indexed(C, i, j);
2140 // -> -2*[[8,4],[2,1]].i.j+[[-1,0],[2,1]].k.j*[[1,2],[3,4]].i.k
2141 cout << e.simplify_indexed() << endl;
2142 // -> [[-13,-6],[1,2]].i.j
2146 Using indices is most useful when working with rectangular matrices and
2147 one-dimensional vectors because you don't have to worry about having to
2148 transpose matrices before multiplying them. @xref{Indexed objects}, for
2149 more information about using matrices with indices, and about indices in
2152 The @code{matrix} class provides a couple of additional methods for
2153 computing determinants, traces, characteristic polynomials and ranks:
2155 @cindex @code{determinant()}
2156 @cindex @code{trace()}
2157 @cindex @code{charpoly()}
2158 @cindex @code{rank()}
2160 ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
2161 ex matrix::trace() const;
2162 ex matrix::charpoly(const ex & lambda) const;
2163 unsigned matrix::rank() const;
2166 The @samp{algo} argument of @code{determinant()} allows to select
2167 between different algorithms for calculating the determinant. The
2168 asymptotic speed (as parametrized by the matrix size) can greatly differ
2169 between those algorithms, depending on the nature of the matrix'
2170 entries. The possible values are defined in the @file{flags.h} header
2171 file. By default, GiNaC uses a heuristic to automatically select an
2172 algorithm that is likely (but not guaranteed) to give the result most
2175 @cindex @code{inverse()} (matrix)
2176 @cindex @code{solve()}
2177 Matrices may also be inverted using the @code{ex matrix::inverse()}
2178 method and linear systems may be solved with:
2181 matrix matrix::solve(const matrix & vars, const matrix & rhs,
2182 unsigned algo=solve_algo::automatic) const;
2185 Assuming the matrix object this method is applied on is an @code{m}
2186 times @code{n} matrix, then @code{vars} must be a @code{n} times
2187 @code{p} matrix of symbolic indeterminates and @code{rhs} a @code{m}
2188 times @code{p} matrix. The returned matrix then has dimension @code{n}
2189 times @code{p} and in the case of an underdetermined system will still
2190 contain some of the indeterminates from @code{vars}. If the system is
2191 overdetermined, an exception is thrown.
2194 @node Indexed objects, Non-commutative objects, Matrices, Basic concepts
2195 @c node-name, next, previous, up
2196 @section Indexed objects
2198 GiNaC allows you to handle expressions containing general indexed objects in
2199 arbitrary spaces. It is also able to canonicalize and simplify such
2200 expressions and perform symbolic dummy index summations. There are a number
2201 of predefined indexed objects provided, like delta and metric tensors.
2203 There are few restrictions placed on indexed objects and their indices and
2204 it is easy to construct nonsense expressions, but our intention is to
2205 provide a general framework that allows you to implement algorithms with
2206 indexed quantities, getting in the way as little as possible.
2208 @cindex @code{idx} (class)
2209 @cindex @code{indexed} (class)
2210 @subsection Indexed quantities and their indices
2212 Indexed expressions in GiNaC are constructed of two special types of objects,
2213 @dfn{index objects} and @dfn{indexed objects}.
2217 @cindex contravariant
2220 @item Index objects are of class @code{idx} or a subclass. Every index has
2221 a @dfn{value} and a @dfn{dimension} (which is the dimension of the space
2222 the index lives in) which can both be arbitrary expressions but are usually
2223 a number or a simple symbol. In addition, indices of class @code{varidx} have
2224 a @dfn{variance} (they can be co- or contravariant), and indices of class
2225 @code{spinidx} have a variance and can be @dfn{dotted} or @dfn{undotted}.
2227 @item Indexed objects are of class @code{indexed} or a subclass. They
2228 contain a @dfn{base expression} (which is the expression being indexed), and
2229 one or more indices.
2233 @strong{Please notice:} when printing expressions, covariant indices and indices
2234 without variance are denoted @samp{.i} while contravariant indices are
2235 denoted @samp{~i}. Dotted indices have a @samp{*} in front of the index
2236 value. In the following, we are going to use that notation in the text so
2237 instead of @math{A^i_jk} we will write @samp{A~i.j.k}. Index dimensions are
2238 not visible in the output.
2240 A simple example shall illustrate the concepts:
2244 #include <ginac/ginac.h>
2245 using namespace std;
2246 using namespace GiNaC;
2250 symbol i_sym("i"), j_sym("j");
2251 idx i(i_sym, 3), j(j_sym, 3);
2254 cout << indexed(A, i, j) << endl;
2256 cout << index_dimensions << indexed(A, i, j) << endl;
2258 cout << dflt; // reset cout to default output format (dimensions hidden)
2262 The @code{idx} constructor takes two arguments, the index value and the
2263 index dimension. First we define two index objects, @code{i} and @code{j},
2264 both with the numeric dimension 3. The value of the index @code{i} is the
2265 symbol @code{i_sym} (which prints as @samp{i}) and the value of the index
2266 @code{j} is the symbol @code{j_sym} (which prints as @samp{j}). Next we
2267 construct an expression containing one indexed object, @samp{A.i.j}. It has
2268 the symbol @code{A} as its base expression and the two indices @code{i} and
2271 The dimensions of indices are normally not visible in the output, but one
2272 can request them to be printed with the @code{index_dimensions} manipulator,
2275 Note the difference between the indices @code{i} and @code{j} which are of
2276 class @code{idx}, and the index values which are the symbols @code{i_sym}
2277 and @code{j_sym}. The indices of indexed objects cannot directly be symbols
2278 or numbers but must be index objects. For example, the following is not
2279 correct and will raise an exception:
2282 symbol i("i"), j("j");
2283 e = indexed(A, i, j); // ERROR: indices must be of type idx
2286 You can have multiple indexed objects in an expression, index values can
2287 be numeric, and index dimensions symbolic:
2291 symbol B("B"), dim("dim");
2292 cout << 4 * indexed(A, i)
2293 + indexed(B, idx(j_sym, 4), idx(2, 3), idx(i_sym, dim)) << endl;
2298 @code{B} has a 4-dimensional symbolic index @samp{k}, a 3-dimensional numeric
2299 index of value 2, and a symbolic index @samp{i} with the symbolic dimension
2300 @samp{dim}. Note that GiNaC doesn't automatically notify you that the free
2301 indices of @samp{A} and @samp{B} in the sum don't match (you have to call
2302 @code{simplify_indexed()} for that, see below).
2304 In fact, base expressions, index values and index dimensions can be
2305 arbitrary expressions:
2309 cout << indexed(A+B, idx(2*i_sym+1, dim/2)) << endl;
2314 It's also possible to construct nonsense like @samp{Pi.sin(x)}. You will not
2315 get an error message from this but you will probably not be able to do
2316 anything useful with it.
2318 @cindex @code{get_value()}
2319 @cindex @code{get_dimension()}
2323 ex idx::get_value();
2324 ex idx::get_dimension();
2327 return the value and dimension of an @code{idx} object. If you have an index
2328 in an expression, such as returned by calling @code{.op()} on an indexed
2329 object, you can get a reference to the @code{idx} object with the function
2330 @code{ex_to<idx>()} on the expression.
2332 There are also the methods
2335 bool idx::is_numeric();
2336 bool idx::is_symbolic();
2337 bool idx::is_dim_numeric();
2338 bool idx::is_dim_symbolic();
2341 for checking whether the value and dimension are numeric or symbolic
2342 (non-numeric). Using the @code{info()} method of an index (see @ref{Information
2343 about expressions}) returns information about the index value.
2345 @cindex @code{varidx} (class)
2346 If you need co- and contravariant indices, use the @code{varidx} class:
2350 symbol mu_sym("mu"), nu_sym("nu");
2351 varidx mu(mu_sym, 4), nu(nu_sym, 4); // default is contravariant ~mu, ~nu
2352 varidx mu_co(mu_sym, 4, true); // covariant index .mu
2354 cout << indexed(A, mu, nu) << endl;
2356 cout << indexed(A, mu_co, nu) << endl;
2358 cout << indexed(A, mu.toggle_variance(), nu) << endl;
2363 A @code{varidx} is an @code{idx} with an additional flag that marks it as
2364 co- or contravariant. The default is a contravariant (upper) index, but
2365 this can be overridden by supplying a third argument to the @code{varidx}
2366 constructor. The two methods
2369 bool varidx::is_covariant();
2370 bool varidx::is_contravariant();
2373 allow you to check the variance of a @code{varidx} object (use @code{ex_to<varidx>()}
2374 to get the object reference from an expression). There's also the very useful
2378 ex varidx::toggle_variance();
2381 which makes a new index with the same value and dimension but the opposite
2382 variance. By using it you only have to define the index once.
2384 @cindex @code{spinidx} (class)
2385 The @code{spinidx} class provides dotted and undotted variant indices, as
2386 used in the Weyl-van-der-Waerden spinor formalism:
2390 symbol K("K"), C_sym("C"), D_sym("D");
2391 spinidx C(C_sym, 2), D(D_sym); // default is 2-dimensional,
2392 // contravariant, undotted
2393 spinidx C_co(C_sym, 2, true); // covariant index
2394 spinidx D_dot(D_sym, 2, false, true); // contravariant, dotted
2395 spinidx D_co_dot(D_sym, 2, true, true); // covariant, dotted
2397 cout << indexed(K, C, D) << endl;
2399 cout << indexed(K, C_co, D_dot) << endl;
2401 cout << indexed(K, D_co_dot, D) << endl;
2406 A @code{spinidx} is a @code{varidx} with an additional flag that marks it as
2407 dotted or undotted. The default is undotted but this can be overridden by
2408 supplying a fourth argument to the @code{spinidx} constructor. The two
2412 bool spinidx::is_dotted();
2413 bool spinidx::is_undotted();
2416 allow you to check whether or not a @code{spinidx} object is dotted (use
2417 @code{ex_to<spinidx>()} to get the object reference from an expression).
2418 Finally, the two methods
2421 ex spinidx::toggle_dot();
2422 ex spinidx::toggle_variance_dot();
2425 create a new index with the same value and dimension but opposite dottedness
2426 and the same or opposite variance.
2428 @subsection Substituting indices
2430 @cindex @code{subs()}
2431 Sometimes you will want to substitute one symbolic index with another
2432 symbolic or numeric index, for example when calculating one specific element
2433 of a tensor expression. This is done with the @code{.subs()} method, as it
2434 is done for symbols (see @ref{Substituting expressions}).
2436 You have two possibilities here. You can either substitute the whole index
2437 by another index or expression:
2441 ex e = indexed(A, mu_co);
2442 cout << e << " becomes " << e.subs(mu_co == nu) << endl;
2443 // -> A.mu becomes A~nu
2444 cout << e << " becomes " << e.subs(mu_co == varidx(0, 4)) << endl;
2445 // -> A.mu becomes A~0
2446 cout << e << " becomes " << e.subs(mu_co == 0) << endl;
2447 // -> A.mu becomes A.0
2451 The third example shows that trying to replace an index with something that
2452 is not an index will substitute the index value instead.
2454 Alternatively, you can substitute the @emph{symbol} of a symbolic index by
2459 ex e = indexed(A, mu_co);
2460 cout << e << " becomes " << e.subs(mu_sym == nu_sym) << endl;
2461 // -> A.mu becomes A.nu
2462 cout << e << " becomes " << e.subs(mu_sym == 0) << endl;
2463 // -> A.mu becomes A.0
2467 As you see, with the second method only the value of the index will get
2468 substituted. Its other properties, including its dimension, remain unchanged.
2469 If you want to change the dimension of an index you have to substitute the
2470 whole index by another one with the new dimension.
2472 Finally, substituting the base expression of an indexed object works as
2477 ex e = indexed(A, mu_co);
2478 cout << e << " becomes " << e.subs(A == A+B) << endl;
2479 // -> A.mu becomes (B+A).mu
2483 @subsection Symmetries
2484 @cindex @code{symmetry} (class)
2485 @cindex @code{sy_none()}
2486 @cindex @code{sy_symm()}
2487 @cindex @code{sy_anti()}
2488 @cindex @code{sy_cycl()}
2490 Indexed objects can have certain symmetry properties with respect to their
2491 indices. Symmetries are specified as a tree of objects of class @code{symmetry}
2492 that is constructed with the helper functions
2495 symmetry sy_none(...);
2496 symmetry sy_symm(...);
2497 symmetry sy_anti(...);
2498 symmetry sy_cycl(...);
2501 @code{sy_none()} stands for no symmetry, @code{sy_symm()} and @code{sy_anti()}
2502 specify fully symmetric or antisymmetric, respectively, and @code{sy_cycl()}
2503 represents a cyclic symmetry. Each of these functions accepts up to four
2504 arguments which can be either symmetry objects themselves or unsigned integer
2505 numbers that represent an index position (counting from 0). A symmetry
2506 specification that consists of only a single @code{sy_symm()}, @code{sy_anti()}
2507 or @code{sy_cycl()} with no arguments specifies the respective symmetry for
2510 Here are some examples of symmetry definitions:
2515 e = indexed(A, i, j);
2516 e = indexed(A, sy_none(), i, j); // equivalent
2517 e = indexed(A, sy_none(0, 1), i, j); // equivalent
2519 // Symmetric in all three indices:
2520 e = indexed(A, sy_symm(), i, j, k);
2521 e = indexed(A, sy_symm(0, 1, 2), i, j, k); // equivalent
2522 e = indexed(A, sy_symm(2, 0, 1), i, j, k); // same symmetry, but yields a
2523 // different canonical order
2525 // Symmetric in the first two indices only:
2526 e = indexed(A, sy_symm(0, 1), i, j, k);
2527 e = indexed(A, sy_none(sy_symm(0, 1), 2), i, j, k); // equivalent
2529 // Antisymmetric in the first and last index only (index ranges need not
2531 e = indexed(A, sy_anti(0, 2), i, j, k);
2532 e = indexed(A, sy_none(sy_anti(0, 2), 1), i, j, k); // equivalent
2534 // An example of a mixed symmetry: antisymmetric in the first two and
2535 // last two indices, symmetric when swapping the first and last index
2536 // pairs (like the Riemann curvature tensor):
2537 e = indexed(A, sy_symm(sy_anti(0, 1), sy_anti(2, 3)), i, j, k, l);
2539 // Cyclic symmetry in all three indices:
2540 e = indexed(A, sy_cycl(), i, j, k);
2541 e = indexed(A, sy_cycl(0, 1, 2), i, j, k); // equivalent
2543 // The following examples are invalid constructions that will throw
2544 // an exception at run time.
2546 // An index may not appear multiple times:
2547 e = indexed(A, sy_symm(0, 0, 1), i, j, k); // ERROR
2548 e = indexed(A, sy_none(sy_symm(0, 1), sy_anti(0, 2)), i, j, k); // ERROR
2550 // Every child of sy_symm(), sy_anti() and sy_cycl() must refer to the
2551 // same number of indices:
2552 e = indexed(A, sy_symm(sy_anti(0, 1), 2), i, j, k); // ERROR
2554 // And of course, you cannot specify indices which are not there:
2555 e = indexed(A, sy_symm(0, 1, 2, 3), i, j, k); // ERROR
2559 If you need to specify more than four indices, you have to use the
2560 @code{.add()} method of the @code{symmetry} class. For example, to specify
2561 full symmetry in the first six indices you would write
2562 @code{sy_symm(0, 1, 2, 3).add(4).add(5)}.
2564 If an indexed object has a symmetry, GiNaC will automatically bring the
2565 indices into a canonical order which allows for some immediate simplifications:
2569 cout << indexed(A, sy_symm(), i, j)
2570 + indexed(A, sy_symm(), j, i) << endl;
2572 cout << indexed(B, sy_anti(), i, j)
2573 + indexed(B, sy_anti(), j, i) << endl;
2575 cout << indexed(B, sy_anti(), i, j, k)
2576 - indexed(B, sy_anti(), j, k, i) << endl;
2581 @cindex @code{get_free_indices()}
2583 @subsection Dummy indices
2585 GiNaC treats certain symbolic index pairs as @dfn{dummy indices} meaning
2586 that a summation over the index range is implied. Symbolic indices which are
2587 not dummy indices are called @dfn{free indices}. Numeric indices are neither
2588 dummy nor free indices.
2590 To be recognized as a dummy index pair, the two indices must be of the same
2591 class and their value must be the same single symbol (an index like
2592 @samp{2*n+1} is never a dummy index). If the indices are of class
2593 @code{varidx} they must also be of opposite variance; if they are of class
2594 @code{spinidx} they must be both dotted or both undotted.
2596 The method @code{.get_free_indices()} returns a vector containing the free
2597 indices of an expression. It also checks that the free indices of the terms
2598 of a sum are consistent:
2602 symbol A("A"), B("B"), C("C");
2604 symbol i_sym("i"), j_sym("j"), k_sym("k"), l_sym("l");
2605 idx i(i_sym, 3), j(j_sym, 3), k(k_sym, 3), l(l_sym, 3);
2607 ex e = indexed(A, i, j) * indexed(B, j, k) + indexed(C, k, l, i, l);
2608 cout << exprseq(e.get_free_indices()) << endl;
2610 // 'j' and 'l' are dummy indices
2612 symbol mu_sym("mu"), nu_sym("nu"), rho_sym("rho"), sigma_sym("sigma");
2613 varidx mu(mu_sym, 4), nu(nu_sym, 4), rho(rho_sym, 4), sigma(sigma_sym, 4);
2615 e = indexed(A, mu, nu) * indexed(B, nu.toggle_variance(), rho)
2616 + indexed(C, mu, sigma, rho, sigma.toggle_variance());
2617 cout << exprseq(e.get_free_indices()) << endl;
2619 // 'nu' is a dummy index, but 'sigma' is not
2621 e = indexed(A, mu, mu);
2622 cout << exprseq(e.get_free_indices()) << endl;
2624 // 'mu' is not a dummy index because it appears twice with the same
2627 e = indexed(A, mu, nu) + 42;
2628 cout << exprseq(e.get_free_indices()) << endl; // ERROR
2629 // this will throw an exception:
2630 // "add::get_free_indices: inconsistent indices in sum"
2634 @cindex @code{expand_dummy_sum()}
2635 A dummy index summation like
2642 can be expanded for indices with numeric
2643 dimensions (e.g. 3) into the explicit sum like
2645 $a_1b^1+a_2b^2+a_3b^3 $.
2648 a.1 b~1 + a.2 b~2 + a.3 b~3.
2650 This is performed by the function
2653 ex expand_dummy_sum(const ex & e, bool subs_idx = false);
2656 which takes an expression @code{e} and returns the expanded sum for all
2657 dummy indices with numeric dimensions. If the parameter @code{subs_idx}
2658 is set to @code{true} then all substitutions are made by @code{idx} class
2659 indices, i.e. without variance. In this case the above sum
2668 $a_1b_1+a_2b_2+a_3b_3 $.
2671 a.1 b.1 + a.2 b.2 + a.3 b.3.
2675 @cindex @code{simplify_indexed()}
2676 @subsection Simplifying indexed expressions
2678 In addition to the few automatic simplifications that GiNaC performs on
2679 indexed expressions (such as re-ordering the indices of symmetric tensors
2680 and calculating traces and convolutions of matrices and predefined tensors)
2684 ex ex::simplify_indexed();
2685 ex ex::simplify_indexed(const scalar_products & sp);
2688 that performs some more expensive operations:
2691 @item it checks the consistency of free indices in sums in the same way
2692 @code{get_free_indices()} does
2693 @item it tries to give dummy indices that appear in different terms of a sum
2694 the same name to allow simplifications like @math{a_i*b_i-a_j*b_j=0}
2695 @item it (symbolically) calculates all possible dummy index summations/contractions
2696 with the predefined tensors (this will be explained in more detail in the
2698 @item it detects contractions that vanish for symmetry reasons, for example
2699 the contraction of a symmetric and a totally antisymmetric tensor
2700 @item as a special case of dummy index summation, it can replace scalar products
2701 of two tensors with a user-defined value
2704 The last point is done with the help of the @code{scalar_products} class
2705 which is used to store scalar products with known values (this is not an
2706 arithmetic class, you just pass it to @code{simplify_indexed()}):
2710 symbol A("A"), B("B"), C("C"), i_sym("i");
2714 sp.add(A, B, 0); // A and B are orthogonal
2715 sp.add(A, C, 0); // A and C are orthogonal
2716 sp.add(A, A, 4); // A^2 = 4 (A has length 2)
2718 e = indexed(A + B, i) * indexed(A + C, i);
2720 // -> (B+A).i*(A+C).i
2722 cout << e.expand(expand_options::expand_indexed).simplify_indexed(sp)
2728 The @code{scalar_products} object @code{sp} acts as a storage for the
2729 scalar products added to it with the @code{.add()} method. This method
2730 takes three arguments: the two expressions of which the scalar product is
2731 taken, and the expression to replace it with.
2733 @cindex @code{expand()}
2734 The example above also illustrates a feature of the @code{expand()} method:
2735 if passed the @code{expand_indexed} option it will distribute indices
2736 over sums, so @samp{(A+B).i} becomes @samp{A.i+B.i}.
2738 @cindex @code{tensor} (class)
2739 @subsection Predefined tensors
2741 Some frequently used special tensors such as the delta, epsilon and metric
2742 tensors are predefined in GiNaC. They have special properties when
2743 contracted with other tensor expressions and some of them have constant
2744 matrix representations (they will evaluate to a number when numeric
2745 indices are specified).
2747 @cindex @code{delta_tensor()}
2748 @subsubsection Delta tensor
2750 The delta tensor takes two indices, is symmetric and has the matrix
2751 representation @code{diag(1, 1, 1, ...)}. It is constructed by the function
2752 @code{delta_tensor()}:
2756 symbol A("A"), B("B");
2758 idx i(symbol("i"), 3), j(symbol("j"), 3),
2759 k(symbol("k"), 3), l(symbol("l"), 3);
2761 ex e = indexed(A, i, j) * indexed(B, k, l)
2762 * delta_tensor(i, k) * delta_tensor(j, l);
2763 cout << e.simplify_indexed() << endl;
2766 cout << delta_tensor(i, i) << endl;
2771 @cindex @code{metric_tensor()}
2772 @subsubsection General metric tensor
2774 The function @code{metric_tensor()} creates a general symmetric metric
2775 tensor with two indices that can be used to raise/lower tensor indices. The
2776 metric tensor is denoted as @samp{g} in the output and if its indices are of
2777 mixed variance it is automatically replaced by a delta tensor:
2783 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
2785 ex e = metric_tensor(mu, nu) * indexed(A, nu.toggle_variance(), rho);
2786 cout << e.simplify_indexed() << endl;
2789 e = delta_tensor(mu, nu.toggle_variance()) * metric_tensor(nu, rho);
2790 cout << e.simplify_indexed() << endl;
2793 e = metric_tensor(mu.toggle_variance(), nu.toggle_variance())
2794 * metric_tensor(nu, rho);
2795 cout << e.simplify_indexed() << endl;
2798 e = metric_tensor(nu.toggle_variance(), rho.toggle_variance())
2799 * metric_tensor(mu, nu) * (delta_tensor(mu.toggle_variance(), rho)
2800 + indexed(A, mu.toggle_variance(), rho));
2801 cout << e.simplify_indexed() << endl;
2806 @cindex @code{lorentz_g()}
2807 @subsubsection Minkowski metric tensor
2809 The Minkowski metric tensor is a special metric tensor with a constant
2810 matrix representation which is either @code{diag(1, -1, -1, ...)} (negative
2811 signature, the default) or @code{diag(-1, 1, 1, ...)} (positive signature).
2812 It is created with the function @code{lorentz_g()} (although it is output as
2817 varidx mu(symbol("mu"), 4);
2819 e = delta_tensor(varidx(0, 4), mu.toggle_variance())
2820 * lorentz_g(mu, varidx(0, 4)); // negative signature
2821 cout << e.simplify_indexed() << endl;
2824 e = delta_tensor(varidx(0, 4), mu.toggle_variance())
2825 * lorentz_g(mu, varidx(0, 4), true); // positive signature
2826 cout << e.simplify_indexed() << endl;
2831 @cindex @code{spinor_metric()}
2832 @subsubsection Spinor metric tensor
2834 The function @code{spinor_metric()} creates an antisymmetric tensor with
2835 two indices that is used to raise/lower indices of 2-component spinors.
2836 It is output as @samp{eps}:
2842 spinidx A(symbol("A")), B(symbol("B")), C(symbol("C"));
2843 ex A_co = A.toggle_variance(), B_co = B.toggle_variance();
2845 e = spinor_metric(A, B) * indexed(psi, B_co);
2846 cout << e.simplify_indexed() << endl;
2849 e = spinor_metric(A, B) * indexed(psi, A_co);
2850 cout << e.simplify_indexed() << endl;
2853 e = spinor_metric(A_co, B_co) * indexed(psi, B);
2854 cout << e.simplify_indexed() << endl;
2857 e = spinor_metric(A_co, B_co) * indexed(psi, A);
2858 cout << e.simplify_indexed() << endl;
2861 e = spinor_metric(A_co, B_co) * spinor_metric(A, B);
2862 cout << e.simplify_indexed() << endl;
2865 e = spinor_metric(A_co, B_co) * spinor_metric(B, C);
2866 cout << e.simplify_indexed() << endl;
2871 The matrix representation of the spinor metric is @code{[[0, 1], [-1, 0]]}.
2873 @cindex @code{epsilon_tensor()}
2874 @cindex @code{lorentz_eps()}
2875 @subsubsection Epsilon tensor
2877 The epsilon tensor is totally antisymmetric, its number of indices is equal
2878 to the dimension of the index space (the indices must all be of the same
2879 numeric dimension), and @samp{eps.1.2.3...} (resp. @samp{eps~0~1~2...}) is
2880 defined to be 1. Its behavior with indices that have a variance also
2881 depends on the signature of the metric. Epsilon tensors are output as
2884 There are three functions defined to create epsilon tensors in 2, 3 and 4
2888 ex epsilon_tensor(const ex & i1, const ex & i2);
2889 ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
2890 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4,
2891 bool pos_sig = false);
2894 The first two functions create an epsilon tensor in 2 or 3 Euclidean
2895 dimensions, the last function creates an epsilon tensor in a 4-dimensional
2896 Minkowski space (the last @code{bool} argument specifies whether the metric
2897 has negative or positive signature, as in the case of the Minkowski metric
2902 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4),
2903 sig(symbol("sig"), 4), lam(symbol("lam"), 4), bet(symbol("bet"), 4);
2904 e = lorentz_eps(mu, nu, rho, sig) *
2905 lorentz_eps(mu.toggle_variance(), nu.toggle_variance(), lam, bet);
2906 cout << simplify_indexed(e) << endl;
2907 // -> 2*eta~bet~rho*eta~sig~lam-2*eta~sig~bet*eta~rho~lam
2909 idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
2910 symbol A("A"), B("B");
2911 e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(B, k);
2912 cout << simplify_indexed(e) << endl;
2913 // -> -B.k*A.j*eps.i.k.j
2914 e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(A, k);
2915 cout << simplify_indexed(e) << endl;
2920 @subsection Linear algebra
2922 The @code{matrix} class can be used with indices to do some simple linear
2923 algebra (linear combinations and products of vectors and matrices, traces
2924 and scalar products):
2928 idx i(symbol("i"), 2), j(symbol("j"), 2);
2929 symbol x("x"), y("y");
2931 // A is a 2x2 matrix, X is a 2x1 vector
2932 matrix A(2, 2), X(2, 1);
2937 cout << indexed(A, i, i) << endl;
2940 ex e = indexed(A, i, j) * indexed(X, j);
2941 cout << e.simplify_indexed() << endl;
2942 // -> [[2*y+x],[4*y+3*x]].i
2944 e = indexed(A, i, j) * indexed(X, i) + indexed(X, j) * 2;
2945 cout << e.simplify_indexed() << endl;
2946 // -> [[3*y+3*x,6*y+2*x]].j
2950 You can of course obtain the same results with the @code{matrix::add()},
2951 @code{matrix::mul()} and @code{matrix::trace()} methods (@pxref{Matrices})
2952 but with indices you don't have to worry about transposing matrices.
2954 Matrix indices always start at 0 and their dimension must match the number
2955 of rows/columns of the matrix. Matrices with one row or one column are
2956 vectors and can have one or two indices (it doesn't matter whether it's a
2957 row or a column vector). Other matrices must have two indices.
2959 You should be careful when using indices with variance on matrices. GiNaC
2960 doesn't look at the variance and doesn't know that @samp{F~mu~nu} and
2961 @samp{F.mu.nu} are different matrices. In this case you should use only
2962 one form for @samp{F} and explicitly multiply it with a matrix representation
2963 of the metric tensor.
2966 @node Non-commutative objects, Hash maps, Indexed objects, Basic concepts
2967 @c node-name, next, previous, up
2968 @section Non-commutative objects
2970 GiNaC is equipped to handle certain non-commutative algebras. Three classes of
2971 non-commutative objects are built-in which are mostly of use in high energy
2975 @item Clifford (Dirac) algebra (class @code{clifford})
2976 @item su(3) Lie algebra (class @code{color})
2977 @item Matrices (unindexed) (class @code{matrix})
2980 The @code{clifford} and @code{color} classes are subclasses of
2981 @code{indexed} because the elements of these algebras usually carry
2982 indices. The @code{matrix} class is described in more detail in
2985 Unlike most computer algebra systems, GiNaC does not primarily provide an
2986 operator (often denoted @samp{&*}) for representing inert products of
2987 arbitrary objects. Rather, non-commutativity in GiNaC is a property of the
2988 classes of objects involved, and non-commutative products are formed with
2989 the usual @samp{*} operator, as are ordinary products. GiNaC is capable of
2990 figuring out by itself which objects commutate and will group the factors
2991 by their class. Consider this example:
2995 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
2996 idx a(symbol("a"), 8), b(symbol("b"), 8);
2997 ex e = -dirac_gamma(mu) * (2*color_T(a)) * 8 * color_T(b) * dirac_gamma(nu);
2999 // -> -16*(gamma~mu*gamma~nu)*(T.a*T.b)
3003 As can be seen, GiNaC pulls out the overall commutative factor @samp{-16} and
3004 groups the non-commutative factors (the gammas and the su(3) generators)
3005 together while preserving the order of factors within each class (because
3006 Clifford objects commutate with color objects). The resulting expression is a
3007 @emph{commutative} product with two factors that are themselves non-commutative
3008 products (@samp{gamma~mu*gamma~nu} and @samp{T.a*T.b}). For clarification,
3009 parentheses are placed around the non-commutative products in the output.
3011 @cindex @code{ncmul} (class)
3012 Non-commutative products are internally represented by objects of the class
3013 @code{ncmul}, as opposed to commutative products which are handled by the
3014 @code{mul} class. You will normally not have to worry about this distinction,
3017 The advantage of this approach is that you never have to worry about using
3018 (or forgetting to use) a special operator when constructing non-commutative
3019 expressions. Also, non-commutative products in GiNaC are more intelligent
3020 than in other computer algebra systems; they can, for example, automatically
3021 canonicalize themselves according to rules specified in the implementation
3022 of the non-commutative classes. The drawback is that to work with other than
3023 the built-in algebras you have to implement new classes yourself. Both
3024 symbols and user-defined functions can be specified as being non-commutative.
3026 @cindex @code{return_type()}
3027 @cindex @code{return_type_tinfo()}
3028 Information about the commutativity of an object or expression can be
3029 obtained with the two member functions
3032 unsigned ex::return_type() const;
3033 unsigned ex::return_type_tinfo() const;
3036 The @code{return_type()} function returns one of three values (defined in
3037 the header file @file{flags.h}), corresponding to three categories of
3038 expressions in GiNaC:
3041 @item @code{return_types::commutative}: Commutates with everything. Most GiNaC
3042 classes are of this kind.
3043 @item @code{return_types::noncommutative}: Non-commutative, belonging to a
3044 certain class of non-commutative objects which can be determined with the
3045 @code{return_type_tinfo()} method. Expressions of this category commutate
3046 with everything except @code{noncommutative} expressions of the same
3048 @item @code{return_types::noncommutative_composite}: Non-commutative, composed
3049 of non-commutative objects of different classes. Expressions of this
3050 category don't commutate with any other @code{noncommutative} or
3051 @code{noncommutative_composite} expressions.
3054 The value returned by the @code{return_type_tinfo()} method is valid only
3055 when the return type of the expression is @code{noncommutative}. It is a
3056 value that is unique to the class of the object and usually one of the
3057 constants in @file{tinfos.h}, or derived therefrom.
3059 Here are a couple of examples:
3062 @multitable @columnfractions 0.33 0.33 0.34
3063 @item @strong{Expression} @tab @strong{@code{return_type()}} @tab @strong{@code{return_type_tinfo()}}
3064 @item @code{42} @tab @code{commutative} @tab -
3065 @item @code{2*x-y} @tab @code{commutative} @tab -
3066 @item @code{dirac_ONE()} @tab @code{noncommutative} @tab @code{TINFO_clifford}
3067 @item @code{dirac_gamma(mu)*dirac_gamma(nu)} @tab @code{noncommutative} @tab @code{TINFO_clifford}
3068 @item @code{2*color_T(a)} @tab @code{noncommutative} @tab @code{TINFO_color}
3069 @item @code{dirac_ONE()*color_T(a)} @tab @code{noncommutative_composite} @tab -
3073 Note: the @code{return_type_tinfo()} of Clifford objects is only equal to
3074 @code{TINFO_clifford} for objects with a representation label of zero.
3075 Other representation labels yield a different @code{return_type_tinfo()},
3076 but it's the same for any two objects with the same label. This is also true
3079 A last note: With the exception of matrices, positive integer powers of
3080 non-commutative objects are automatically expanded in GiNaC. For example,
3081 @code{pow(a*b, 2)} becomes @samp{a*b*a*b} if @samp{a} and @samp{b} are
3082 non-commutative expressions).
3085 @cindex @code{clifford} (class)
3086 @subsection Clifford algebra
3089 Clifford algebras are supported in two flavours: Dirac gamma
3090 matrices (more physical) and generic Clifford algebras (more
3093 @cindex @code{dirac_gamma()}
3094 @subsubsection Dirac gamma matrices
3095 Dirac gamma matrices (note that GiNaC doesn't treat them
3096 as matrices) are designated as @samp{gamma~mu} and satisfy
3097 @samp{gamma~mu*gamma~nu + gamma~nu*gamma~mu = 2*eta~mu~nu} where
3098 @samp{eta~mu~nu} is the Minkowski metric tensor. Dirac gammas are
3099 constructed by the function
3102 ex dirac_gamma(const ex & mu, unsigned char rl = 0);
3105 which takes two arguments: the index and a @dfn{representation label} in the
3106 range 0 to 255 which is used to distinguish elements of different Clifford
3107 algebras (this is also called a @dfn{spin line index}). Gammas with different
3108 labels commutate with each other. The dimension of the index can be 4 or (in
3109 the framework of dimensional regularization) any symbolic value. Spinor
3110 indices on Dirac gammas are not supported in GiNaC.
3112 @cindex @code{dirac_ONE()}
3113 The unity element of a Clifford algebra is constructed by
3116 ex dirac_ONE(unsigned char rl = 0);
3119 @strong{Please notice:} You must always use @code{dirac_ONE()} when referring to
3120 multiples of the unity element, even though it's customary to omit it.
3121 E.g. instead of @code{dirac_gamma(mu)*(dirac_slash(q,4)+m)} you have to
3122 write @code{dirac_gamma(mu)*(dirac_slash(q,4)+m*dirac_ONE())}. Otherwise,
3123 GiNaC will complain and/or produce incorrect results.
3125 @cindex @code{dirac_gamma5()}
3126 There is a special element @samp{gamma5} that commutates with all other
3127 gammas, has a unit square, and in 4 dimensions equals
3128 @samp{gamma~0 gamma~1 gamma~2 gamma~3}, provided by
3131 ex dirac_gamma5(unsigned char rl = 0);
3134 @cindex @code{dirac_gammaL()}
3135 @cindex @code{dirac_gammaR()}
3136 The chiral projectors @samp{(1+/-gamma5)/2} are also available as proper
3137 objects, constructed by
3140 ex dirac_gammaL(unsigned char rl = 0);
3141 ex dirac_gammaR(unsigned char rl = 0);
3144 They observe the relations @samp{gammaL^2 = gammaL}, @samp{gammaR^2 = gammaR},
3145 and @samp{gammaL gammaR = gammaR gammaL = 0}.
3147 @cindex @code{dirac_slash()}
3148 Finally, the function
3151 ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
3154 creates a term that represents a contraction of @samp{e} with the Dirac
3155 Lorentz vector (it behaves like a term of the form @samp{e.mu gamma~mu}
3156 with a unique index whose dimension is given by the @code{dim} argument).
3157 Such slashed expressions are printed with a trailing backslash, e.g. @samp{e\}.
3159 In products of dirac gammas, superfluous unity elements are automatically
3160 removed, squares are replaced by their values, and @samp{gamma5}, @samp{gammaL}
3161 and @samp{gammaR} are moved to the front.
3163 The @code{simplify_indexed()} function performs contractions in gamma strings,
3169 symbol a("a"), b("b"), D("D");
3170 varidx mu(symbol("mu"), D);
3171 ex e = dirac_gamma(mu) * dirac_slash(a, D)
3172 * dirac_gamma(mu.toggle_variance());
3174 // -> gamma~mu*a\*gamma.mu
3175 e = e.simplify_indexed();
3178 cout << e.subs(D == 4) << endl;
3184 @cindex @code{dirac_trace()}
3185 To calculate the trace of an expression containing strings of Dirac gammas
3186 you use one of the functions
3189 ex dirac_trace(const ex & e, const std::set<unsigned char> & rls,
3190 const ex & trONE = 4);
3191 ex dirac_trace(const ex & e, const lst & rll, const ex & trONE = 4);
3192 ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4);
3195 These functions take the trace over all gammas in the specified set @code{rls}
3196 or list @code{rll} of representation labels, or the single label @code{rl};
3197 gammas with other labels are left standing. The last argument to
3198 @code{dirac_trace()} is the value to be returned for the trace of the unity
3199 element, which defaults to 4.
3201 The @code{dirac_trace()} function is a linear functional that is equal to the
3202 ordinary matrix trace only in @math{D = 4} dimensions. In particular, the
3203 functional is not cyclic in
3206 dimensions when acting on
3207 expressions containing @samp{gamma5}, so it's not a proper trace. This
3208 @samp{gamma5} scheme is described in greater detail in
3209 @cite{The Role of gamma5 in Dimensional Regularization}.
3211 The value of the trace itself is also usually different in 4 and in
3219 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
3220 ex e = dirac_gamma(mu) * dirac_gamma(nu) *
3221 dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
3222 cout << dirac_trace(e).simplify_indexed() << endl;
3229 varidx mu(symbol("mu"), D), nu(symbol("nu"), D), rho(symbol("rho"), D);
3230 ex e = dirac_gamma(mu) * dirac_gamma(nu) *
3231 dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
3232 cout << dirac_trace(e).simplify_indexed() << endl;
3233 // -> 8*eta~rho~nu-4*eta~rho~nu*D
3237 Here is an example for using @code{dirac_trace()} to compute a value that
3238 appears in the calculation of the one-loop vacuum polarization amplitude in
3243 symbol q("q"), l("l"), m("m"), ldotq("ldotq"), D("D");
3244 varidx mu(symbol("mu"), D), nu(symbol("nu"), D);
3247 sp.add(l, l, pow(l, 2));
3248 sp.add(l, q, ldotq);
3250 ex e = dirac_gamma(mu) *
3251 (dirac_slash(l, D) + dirac_slash(q, D) + m * dirac_ONE()) *
3252 dirac_gamma(mu.toggle_variance()) *
3253 (dirac_slash(l, D) + m * dirac_ONE());
3254 e = dirac_trace(e).simplify_indexed(sp);
3255 e = e.collect(lst(l, ldotq, m));
3257 // -> (8-4*D)*l^2+(8-4*D)*ldotq+4*D*m^2
3261 The @code{canonicalize_clifford()} function reorders all gamma products that
3262 appear in an expression to a canonical (but not necessarily simple) form.
3263 You can use this to compare two expressions or for further simplifications:
3267 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
3268 ex e = dirac_gamma(mu) * dirac_gamma(nu) + dirac_gamma(nu) * dirac_gamma(mu);
3270 // -> gamma~mu*gamma~nu+gamma~nu*gamma~mu
3272 e = canonicalize_clifford(e);
3274 // -> 2*ONE*eta~mu~nu
3278 @cindex @code{clifford_unit()}
3279 @subsubsection A generic Clifford algebra
3281 A generic Clifford algebra, i.e. a
3285 dimensional algebra with
3289 satisfying the identities
3291 $e_i e_j + e_j e_i = M(i, j) + M(j, i) $
3294 e~i e~j + e~j e~i = M(i, j) + M(j, i)
3296 for some bilinear form (@code{metric})
3297 @math{M(i, j)}, which may be non-symmetric (see arXiv:math.QA/9911180)
3298 and contain symbolic entries. Such generators are created by the
3302 ex clifford_unit(const ex & mu, const ex & metr, unsigned char rl = 0);
3305 where @code{mu} should be a @code{idx} (or descendant) class object
3306 indexing the generators.
3307 Parameter @code{metr} defines the metric @math{M(i, j)} and can be
3308 represented by a square @code{matrix}, @code{tensormetric} or @code{indexed} class
3309 object. In fact, any expression either with two free indices or without
3310 indices at all is admitted as @code{metr}. In the later case an @code{indexed}
3311 object with two newly created indices with @code{metr} as its
3312 @code{op(0)} will be used.
3313 Optional parameter @code{rl} allows to distinguish different
3314 Clifford algebras, which will commute with each other.
3316 Note that the call @code{clifford_unit(mu, minkmetric())} creates
3317 something very close to @code{dirac_gamma(mu)}, although
3318 @code{dirac_gamma} have more efficient simplification mechanism.
3319 @cindex @code{clifford::get_metric()}
3320 The method @code{clifford::get_metric()} returns a metric defining this
3323 If the matrix @math{M(i, j)} is in fact symmetric you may prefer to create
3324 the Clifford algebra units with a call like that
3327 ex e = clifford_unit(mu, indexed(M, sy_symm(), i, j));
3330 since this may yield some further automatic simplifications. Again, for a
3331 metric defined through a @code{matrix} such a symmetry is detected
3334 Individual generators of a Clifford algebra can be accessed in several
3340 idx i(symbol("i"), 4);
3342 ex M = diag_matrix(lst(1, -1, 0, s));
3343 ex e = clifford_unit(i, M);
3344 ex e0 = e.subs(i == 0);
3345 ex e1 = e.subs(i == 1);
3346 ex e2 = e.subs(i == 2);
3347 ex e3 = e.subs(i == 3);
3352 will produce four anti-commuting generators of a Clifford algebra with properties
3354 $e_0^2=1 $, $e_1^2=-1$, $e_2^2=0$ and $e_3^2=s$.
3357 @code{pow(e0, 2) = 1}, @code{pow(e1, 2) = -1}, @code{pow(e2, 2) = 0} and
3358 @code{pow(e3, 2) = s}.
3361 @cindex @code{lst_to_clifford()}
3362 A similar effect can be achieved from the function
3365 ex lst_to_clifford(const ex & v, const ex & mu, const ex & metr,
3366 unsigned char rl = 0);
3367 ex lst_to_clifford(const ex & v, const ex & e);
3370 which converts a list or vector
3372 $v = (v^0, v^1, ..., v^n)$
3375 @samp{v = (v~0, v~1, ..., v~n)}
3380 $v^0 e_0 + v^1 e_1 + ... + v^n e_n$
3383 @samp{v~0 e.0 + v~1 e.1 + ... + v~n e.n}
3386 directly supplied in the second form of the procedure. In the first form
3387 the Clifford unit @samp{e.k} is generated by the call of
3388 @code{clifford_unit(mu, metr, rl)}. The previous code may be rewritten
3389 with the help of @code{lst_to_clifford()} as follows
3394 idx i(symbol("i"), 4);
3396 ex M = diag_matrix(lst(1, -1, 0, s));
3397 ex e0 = lst_to_clifford(lst(1, 0, 0, 0), i, M);
3398 ex e1 = lst_to_clifford(lst(0, 1, 0, 0), i, M);
3399 ex e2 = lst_to_clifford(lst(0, 0, 1, 0), i, M);
3400 ex e3 = lst_to_clifford(lst(0, 0, 0, 1), i, M);
3405 @cindex @code{clifford_to_lst()}
3406 There is the inverse function
3409 lst clifford_to_lst(const ex & e, const ex & c, bool algebraic = true);
3412 which takes an expression @code{e} and tries to find a list
3414 $v = (v^0, v^1, ..., v^n)$
3417 @samp{v = (v~0, v~1, ..., v~n)}
3421 $e = v^0 c_0 + v^1 c_1 + ... + v^n c_n$
3424 @samp{e = v~0 c.0 + v~1 c.1 + ... + v~n c.n}
3426 with respect to the given Clifford units @code{c} and with none of the
3427 @samp{v~k} containing Clifford units @code{c} (of course, this
3428 may be impossible). This function can use an @code{algebraic} method
3429 (default) or a symbolic one. With the @code{algebraic} method the @samp{v~k} are calculated as
3431 $(e c_k + c_k e)/c_k^2$. If $c_k^2$
3434 @samp{(e c.k + c.k e)/pow(c.k, 2)}. If @samp{pow(c.k, 2)}
3436 is zero or is not @code{numeric} for some @samp{k}
3437 then the method will be automatically changed to symbolic. The same effect
3438 is obtained by the assignment (@code{algebraic = false}) in the procedure call.
3440 @cindex @code{clifford_prime()}
3441 @cindex @code{clifford_star()}
3442 @cindex @code{clifford_bar()}
3443 There are several functions for (anti-)automorphisms of Clifford algebras:
3446 ex clifford_prime(const ex & e)
3447 inline ex clifford_star(const ex & e) @{ return e.conjugate(); @}
3448 inline ex clifford_bar(const ex & e) @{ return clifford_prime(e.conjugate()); @}
3451 The automorphism of a Clifford algebra @code{clifford_prime()} simply
3452 changes signs of all Clifford units in the expression. The reversion
3453 of a Clifford algebra @code{clifford_star()} coincides with the
3454 @code{conjugate()} method and effectively reverses the order of Clifford
3455 units in any product. Finally the main anti-automorphism
3456 of a Clifford algebra @code{clifford_bar()} is the composition of the
3457 previous two, i.e. it makes the reversion and changes signs of all Clifford units
3458 in a product. These functions correspond to the notations
3473 used in Clifford algebra textbooks.
3475 @cindex @code{clifford_norm()}
3479 ex clifford_norm(const ex & e);
3482 @cindex @code{clifford_inverse()}
3483 calculates the norm of a Clifford number from the expression
3485 $||e||^2 = e\overline{e}$.
3488 @code{||e||^2 = e \bar@{e@}}
3490 The inverse of a Clifford expression is returned by the function
3493 ex clifford_inverse(const ex & e);
3496 which calculates it as
3498 $e^{-1} = \overline{e}/||e||^2$.
3501 @math{e^@{-1@} = \bar@{e@}/||e||^2}
3510 then an exception is raised.
3512 @cindex @code{remove_dirac_ONE()}
3513 If a Clifford number happens to be a factor of
3514 @code{dirac_ONE()} then we can convert it to a ``real'' (non-Clifford)
3515 expression by the function
3518 ex remove_dirac_ONE(const ex & e);
3521 @cindex @code{canonicalize_clifford()}
3522 The function @code{canonicalize_clifford()} works for a
3523 generic Clifford algebra in a similar way as for Dirac gammas.
3525 The next provided function is
3527 @cindex @code{clifford_moebius_map()}
3529 ex clifford_moebius_map(const ex & a, const ex & b, const ex & c,
3530 const ex & d, const ex & v, const ex & G,
3531 unsigned char rl = 0);
3532 ex clifford_moebius_map(const ex & M, const ex & v, const ex & G,
3533 unsigned char rl = 0);
3536 It takes a list or vector @code{v} and makes the Moebius (conformal or
3537 linear-fractional) transformation @samp{v -> (av+b)/(cv+d)} defined by
3538 the matrix @samp{M = [[a, b], [c, d]]}. The parameter @code{G} defines
3539 the metric of the surrounding (pseudo-)Euclidean space. This can be an
3540 indexed object, tensormetric, matrix or a Clifford unit, in the later
3541 case the optional parameter @code{rl} is ignored even if supplied.
3542 Depending from the type of @code{v} the returned value of this function
3543 is either a vector or a list holding vector's components.
3545 @cindex @code{clifford_max_label()}
3546 Finally the function
3549 char clifford_max_label(const ex & e, bool ignore_ONE = false);
3552 can detect a presence of Clifford objects in the expression @code{e}: if
3553 such objects are found it returns the maximal
3554 @code{representation_label} of them, otherwise @code{-1}. The optional
3555 parameter @code{ignore_ONE} indicates if @code{dirac_ONE} objects should
3556 be ignored during the search.
3558 LaTeX output for Clifford units looks like
3559 @code{\clifford[1]@{e@}^@{@{\nu@}@}}, where @code{1} is the
3560 @code{representation_label} and @code{\nu} is the index of the
3561 corresponding unit. This provides a flexible typesetting with a suitable
3562 defintion of the @code{\clifford} command. For example, the definition
3564 \newcommand@{\clifford@}[1][]@{@}
3566 typesets all Clifford units identically, while the alternative definition
3568 \newcommand@{\clifford@}[2][]@{\ifcase #1 #2\or \tilde@{#2@} \or \breve@{#2@} \fi@}
3570 prints units with @code{representation_label=0} as
3577 with @code{representation_label=1} as
3584 and with @code{representation_label=2} as
3592 @cindex @code{color} (class)
3593 @subsection Color algebra
3595 @cindex @code{color_T()}
3596 For computations in quantum chromodynamics, GiNaC implements the base elements
3597 and structure constants of the su(3) Lie algebra (color algebra). The base
3598 elements @math{T_a} are constructed by the function
3601 ex color_T(const ex & a, unsigned char rl = 0);
3604 which takes two arguments: the index and a @dfn{representation label} in the
3605 range 0 to 255 which is used to distinguish elements of different color
3606 algebras. Objects with different labels commutate with each other. The
3607 dimension of the index must be exactly 8 and it should be of class @code{idx},
3610 @cindex @code{color_ONE()}
3611 The unity element of a color algebra is constructed by
3614 ex color_ONE(unsigned char rl = 0);
3617 @strong{Please notice:} You must always use @code{color_ONE()} when referring to
3618 multiples of the unity element, even though it's customary to omit it.
3619 E.g. instead of @code{color_T(a)*(color_T(b)*indexed(X,b)+1)} you have to
3620 write @code{color_T(a)*(color_T(b)*indexed(X,b)+color_ONE())}. Otherwise,
3621 GiNaC may produce incorrect results.
3623 @cindex @code{color_d()}
3624 @cindex @code{color_f()}
3628 ex color_d(const ex & a, const ex & b, const ex & c);
3629 ex color_f(const ex & a, const ex & b, const ex & c);
3632 create the symmetric and antisymmetric structure constants @math{d_abc} and
3633 @math{f_abc} which satisfy @math{@{T_a, T_b@} = 1/3 delta_ab + d_abc T_c}
3634 and @math{[T_a, T_b] = i f_abc T_c}.
3636 These functions evaluate to their numerical values,
3637 if you supply numeric indices to them. The index values should be in
3638 the range from 1 to 8, not from 0 to 7. This departure from usual conventions
3639 goes along better with the notations used in physical literature.
3641 @cindex @code{color_h()}
3642 There's an additional function
3645 ex color_h(const ex & a, const ex & b, const ex & c);
3648 which returns the linear combination @samp{color_d(a, b, c)+I*color_f(a, b, c)}.
3650 The function @code{simplify_indexed()} performs some simplifications on
3651 expressions containing color objects:
3656 idx a(symbol("a"), 8), b(symbol("b"), 8), c(symbol("c"), 8),
3657 k(symbol("k"), 8), l(symbol("l"), 8);
3659 e = color_d(a, b, l) * color_f(a, b, k);
3660 cout << e.simplify_indexed() << endl;
3663 e = color_d(a, b, l) * color_d(a, b, k);
3664 cout << e.simplify_indexed() << endl;
3667 e = color_f(l, a, b) * color_f(a, b, k);
3668 cout << e.simplify_indexed() << endl;
3671 e = color_h(a, b, c) * color_h(a, b, c);
3672 cout << e.simplify_indexed() << endl;
3675 e = color_h(a, b, c) * color_T(b) * color_T(c);
3676 cout << e.simplify_indexed() << endl;
3679 e = color_h(a, b, c) * color_T(a) * color_T(b) * color_T(c);
3680 cout << e.simplify_indexed() << endl;
3683 e = color_T(k) * color_T(a) * color_T(b) * color_T(k);
3684 cout << e.simplify_indexed() << endl;
3685 // -> 1/4*delta.b.a*ONE-1/6*T.a*T.b
3689 @cindex @code{color_trace()}
3690 To calculate the trace of an expression containing color objects you use one
3694 ex color_trace(const ex & e, const std::set<unsigned char> & rls);
3695 ex color_trace(const ex & e, const lst & rll);
3696 ex color_trace(const ex & e, unsigned char rl = 0);
3699 These functions take the trace over all color @samp{T} objects in the
3700 specified set @code{rls} or list @code{rll} of representation labels, or the
3701 single label @code{rl}; @samp{T}s with other labels are left standing. For
3706 e = color_trace(4 * color_T(a) * color_T(b) * color_T(c));
3708 // -> -I*f.a.c.b+d.a.c.b
3713 @node Hash maps, Methods and functions, Non-commutative objects, Basic concepts
3714 @c node-name, next, previous, up
3717 @cindex @code{exhashmap} (class)
3719 For your convenience, GiNaC offers the container template @code{exhashmap<T>}
3720 that can be used as a drop-in replacement for the STL
3721 @code{std::map<ex, T, ex_is_less>}, using hash tables to provide faster,
3722 typically constant-time, element look-up than @code{map<>}.
3724 @code{exhashmap<>} supports all @code{map<>} members and operations, with the
3725 following differences:
3729 no @code{lower_bound()} and @code{upper_bound()} methods
3731 no reverse iterators, no @code{rbegin()}/@code{rend()}
3733 no @code{operator<(exhashmap, exhashmap)}
3735 the comparison function object @code{key_compare} is hardcoded to
3738 the constructor @code{exhashmap(size_t n)} allows specifying the minimum
3739 initial hash table size (the actual table size after construction may be
3740 larger than the specified value)
3742 the method @code{size_t bucket_count()} returns the current size of the hash
3745 @code{insert()} and @code{erase()} operations invalidate all iterators
3749 @node Methods and functions, Information about expressions, Hash maps, Top
3750 @c node-name, next, previous, up
3751 @chapter Methods and functions
3754 In this chapter the most important algorithms provided by GiNaC will be
3755 described. Some of them are implemented as functions on expressions,
3756 others are implemented as methods provided by expression objects. If
3757 they are methods, there exists a wrapper function around it, so you can
3758 alternatively call it in a functional way as shown in the simple
3763 cout << "As method: " << sin(1).evalf() << endl;
3764 cout << "As function: " << evalf(sin(1)) << endl;
3768 @cindex @code{subs()}
3769 The general rule is that wherever methods accept one or more parameters
3770 (@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
3771 wrapper accepts is the same but preceded by the object to act on
3772 (@var{object}, @var{arg1}, @var{arg2}, @dots{}). This approach is the
3773 most natural one in an OO model but it may lead to confusion for MapleV
3774 users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
3775 would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
3776 @code{A} and @code{x}). On the other hand, since MapleV returns 3 on
3777 @code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
3778 coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
3779 here. Also, users of MuPAD will in most cases feel more comfortable
3780 with GiNaC's convention. All function wrappers are implemented
3781 as simple inline functions which just call the corresponding method and
3782 are only provided for users uncomfortable with OO who are dead set to
3783 avoid method invocations. Generally, nested function wrappers are much
3784 harder to read than a sequence of methods and should therefore be
3785 avoided if possible. On the other hand, not everything in GiNaC is a
3786 method on class @code{ex} and sometimes calling a function cannot be
3790 * Information about expressions::
3791 * Numerical evaluation::
3792 * Substituting expressions::
3793 * Pattern matching and advanced substitutions::
3794 * Applying a function on subexpressions::
3795 * Visitors and tree traversal::
3796 * Polynomial arithmetic:: Working with polynomials.
3797 * Rational expressions:: Working with rational functions.
3798 * Symbolic differentiation::
3799 * Series expansion:: Taylor and Laurent expansion.
3801 * Built-in functions:: List of predefined mathematical functions.
3802 * Multiple polylogarithms::
3803 * Complex expressions::
3804 * Solving linear systems of equations::
3805 * Input/output:: Input and output of expressions.
3809 @node Information about expressions, Numerical evaluation, Methods and functions, Methods and functions
3810 @c node-name, next, previous, up
3811 @section Getting information about expressions
3813 @subsection Checking expression types
3814 @cindex @code{is_a<@dots{}>()}
3815 @cindex @code{is_exactly_a<@dots{}>()}
3816 @cindex @code{ex_to<@dots{}>()}
3817 @cindex Converting @code{ex} to other classes
3818 @cindex @code{info()}
3819 @cindex @code{return_type()}
3820 @cindex @code{return_type_tinfo()}
3822 Sometimes it's useful to check whether a given expression is a plain number,
3823 a sum, a polynomial with integer coefficients, or of some other specific type.
3824 GiNaC provides a couple of functions for this:
3827 bool is_a<T>(const ex & e);
3828 bool is_exactly_a<T>(const ex & e);
3829 bool ex::info(unsigned flag);
3830 unsigned ex::return_type() const;
3831 unsigned ex::return_type_tinfo() const;
3834 When the test made by @code{is_a<T>()} returns true, it is safe to call
3835 one of the functions @code{ex_to<T>()}, where @code{T} is one of the
3836 class names (@xref{The class hierarchy}, for a list of all classes). For
3837 example, assuming @code{e} is an @code{ex}:
3842 if (is_a<numeric>(e))
3843 numeric n = ex_to<numeric>(e);
3848 @code{is_a<T>(e)} allows you to check whether the top-level object of
3849 an expression @samp{e} is an instance of the GiNaC class @samp{T}
3850 (@xref{The class hierarchy}, for a list of all classes). This is most useful,
3851 e.g., for checking whether an expression is a number, a sum, or a product:
3858 is_a<numeric>(e1); // true
3859 is_a<numeric>(e2); // false
3860 is_a<add>(e1); // false
3861 is_a<add>(e2); // true
3862 is_a<mul>(e1); // false
3863 is_a<mul>(e2); // false
3867 In contrast, @code{is_exactly_a<T>(e)} allows you to check whether the
3868 top-level object of an expression @samp{e} is an instance of the GiNaC
3869 class @samp{T}, not including parent classes.
3871 The @code{info()} method is used for checking certain attributes of
3872 expressions. The possible values for the @code{flag} argument are defined
3873 in @file{ginac/flags.h}, the most important being explained in the following
3877 @multitable @columnfractions .30 .70
3878 @item @strong{Flag} @tab @strong{Returns true if the object is@dots{}}
3879 @item @code{numeric}
3880 @tab @dots{}a number (same as @code{is_a<numeric>(...)})
3882 @tab @dots{}a real number, symbol or constant (i.e. is not complex)
3883 @item @code{rational}
3884 @tab @dots{}an exact rational number (integers are rational, too)
3885 @item @code{integer}
3886 @tab @dots{}a (non-complex) integer
3887 @item @code{crational}
3888 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
3889 @item @code{cinteger}
3890 @tab @dots{}a (complex) integer (such as @math{2-3*I})
3891 @item @code{positive}
3892 @tab @dots{}not complex and greater than 0
3893 @item @code{negative}
3894 @tab @dots{}not complex and less than 0
3895 @item @code{nonnegative}
3896 @tab @dots{}not complex and greater than or equal to 0
3898 @tab @dots{}an integer greater than 0
3900 @tab @dots{}an integer less than 0
3901 @item @code{nonnegint}
3902 @tab @dots{}an integer greater than or equal to 0
3904 @tab @dots{}an even integer
3906 @tab @dots{}an odd integer
3908 @tab @dots{}a prime integer (probabilistic primality test)
3909 @item @code{relation}
3910 @tab @dots{}a relation (same as @code{is_a<relational>(...)})
3911 @item @code{relation_equal}
3912 @tab @dots{}a @code{==} relation
3913 @item @code{relation_not_equal}
3914 @tab @dots{}a @code{!=} relation
3915 @item @code{relation_less}
3916 @tab @dots{}a @code{<} relation
3917 @item @code{relation_less_or_equal}
3918 @tab @dots{}a @code{<=} relation
3919 @item @code{relation_greater}
3920 @tab @dots{}a @code{>} relation
3921 @item @code{relation_greater_or_equal}
3922 @tab @dots{}a @code{>=} relation
3924 @tab @dots{}a symbol (same as @code{is_a<symbol>(...)})
3926 @tab @dots{}a list (same as @code{is_a<lst>(...)})
3927 @item @code{polynomial}
3928 @tab @dots{}a polynomial (i.e. only consists of sums and products of numbers and symbols with positive integer powers)
3929 @item @code{integer_polynomial}
3930 @tab @dots{}a polynomial with (non-complex) integer coefficients
3931 @item @code{cinteger_polynomial}
3932 @tab @dots{}a polynomial with (possibly complex) integer coefficients (such as @math{2-3*I})
3933 @item @code{rational_polynomial}
3934 @tab @dots{}a polynomial with (non-complex) rational coefficients
3935 @item @code{crational_polynomial}
3936 @tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I})
3937 @item @code{rational_function}
3938 @tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)})
3939 @item @code{algebraic}
3940 @tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1})
3944 To determine whether an expression is commutative or non-commutative and if
3945 so, with which other expressions it would commutate, you use the methods
3946 @code{return_type()} and @code{return_type_tinfo()}. @xref{Non-commutative objects},
3947 for an explanation of these.
3950 @subsection Accessing subexpressions
3953 Many GiNaC classes, like @code{add}, @code{mul}, @code{lst}, and
3954 @code{function}, act as containers for subexpressions. For example, the
3955 subexpressions of a sum (an @code{add} object) are the individual terms,
3956 and the subexpressions of a @code{function} are the function's arguments.
3958 @cindex @code{nops()}
3960 GiNaC provides several ways of accessing subexpressions. The first way is to
3965 ex ex::op(size_t i);
3968 @code{nops()} determines the number of subexpressions (operands) contained
3969 in the expression, while @code{op(i)} returns the @code{i}-th
3970 (0..@code{nops()-1}) subexpression. In the case of a @code{power} object,
3971 @code{op(0)} will return the basis and @code{op(1)} the exponent. For
3972 @code{indexed} objects, @code{op(0)} is the base expression and @code{op(i)},
3973 @math{i>0} are the indices.
3976 @cindex @code{const_iterator}
3977 The second way to access subexpressions is via the STL-style random-access
3978 iterator class @code{const_iterator} and the methods
3981 const_iterator ex::begin();
3982 const_iterator ex::end();
3985 @code{begin()} returns an iterator referring to the first subexpression;
3986 @code{end()} returns an iterator which is one-past the last subexpression.
3987 If the expression has no subexpressions, then @code{begin() == end()}. These
3988 iterators can also be used in conjunction with non-modifying STL algorithms.
3990 Here is an example that (non-recursively) prints the subexpressions of a
3991 given expression in three different ways:
3998 for (size_t i = 0; i != e.nops(); ++i)
3999 cout << e.op(i) << endl;
4002 for (const_iterator i = e.begin(); i != e.end(); ++i)
4005 // with iterators and STL copy()
4006 std::copy(e.begin(), e.end(), std::ostream_iterator<ex>(cout, "\n"));
4010 @cindex @code{const_preorder_iterator}
4011 @cindex @code{const_postorder_iterator}
4012 @code{op()}/@code{nops()} and @code{const_iterator} only access an
4013 expression's immediate children. GiNaC provides two additional iterator
4014 classes, @code{const_preorder_iterator} and @code{const_postorder_iterator},
4015 that iterate over all objects in an expression tree, in preorder or postorder,
4016 respectively. They are STL-style forward iterators, and are created with the
4020 const_preorder_iterator ex::preorder_begin();
4021 const_preorder_iterator ex::preorder_end();
4022 const_postorder_iterator ex::postorder_begin();
4023 const_postorder_iterator ex::postorder_end();
4026 The following example illustrates the differences between
4027 @code{const_iterator}, @code{const_preorder_iterator}, and
4028 @code{const_postorder_iterator}:
4032 symbol A("A"), B("B"), C("C");
4033 ex e = lst(lst(A, B), C);
4035 std::copy(e.begin(), e.end(),
4036 std::ostream_iterator<ex>(cout, "\n"));
4040 std::copy(e.preorder_begin(), e.preorder_end(),
4041 std::ostream_iterator<ex>(cout, "\n"));
4048 std::copy(e.postorder_begin(), e.postorder_end(),
4049 std::ostream_iterator<ex>(cout, "\n"));
4058 @cindex @code{relational} (class)
4059 Finally, the left-hand side and right-hand side expressions of objects of
4060 class @code{relational} (and only of these) can also be accessed with the
4069 @subsection Comparing expressions
4070 @cindex @code{is_equal()}
4071 @cindex @code{is_zero()}
4073 Expressions can be compared with the usual C++ relational operators like
4074 @code{==}, @code{>}, and @code{<} but if the expressions contain symbols,
4075 the result is usually not determinable and the result will be @code{false},
4076 except in the case of the @code{!=} operator. You should also be aware that
4077 GiNaC will only do the most trivial test for equality (subtracting both
4078 expressions), so something like @code{(pow(x,2)+x)/x==x+1} will return
4081 Actually, if you construct an expression like @code{a == b}, this will be
4082 represented by an object of the @code{relational} class (@pxref{Relations})
4083 which is not evaluated until (explicitly or implicitly) cast to a @code{bool}.
4085 There are also two methods
4088 bool ex::is_equal(const ex & other);
4092 for checking whether one expression is equal to another, or equal to zero,
4093 respectively. See also the method @code{ex::is_zero_matrix()},
4097 @subsection Ordering expressions
4098 @cindex @code{ex_is_less} (class)
4099 @cindex @code{ex_is_equal} (class)
4100 @cindex @code{compare()}
4102 Sometimes it is necessary to establish a mathematically well-defined ordering
4103 on a set of arbitrary expressions, for example to use expressions as keys
4104 in a @code{std::map<>} container, or to bring a vector of expressions into
4105 a canonical order (which is done internally by GiNaC for sums and products).
4107 The operators @code{<}, @code{>} etc. described in the last section cannot
4108 be used for this, as they don't implement an ordering relation in the
4109 mathematical sense. In particular, they are not guaranteed to be
4110 antisymmetric: if @samp{a} and @samp{b} are different expressions, and
4111 @code{a < b} yields @code{false}, then @code{b < a} doesn't necessarily
4114 By default, STL classes and algorithms use the @code{<} and @code{==}
4115 operators to compare objects, which are unsuitable for expressions, but GiNaC
4116 provides two functors that can be supplied as proper binary comparison
4117 predicates to the STL:
4120 class ex_is_less : public std::binary_function<ex, ex, bool> @{
4122 bool operator()(const ex &lh, const ex &rh) const;
4125 class ex_is_equal : public std::binary_function<ex, ex, bool> @{
4127 bool operator()(const ex &lh, const ex &rh) const;
4131 For example, to define a @code{map} that maps expressions to strings you
4135 std::map<ex, std::string, ex_is_less> myMap;
4138 Omitting the @code{ex_is_less} template parameter will introduce spurious
4139 bugs because the map operates improperly.
4141 Other examples for the use of the functors:
4149 std::sort(v.begin(), v.end(), ex_is_less());
4151 // count the number of expressions equal to '1'
4152 unsigned num_ones = std::count_if(v.begin(), v.end(),
4153 std::bind2nd(ex_is_equal(), 1));
4156 The implementation of @code{ex_is_less} uses the member function
4159 int ex::compare(const ex & other) const;
4162 which returns @math{0} if @code{*this} and @code{other} are equal, @math{-1}
4163 if @code{*this} sorts before @code{other}, and @math{1} if @code{*this} sorts
4167 @node Numerical evaluation, Substituting expressions, Information about expressions, Methods and functions
4168 @c node-name, next, previous, up
4169 @section Numerical evaluation
4170 @cindex @code{evalf()}
4172 GiNaC keeps algebraic expressions, numbers and constants in their exact form.
4173 To evaluate them using floating-point arithmetic you need to call
4176 ex ex::evalf(int level = 0) const;
4179 @cindex @code{Digits}
4180 The accuracy of the evaluation is controlled by the global object @code{Digits}
4181 which can be assigned an integer value. The default value of @code{Digits}
4182 is 17. @xref{Numbers}, for more information and examples.
4184 To evaluate an expression to a @code{double} floating-point number you can
4185 call @code{evalf()} followed by @code{numeric::to_double()}, like this:
4189 // Approximate sin(x/Pi)
4191 ex e = series(sin(x/Pi), x == 0, 6);
4193 // Evaluate numerically at x=0.1
4194 ex f = evalf(e.subs(x == 0.1));
4196 // ex_to<numeric> is an unsafe cast, so check the type first
4197 if (is_a<numeric>(f)) @{
4198 double d = ex_to<numeric>(f).to_double();
4207 @node Substituting expressions, Pattern matching and advanced substitutions, Numerical evaluation, Methods and functions
4208 @c node-name, next, previous, up
4209 @section Substituting expressions
4210 @cindex @code{subs()}
4212 Algebraic objects inside expressions can be replaced with arbitrary
4213 expressions via the @code{.subs()} method:
4216 ex ex::subs(const ex & e, unsigned options = 0);
4217 ex ex::subs(const exmap & m, unsigned options = 0);
4218 ex ex::subs(const lst & syms, const lst & repls, unsigned options = 0);
4221 In the first form, @code{subs()} accepts a relational of the form
4222 @samp{object == expression} or a @code{lst} of such relationals:
4226 symbol x("x"), y("y");
4228 ex e1 = 2*x^2-4*x+3;
4229 cout << "e1(7) = " << e1.subs(x == 7) << endl;
4233 cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
4238 If you specify multiple substitutions, they are performed in parallel, so e.g.
4239 @code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
4241 The second form of @code{subs()} takes an @code{exmap} object which is a
4242 pair associative container that maps expressions to expressions (currently
4243 implemented as a @code{std::map}). This is the most efficient one of the
4244 three @code{subs()} forms and should be used when the number of objects to
4245 be substituted is large or unknown.
4247 Using this form, the second example from above would look like this:
4251 symbol x("x"), y("y");
4257 cout << "e2(-2, 4) = " << e2.subs(m) << endl;
4261 The third form of @code{subs()} takes two lists, one for the objects to be
4262 replaced and one for the expressions to be substituted (both lists must
4263 contain the same number of elements). Using this form, you would write
4267 symbol x("x"), y("y");
4270 cout << "e2(-2, 4) = " << e2.subs(lst(x, y), lst(-2, 4)) << endl;
4274 The optional last argument to @code{subs()} is a combination of
4275 @code{subs_options} flags. There are three options available:
4276 @code{subs_options::no_pattern} disables pattern matching, which makes
4277 large @code{subs()} operations significantly faster if you are not using
4278 patterns. The second option, @code{subs_options::algebraic} enables
4279 algebraic substitutions in products and powers.
4280 @ref{Pattern matching and advanced substitutions}, for more information
4281 about patterns and algebraic substitutions. The third option,
4282 @code{subs_options::no_index_renaming} disables the feature that dummy
4283 indices are renamed if the subsitution could give a result in which a
4284 dummy index occurs more than two times. This is sometimes necessary if
4285 you want to use @code{subs()} to rename your dummy indices.
4287 @code{subs()} performs syntactic substitution of any complete algebraic
4288 object; it does not try to match sub-expressions as is demonstrated by the
4293 symbol x("x"), y("y"), z("z");
4295 ex e1 = pow(x+y, 2);
4296 cout << e1.subs(x+y == 4) << endl;
4299 ex e2 = sin(x)*sin(y)*cos(x);
4300 cout << e2.subs(sin(x) == cos(x)) << endl;
4301 // -> cos(x)^2*sin(y)
4304 cout << e3.subs(x+y == 4) << endl;
4306 // (and not 4+z as one might expect)
4310 A more powerful form of substitution using wildcards is described in the
4314 @node Pattern matching and advanced substitutions, Applying a function on subexpressions, Substituting expressions, Methods and functions
4315 @c node-name, next, previous, up
4316 @section Pattern matching and advanced substitutions
4317 @cindex @code{wildcard} (class)
4318 @cindex Pattern matching
4320 GiNaC allows the use of patterns for checking whether an expression is of a
4321 certain form or contains subexpressions of a certain form, and for
4322 substituting expressions in a more general way.
4324 A @dfn{pattern} is an algebraic expression that optionally contains wildcards.
4325 A @dfn{wildcard} is a special kind of object (of class @code{wildcard}) that
4326 represents an arbitrary expression. Every wildcard has a @dfn{label} which is
4327 an unsigned integer number to allow having multiple different wildcards in a
4328 pattern. Wildcards are printed as @samp{$label} (this is also the way they
4329 are specified in @command{ginsh}). In C++ code, wildcard objects are created
4333 ex wild(unsigned label = 0);
4336 which is simply a wrapper for the @code{wildcard()} constructor with a shorter
4339 Some examples for patterns:
4341 @multitable @columnfractions .5 .5
4342 @item @strong{Constructed as} @tab @strong{Output as}
4343 @item @code{wild()} @tab @samp{$0}
4344 @item @code{pow(x,wild())} @tab @samp{x^$0}
4345 @item @code{atan2(wild(1),wild(2))} @tab @samp{atan2($1,$2)}
4346 @item @code{indexed(A,idx(wild(),3))} @tab @samp{A.$0}
4352 @item Wildcards behave like symbols and are subject to the same algebraic
4353 rules. E.g., @samp{$0+2*$0} is automatically transformed to @samp{3*$0}.
4354 @item As shown in the last example, to use wildcards for indices you have to
4355 use them as the value of an @code{idx} object. This is because indices must
4356 always be of class @code{idx} (or a subclass).
4357 @item Wildcards only represent expressions or subexpressions. It is not
4358 possible to use them as placeholders for other properties like index
4359 dimension or variance, representation labels, symmetry of indexed objects
4361 @item Because wildcards are commutative, it is not possible to use wildcards
4362 as part of noncommutative products.
4363 @item A pattern does not have to contain wildcards. @samp{x} and @samp{x+y}
4364 are also valid patterns.
4367 @subsection Matching expressions
4368 @cindex @code{match()}
4369 The most basic application of patterns is to check whether an expression
4370 matches a given pattern. This is done by the function
4373 bool ex::match(const ex & pattern);
4374 bool ex::match(const ex & pattern, lst & repls);
4377 This function returns @code{true} when the expression matches the pattern
4378 and @code{false} if it doesn't. If used in the second form, the actual
4379 subexpressions matched by the wildcards get returned in the @code{repls}
4380 object as a list of relations of the form @samp{wildcard == expression}.
4381 If @code{match()} returns false, the state of @code{repls} is undefined.
4382 For reproducible results, the list should be empty when passed to
4383 @code{match()}, but it is also possible to find similarities in multiple
4384 expressions by passing in the result of a previous match.
4386 The matching algorithm works as follows:
4389 @item A single wildcard matches any expression. If one wildcard appears
4390 multiple times in a pattern, it must match the same expression in all
4391 places (e.g. @samp{$0} matches anything, and @samp{$0*($0+1)} matches
4392 @samp{x*(x+1)} but not @samp{x*(y+1)}).
4393 @item If the expression is not of the same class as the pattern, the match
4394 fails (i.e. a sum only matches a sum, a function only matches a function,
4396 @item If the pattern is a function, it only matches the same function
4397 (i.e. @samp{sin($0)} matches @samp{sin(x)} but doesn't match @samp{exp(x)}).
4398 @item Except for sums and products, the match fails if the number of
4399 subexpressions (@code{nops()}) is not equal to the number of subexpressions
4401 @item If there are no subexpressions, the expressions and the pattern must
4402 be equal (in the sense of @code{is_equal()}).
4403 @item Except for sums and products, each subexpression (@code{op()}) must
4404 match the corresponding subexpression of the pattern.
4407 Sums (@code{add}) and products (@code{mul}) are treated in a special way to
4408 account for their commutativity and associativity:
4411 @item If the pattern contains a term or factor that is a single wildcard,
4412 this one is used as the @dfn{global wildcard}. If there is more than one
4413 such wildcard, one of them is chosen as the global wildcard in a random
4415 @item Every term/factor of the pattern, except the global wildcard, is
4416 matched against every term of the expression in sequence. If no match is
4417 found, the whole match fails. Terms that did match are not considered in
4419 @item If there are no unmatched terms left, the match succeeds. Otherwise
4420 the match fails unless there is a global wildcard in the pattern, in
4421 which case this wildcard matches the remaining terms.
4424 In general, having more than one single wildcard as a term of a sum or a
4425 factor of a product (such as @samp{a+$0+$1}) will lead to unpredictable or
4428 Here are some examples in @command{ginsh} to demonstrate how it works (the
4429 @code{match()} function in @command{ginsh} returns @samp{FAIL} if the
4430 match fails, and the list of wildcard replacements otherwise):
4433 > match((x+y)^a,(x+y)^a);
4435 > match((x+y)^a,(x+y)^b);
4437 > match((x+y)^a,$1^$2);
4439 > match((x+y)^a,$1^$1);
4441 > match((x+y)^(x+y),$1^$1);
4443 > match((x+y)^(x+y),$1^$2);
4445 > match((a+b)*(a+c),($1+b)*($1+c));
4447 > match((a+b)*(a+c),(a+$1)*(a+$2));
4449 (Unpredictable. The result might also be [$1==c,$2==b].)
4450 > match((a+b)*(a+c),($1+$2)*($1+$3));
4451 (The result is undefined. Due to the sequential nature of the algorithm
4452 and the re-ordering of terms in GiNaC, the match for the first factor
4453 may be @{$1==a,$2==b@} in which case the match for the second factor
4454 succeeds, or it may be @{$1==b,$2==a@} which causes the second match to
4456 > match(a*(x+y)+a*z+b,a*$1+$2);
4457 (This is also ambiguous and may return either @{$1==z,$2==a*(x+y)+b@} or
4458 @{$1=x+y,$2=a*z+b@}.)
4459 > match(a+b+c+d+e+f,c);
4461 > match(a+b+c+d+e+f,c+$0);
4463 > match(a+b+c+d+e+f,c+e+$0);
4465 > match(a+b,a+b+$0);
4467 > match(a*b^2,a^$1*b^$2);
4469 (The matching is syntactic, not algebraic, and "a" doesn't match "a^$1"
4470 even though a==a^1.)
4471 > match(x*atan2(x,x^2),$0*atan2($0,$0^2));
4473 > match(atan2(y,x^2),atan2(y,$0));
4477 @subsection Matching parts of expressions
4478 @cindex @code{has()}
4479 A more general way to look for patterns in expressions is provided by the
4483 bool ex::has(const ex & pattern);
4486 This function checks whether a pattern is matched by an expression itself or
4487 by any of its subexpressions.
4489 Again some examples in @command{ginsh} for illustration (in @command{ginsh},
4490 @code{has()} returns @samp{1} for @code{true} and @samp{0} for @code{false}):
4493 > has(x*sin(x+y+2*a),y);
4495 > has(x*sin(x+y+2*a),x+y);
4497 (This is because in GiNaC, "x+y" is not a subexpression of "x+y+2*a" (which
4498 has the subexpressions "x", "y" and "2*a".)
4499 > has(x*sin(x+y+2*a),x+y+$1);
4501 (But this is possible.)
4502 > has(x*sin(2*(x+y)+2*a),x+y);
4504 (This fails because "2*(x+y)" automatically gets converted to "2*x+2*y" of
4505 which "x+y" is not a subexpression.)
4508 (Although x^1==x and x^0==1, neither "x" nor "1" are actually of the form
4510 > has(4*x^2-x+3,$1*x);
4512 > has(4*x^2+x+3,$1*x);
4514 (Another possible pitfall. The first expression matches because the term
4515 "-x" has the form "(-1)*x" in GiNaC. To check whether a polynomial
4516 contains a linear term you should use the coeff() function instead.)
4519 @cindex @code{find()}
4523 bool ex::find(const ex & pattern, lst & found);
4526 works a bit like @code{has()} but it doesn't stop upon finding the first
4527 match. Instead, it appends all found matches to the specified list. If there
4528 are multiple occurrences of the same expression, it is entered only once to
4529 the list. @code{find()} returns false if no matches were found (in
4530 @command{ginsh}, it returns an empty list):
4533 > find(1+x+x^2+x^3,x);
4535 > find(1+x+x^2+x^3,y);
4537 > find(1+x+x^2+x^3,x^$1);
4539 (Note the absence of "x".)
4540 > expand((sin(x)+sin(y))*(a+b));
4541 sin(y)*a+sin(x)*b+sin(x)*a+sin(y)*b
4546 @subsection Substituting expressions
4547 @cindex @code{subs()}
4548 Probably the most useful application of patterns is to use them for
4549 substituting expressions with the @code{subs()} method. Wildcards can be
4550 used in the search patterns as well as in the replacement expressions, where
4551 they get replaced by the expressions matched by them. @code{subs()} doesn't
4552 know anything about algebra; it performs purely syntactic substitutions.
4557 > subs(a^2+b^2+(x+y)^2,$1^2==$1^3);
4559 > subs(a^4+b^4+(x+y)^4,$1^2==$1^3);
4561 > subs((a+b+c)^2,a+b==x);
4563 > subs((a+b+c)^2,a+b+$1==x+$1);
4565 > subs(a+2*b,a+b==x);
4567 > subs(4*x^3-2*x^2+5*x-1,x==a);
4569 > subs(4*x^3-2*x^2+5*x-1,x^$0==a^$0);
4571 > subs(sin(1+sin(x)),sin($1)==cos($1));
4573 > expand(subs(a*sin(x+y)^2+a*cos(x+y)^2+b,cos($1)^2==1-sin($1)^2));
4577 The last example would be written in C++ in this way:
4581 symbol a("a"), b("b"), x("x"), y("y");
4582 e = a*pow(sin(x+y), 2) + a*pow(cos(x+y), 2) + b;
4583 e = e.subs(pow(cos(wild()), 2) == 1-pow(sin(wild()), 2));
4584 cout << e.expand() << endl;
4589 @subsection The option algebraic
4590 Both @code{has()} and @code{subs()} take an optional argument to pass them
4591 extra options. This section describes what happens if you give the former
4592 the option @code{has_options::algebraic} or the latter
4593 @code{subs:options::algebraic}. In that case the matching condition for
4594 powers and multiplications is changed in such a way that they become
4595 more intuitive. Intuition says that @code{x*y} is a part of @code{x*y*z}.
4596 If you use these options you will find that
4597 @code{(x*y*z).has(x*y, has_options::algebraic)} indeed returns true.
4598 Besides matching some of the factors of a product also powers match as
4599 often as is possible without getting negative exponents. For example
4600 @code{(x^5*y^2*z).subs(x^2*y^2==c, subs_options::algebraic)} will return
4601 @code{x*c^2*z}. This also works with negative powers:
4602 @code{(x^(-3)*y^(-2)*z).subs(1/(x*y)==c, subs_options::algebraic)} will
4603 return @code{x^(-1)*c^2*z}. Note that this only works for multiplications
4604 and not for locating @code{x+y} within @code{x+y+z}.
4607 @node Applying a function on subexpressions, Visitors and tree traversal, Pattern matching and advanced substitutions, Methods and functions
4608 @c node-name, next, previous, up
4609 @section Applying a function on subexpressions
4610 @cindex tree traversal
4611 @cindex @code{map()}
4613 Sometimes you may want to perform an operation on specific parts of an
4614 expression while leaving the general structure of it intact. An example
4615 of this would be a matrix trace operation: the trace of a sum is the sum
4616 of the traces of the individual terms. That is, the trace should @dfn{map}
4617 on the sum, by applying itself to each of the sum's operands. It is possible
4618 to do this manually which usually results in code like this:
4623 if (is_a<matrix>(e))
4624 return ex_to<matrix>(e).trace();
4625 else if (is_a<add>(e)) @{
4627 for (size_t i=0; i<e.nops(); i++)
4628 sum += calc_trace(e.op(i));
4630 @} else if (is_a<mul>)(e)) @{
4638 This is, however, slightly inefficient (if the sum is very large it can take
4639 a long time to add the terms one-by-one), and its applicability is limited to
4640 a rather small class of expressions. If @code{calc_trace()} is called with
4641 a relation or a list as its argument, you will probably want the trace to
4642 be taken on both sides of the relation or of all elements of the list.
4644 GiNaC offers the @code{map()} method to aid in the implementation of such
4648 ex ex::map(map_function & f) const;
4649 ex ex::map(ex (*f)(const ex & e)) const;
4652 In the first (preferred) form, @code{map()} takes a function object that
4653 is subclassed from the @code{map_function} class. In the second form, it
4654 takes a pointer to a function that accepts and returns an expression.
4655 @code{map()} constructs a new expression of the same type, applying the
4656 specified function on all subexpressions (in the sense of @code{op()}),
4659 The use of a function object makes it possible to supply more arguments to
4660 the function that is being mapped, or to keep local state information.
4661 The @code{map_function} class declares a virtual function call operator
4662 that you can overload. Here is a sample implementation of @code{calc_trace()}
4663 that uses @code{map()} in a recursive fashion:
4666 struct calc_trace : public map_function @{
4667 ex operator()(const ex &e)
4669 if (is_a<matrix>(e))
4670 return ex_to<matrix>(e).trace();
4671 else if (is_a<mul>(e)) @{
4674 return e.map(*this);
4679 This function object could then be used like this:
4683 ex M = ... // expression with matrices
4684 calc_trace do_trace;
4685 ex tr = do_trace(M);
4689 Here is another example for you to meditate over. It removes quadratic
4690 terms in a variable from an expanded polynomial:
4693 struct map_rem_quad : public map_function @{
4695 map_rem_quad(const ex & var_) : var(var_) @{@}
4697 ex operator()(const ex & e)
4699 if (is_a<add>(e) || is_a<mul>(e))
4700 return e.map(*this);
4701 else if (is_a<power>(e) &&
4702 e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
4712 symbol x("x"), y("y");
4715 for (int i=0; i<8; i++)
4716 e += pow(x, i) * pow(y, 8-i) * (i+1);
4718 // -> 4*y^5*x^3+5*y^4*x^4+8*y*x^7+7*y^2*x^6+2*y^7*x+6*y^3*x^5+3*y^6*x^2+y^8
4720 map_rem_quad rem_quad(x);
4721 cout << rem_quad(e) << endl;
4722 // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
4726 @command{ginsh} offers a slightly different implementation of @code{map()}
4727 that allows applying algebraic functions to operands. The second argument
4728 to @code{map()} is an expression containing the wildcard @samp{$0} which
4729 acts as the placeholder for the operands:
4734 > map(a+2*b,sin($0));
4736 > map(@{a,b,c@},$0^2+$0);
4737 @{a^2+a,b^2+b,c^2+c@}
4740 Note that it is only possible to use algebraic functions in the second
4741 argument. You can not use functions like @samp{diff()}, @samp{op()},
4742 @samp{subs()} etc. because these are evaluated immediately:
4745 > map(@{a,b,c@},diff($0,a));
4747 This is because "diff($0,a)" evaluates to "0", so the command is equivalent
4748 to "map(@{a,b,c@},0)".
4752 @node Visitors and tree traversal, Polynomial arithmetic, Applying a function on subexpressions, Methods and functions
4753 @c node-name, next, previous, up
4754 @section Visitors and tree traversal
4755 @cindex tree traversal
4756 @cindex @code{visitor} (class)
4757 @cindex @code{accept()}
4758 @cindex @code{visit()}
4759 @cindex @code{traverse()}
4760 @cindex @code{traverse_preorder()}
4761 @cindex @code{traverse_postorder()}
4763 Suppose that you need a function that returns a list of all indices appearing
4764 in an arbitrary expression. The indices can have any dimension, and for
4765 indices with variance you always want the covariant version returned.
4767 You can't use @code{get_free_indices()} because you also want to include
4768 dummy indices in the list, and you can't use @code{find()} as it needs
4769 specific index dimensions (and it would require two passes: one for indices
4770 with variance, one for plain ones).
4772 The obvious solution to this problem is a tree traversal with a type switch,
4773 such as the following:
4776 void gather_indices_helper(const ex & e, lst & l)
4778 if (is_a<varidx>(e)) @{
4779 const varidx & vi = ex_to<varidx>(e);
4780 l.append(vi.is_covariant() ? vi : vi.toggle_variance());
4781 @} else if (is_a<idx>(e)) @{
4784 size_t n = e.nops();
4785 for (size_t i = 0; i < n; ++i)
4786 gather_indices_helper(e.op(i), l);
4790 lst gather_indices(const ex & e)
4793 gather_indices_helper(e, l);
4800 This works fine but fans of object-oriented programming will feel
4801 uncomfortable with the type switch. One reason is that there is a possibility
4802 for subtle bugs regarding derived classes. If we had, for example, written
4805 if (is_a<idx>(e)) @{
4807 @} else if (is_a<varidx>(e)) @{
4811 in @code{gather_indices_helper}, the code wouldn't have worked because the
4812 first line "absorbs" all classes derived from @code{idx}, including
4813 @code{varidx}, so the special case for @code{varidx} would never have been
4816 Also, for a large number of classes, a type switch like the above can get
4817 unwieldy and inefficient (it's a linear search, after all).
4818 @code{gather_indices_helper} only checks for two classes, but if you had to
4819 write a function that required a different implementation for nearly
4820 every GiNaC class, the result would be very hard to maintain and extend.
4822 The cleanest approach to the problem would be to add a new virtual function
4823 to GiNaC's class hierarchy. In our example, there would be specializations
4824 for @code{idx} and @code{varidx} while the default implementation in
4825 @code{basic} performed the tree traversal. Unfortunately, in C++ it's
4826 impossible to add virtual member functions to existing classes without
4827 changing their source and recompiling everything. GiNaC comes with source,
4828 so you could actually do this, but for a small algorithm like the one
4829 presented this would be impractical.
4831 One solution to this dilemma is the @dfn{Visitor} design pattern,
4832 which is implemented in GiNaC (actually, Robert Martin's Acyclic Visitor
4833 variation, described in detail in
4834 @uref{http://objectmentor.com/publications/acv.pdf}). Instead of adding
4835 virtual functions to the class hierarchy to implement operations, GiNaC
4836 provides a single "bouncing" method @code{accept()} that takes an instance
4837 of a special @code{visitor} class and redirects execution to the one
4838 @code{visit()} virtual function of the visitor that matches the type of
4839 object that @code{accept()} was being invoked on.
4841 Visitors in GiNaC must derive from the global @code{visitor} class as well
4842 as from the class @code{T::visitor} of each class @code{T} they want to
4843 visit, and implement the member functions @code{void visit(const T &)} for
4849 void ex::accept(visitor & v) const;
4852 will then dispatch to the correct @code{visit()} member function of the
4853 specified visitor @code{v} for the type of GiNaC object at the root of the
4854 expression tree (e.g. a @code{symbol}, an @code{idx} or a @code{mul}).
4856 Here is an example of a visitor:
4860 : public visitor, // this is required
4861 public add::visitor, // visit add objects
4862 public numeric::visitor, // visit numeric objects
4863 public basic::visitor // visit basic objects
4865 void visit(const add & x)
4866 @{ cout << "called with an add object" << endl; @}
4868 void visit(const numeric & x)
4869 @{ cout << "called with a numeric object" << endl; @}
4871 void visit(const basic & x)
4872 @{ cout << "called with a basic object" << endl; @}
4876 which can be used as follows:
4887 // prints "called with a numeric object"
4889 // prints "called with an add object"
4891 // prints "called with a basic object"
4895 The @code{visit(const basic &)} method gets called for all objects that are
4896 not @code{numeric} or @code{add} and acts as an (optional) default.
4898 From a conceptual point of view, the @code{visit()} methods of the visitor
4899 behave like a newly added virtual function of the visited hierarchy.
4900 In addition, visitors can store state in member variables, and they can
4901 be extended by deriving a new visitor from an existing one, thus building
4902 hierarchies of visitors.
4904 We can now rewrite our index example from above with a visitor:
4907 class gather_indices_visitor
4908 : public visitor, public idx::visitor, public varidx::visitor
4912 void visit(const idx & i)
4917 void visit(const varidx & vi)
4919 l.append(vi.is_covariant() ? vi : vi.toggle_variance());
4923 const lst & get_result() // utility function
4932 What's missing is the tree traversal. We could implement it in
4933 @code{visit(const basic &)}, but GiNaC has predefined methods for this:
4936 void ex::traverse_preorder(visitor & v) const;
4937 void ex::traverse_postorder(visitor & v) const;
4938 void ex::traverse(visitor & v) const;
4941 @code{traverse_preorder()} visits a node @emph{before} visiting its
4942 subexpressions, while @code{traverse_postorder()} visits a node @emph{after}
4943 visiting its subexpressions. @code{traverse()} is a synonym for
4944 @code{traverse_preorder()}.
4946 Here is a new implementation of @code{gather_indices()} that uses the visitor
4947 and @code{traverse()}:
4950 lst gather_indices(const ex & e)
4952 gather_indices_visitor v;
4954 return v.get_result();
4958 Alternatively, you could use pre- or postorder iterators for the tree
4962 lst gather_indices(const ex & e)
4964 gather_indices_visitor v;
4965 for (const_preorder_iterator i = e.preorder_begin();
4966 i != e.preorder_end(); ++i) @{
4969 return v.get_result();
4974 @node Polynomial arithmetic, Rational expressions, Visitors and tree traversal, Methods and functions
4975 @c node-name, next, previous, up
4976 @section Polynomial arithmetic
4978 @subsection Testing whether an expression is a polynomial
4979 @cindex @code{is_polynomial()}
4981 Testing whether an expression is a polynomial in one or more variables
4982 can be done with the method
4984 bool ex::is_polynomial(const ex & vars) const;
4986 In the case of more than
4987 one variable, the variables are given as a list.
4990 (x*y*sin(y)).is_polynomial(x) // Returns true.
4991 (x*y*sin(y)).is_polynomial(lst(x,y)) // Returns false.
4994 @subsection Expanding and collecting
4995 @cindex @code{expand()}
4996 @cindex @code{collect()}
4997 @cindex @code{collect_common_factors()}
4999 A polynomial in one or more variables has many equivalent
5000 representations. Some useful ones serve a specific purpose. Consider
5001 for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
5002 21*y*z + 4*z^2} (written down here in output-style). It is equivalent
5003 to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
5004 representations are the recursive ones where one collects for exponents
5005 in one of the three variable. Since the factors are themselves
5006 polynomials in the remaining two variables the procedure can be
5007 repeated. In our example, two possibilities would be @math{(4*y + z)*x
5008 + 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
5011 To bring an expression into expanded form, its method
5014 ex ex::expand(unsigned options = 0);
5017 may be called. In our example above, this corresponds to @math{4*x*y +
5018 x*z + 20*y^2 + 21*y*z + 4*z^2}. Again, since the canonical form in
5019 GiNaC is not easy to guess you should be prepared to see different
5020 orderings of terms in such sums!
5022 Another useful representation of multivariate polynomials is as a
5023 univariate polynomial in one of the variables with the coefficients
5024 being polynomials in the remaining variables. The method
5025 @code{collect()} accomplishes this task:
5028 ex ex::collect(const ex & s, bool distributed = false);
5031 The first argument to @code{collect()} can also be a list of objects in which
5032 case the result is either a recursively collected polynomial, or a polynomial
5033 in a distributed form with terms like @math{c*x1^e1*...*xn^en}, as specified
5034 by the @code{distributed} flag.
5036 Note that the original polynomial needs to be in expanded form (for the
5037 variables concerned) in order for @code{collect()} to be able to find the
5038 coefficients properly.
5040 The following @command{ginsh} transcript shows an application of @code{collect()}
5041 together with @code{find()}:
5044 > a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
5045 d*p*sin(x)+p*sin(x)+q*d*sin(x)+q*sin(y)+d*sin(x)+q*d*sin(y)+sin(y)+d*sin(y)
5046 +q*sin(x)+d*sin(y)*p+sin(x)+sin(y)*p
5047 > collect(a,@{p,q@});
5048 d*sin(x)+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*p
5049 +(d*sin(x)+sin(y)+d*sin(y)+sin(x))*q+sin(y)+d*sin(y)+sin(x)
5050 > collect(a,find(a,sin($1)));
5051 (1+q+d+q*d+d*p+p)*sin(y)+(1+q+d+q*d+d*p+p)*sin(x)
5052 > collect(a,@{find(a,sin($1)),p,q@});
5053 (1+(1+d)*p+d+q*(1+d))*sin(x)+(1+(1+d)*p+d+q*(1+d))*sin(y)
5054 > collect(a,@{find(a,sin($1)),d@});
5055 (1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
5058 Polynomials can often be brought into a more compact form by collecting
5059 common factors from the terms of sums. This is accomplished by the function
5062 ex collect_common_factors(const ex & e);
5065 This function doesn't perform a full factorization but only looks for
5066 factors which are already explicitly present:
5069 > collect_common_factors(a*x+a*y);
5071 > collect_common_factors(a*x^2+2*a*x*y+a*y^2);
5073 > collect_common_factors(a*(b*(a+c)*x+b*((a+c)*x+(a+c)*y)*y));
5074 (c+a)*a*(x*y+y^2+x)*b
5077 @subsection Degree and coefficients
5078 @cindex @code{degree()}
5079 @cindex @code{ldegree()}
5080 @cindex @code{coeff()}
5082 The degree and low degree of a polynomial can be obtained using the two
5086 int ex::degree(const ex & s);
5087 int ex::ldegree(const ex & s);
5090 which also work reliably on non-expanded input polynomials (they even work
5091 on rational functions, returning the asymptotic degree). By definition, the
5092 degree of zero is zero. To extract a coefficient with a certain power from
5093 an expanded polynomial you use
5096 ex ex::coeff(const ex & s, int n);
5099 You can also obtain the leading and trailing coefficients with the methods
5102 ex ex::lcoeff(const ex & s);
5103 ex ex::tcoeff(const ex & s);
5106 which are equivalent to @code{coeff(s, degree(s))} and @code{coeff(s, ldegree(s))},
5109 An application is illustrated in the next example, where a multivariate
5110 polynomial is analyzed:
5114 symbol x("x"), y("y");
5115 ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
5116 - pow(x+y,2) + 2*pow(y+2,2) - 8;
5117 ex Poly = PolyInp.expand();
5119 for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) @{
5120 cout << "The x^" << i << "-coefficient is "
5121 << Poly.coeff(x,i) << endl;
5123 cout << "As polynomial in y: "
5124 << Poly.collect(y) << endl;
5128 When run, it returns an output in the following fashion:
5131 The x^0-coefficient is y^2+11*y
5132 The x^1-coefficient is 5*y^2-2*y
5133 The x^2-coefficient is -1
5134 The x^3-coefficient is 4*y
5135 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
5138 As always, the exact output may vary between different versions of GiNaC
5139 or even from run to run since the internal canonical ordering is not
5140 within the user's sphere of influence.
5142 @code{degree()}, @code{ldegree()}, @code{coeff()}, @code{lcoeff()},
5143 @code{tcoeff()} and @code{collect()} can also be used to a certain degree
5144 with non-polynomial expressions as they not only work with symbols but with
5145 constants, functions and indexed objects as well:
5149 symbol a("a"), b("b"), c("c"), x("x");
5150 idx i(symbol("i"), 3);
5152 ex e = pow(sin(x) - cos(x), 4);
5153 cout << e.degree(cos(x)) << endl;
5155 cout << e.expand().coeff(sin(x), 3) << endl;
5158 e = indexed(a+b, i) * indexed(b+c, i);
5159 e = e.expand(expand_options::expand_indexed);
5160 cout << e.collect(indexed(b, i)) << endl;
5161 // -> a.i*c.i+(a.i+c.i)*b.i+b.i^2
5166 @subsection Polynomial division
5167 @cindex polynomial division
5170 @cindex pseudo-remainder
5171 @cindex @code{quo()}
5172 @cindex @code{rem()}
5173 @cindex @code{prem()}
5174 @cindex @code{divide()}
5179 ex quo(const ex & a, const ex & b, const ex & x);
5180 ex rem(const ex & a, const ex & b, const ex & x);
5183 compute the quotient and remainder of univariate polynomials in the variable
5184 @samp{x}. The results satisfy @math{a = b*quo(a, b, x) + rem(a, b, x)}.
5186 The additional function
5189 ex prem(const ex & a, const ex & b, const ex & x);
5192 computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
5193 @math{c*a = b*q + prem(a, b, x)}, where @math{c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1)}.
5195 Exact division of multivariate polynomials is performed by the function
5198 bool divide(const ex & a, const ex & b, ex & q);
5201 If @samp{b} divides @samp{a} over the rationals, this function returns @code{true}
5202 and returns the quotient in the variable @code{q}. Otherwise it returns @code{false}
5203 in which case the value of @code{q} is undefined.
5206 @subsection Unit, content and primitive part
5207 @cindex @code{unit()}
5208 @cindex @code{content()}
5209 @cindex @code{primpart()}
5210 @cindex @code{unitcontprim()}
5215 ex ex::unit(const ex & x);
5216 ex ex::content(const ex & x);
5217 ex ex::primpart(const ex & x);
5218 ex ex::primpart(const ex & x, const ex & c);
5221 return the unit part, content part, and primitive polynomial of a multivariate
5222 polynomial with respect to the variable @samp{x} (the unit part being the sign
5223 of the leading coefficient, the content part being the GCD of the coefficients,
5224 and the primitive polynomial being the input polynomial divided by the unit and
5225 content parts). The second variant of @code{primpart()} expects the previously
5226 calculated content part of the polynomial in @code{c}, which enables it to
5227 work faster in the case where the content part has already been computed. The
5228 product of unit, content, and primitive part is the original polynomial.
5230 Additionally, the method
5233 void ex::unitcontprim(const ex & x, ex & u, ex & c, ex & p);
5236 computes the unit, content, and primitive parts in one go, returning them
5237 in @code{u}, @code{c}, and @code{p}, respectively.
5240 @subsection GCD, LCM and resultant
5243 @cindex @code{gcd()}
5244 @cindex @code{lcm()}
5246 The functions for polynomial greatest common divisor and least common
5247 multiple have the synopsis
5250 ex gcd(const ex & a, const ex & b);
5251 ex lcm(const ex & a, const ex & b);
5254 The functions @code{gcd()} and @code{lcm()} accept two expressions
5255 @code{a} and @code{b} as arguments and return a new expression, their
5256 greatest common divisor or least common multiple, respectively. If the
5257 polynomials @code{a} and @code{b} are coprime @code{gcd(a,b)} returns 1
5258 and @code{lcm(a,b)} returns the product of @code{a} and @code{b}. Note that all
5259 the coefficients must be rationals.
5262 #include <ginac/ginac.h>
5263 using namespace GiNaC;
5267 symbol x("x"), y("y"), z("z");
5268 ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
5269 ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
5271 ex P_gcd = gcd(P_a, P_b);
5273 ex P_lcm = lcm(P_a, P_b);
5274 // 4*x*y^2 + 13*y*x*z + 20*y^3 + 81*y^2*z + 67*y*z^2 + 3*x*z^2 + 12*z^3
5279 @cindex @code{resultant()}
5281 The resultant of two expressions only makes sense with polynomials.
5282 It is always computed with respect to a specific symbol within the
5283 expressions. The function has the interface
5286 ex resultant(const ex & a, const ex & b, const ex & s);
5289 Resultants are symmetric in @code{a} and @code{b}. The following example
5290 computes the resultant of two expressions with respect to @code{x} and
5291 @code{y}, respectively:
5294 #include <ginac/ginac.h>
5295 using namespace GiNaC;
5299 symbol x("x"), y("y");
5301 ex e1 = x+pow(y,2), e2 = 2*pow(x,3)-1; // x+y^2, 2*x^3-1
5304 r = resultant(e1, e2, x);
5306 r = resultant(e1, e2, y);
5311 @subsection Square-free decomposition
5312 @cindex square-free decomposition
5313 @cindex factorization
5314 @cindex @code{sqrfree()}
5316 GiNaC still lacks proper factorization support. Some form of
5317 factorization is, however, easily implemented by noting that factors
5318 appearing in a polynomial with power two or more also appear in the
5319 derivative and hence can easily be found by computing the GCD of the
5320 original polynomial and its derivatives. Any decent system has an
5321 interface for this so called square-free factorization. So we provide
5324 ex sqrfree(const ex & a, const lst & l = lst());
5326 Here is an example that by the way illustrates how the exact form of the
5327 result may slightly depend on the order of differentiation, calling for
5328 some care with subsequent processing of the result:
5331 symbol x("x"), y("y");
5332 ex BiVarPol = expand(pow(2-2*y,3) * pow(1+x*y,2) * pow(x-2*y,2) * (x+y));
5334 cout << sqrfree(BiVarPol, lst(x,y)) << endl;
5335 // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)
5337 cout << sqrfree(BiVarPol, lst(y,x)) << endl;
5338 // -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x)
5340 cout << sqrfree(BiVarPol) << endl;
5341 // -> depending on luck, any of the above
5344 Note also, how factors with the same exponents are not fully factorized
5348 @node Rational expressions, Symbolic differentiation, Polynomial arithmetic, Methods and functions
5349 @c node-name, next, previous, up
5350 @section Rational expressions
5352 @subsection The @code{normal} method
5353 @cindex @code{normal()}
5354 @cindex simplification
5355 @cindex temporary replacement
5357 Some basic form of simplification of expressions is called for frequently.
5358 GiNaC provides the method @code{.normal()}, which converts a rational function
5359 into an equivalent rational function of the form @samp{numerator/denominator}
5360 where numerator and denominator are coprime. If the input expression is already
5361 a fraction, it just finds the GCD of numerator and denominator and cancels it,
5362 otherwise it performs fraction addition and multiplication.
5364 @code{.normal()} can also be used on expressions which are not rational functions
5365 as it will replace all non-rational objects (like functions or non-integer
5366 powers) by temporary symbols to bring the expression to the domain of rational
5367 functions before performing the normalization, and re-substituting these
5368 symbols afterwards. This algorithm is also available as a separate method
5369 @code{.to_rational()}, described below.
5371 This means that both expressions @code{t1} and @code{t2} are indeed
5372 simplified in this little code snippet:
5377 ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
5378 ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
5379 std::cout << "t1 is " << t1.normal() << std::endl;
5380 std::cout << "t2 is " << t2.normal() << std::endl;
5384 Of course this works for multivariate polynomials too, so the ratio of
5385 the sample-polynomials from the section about GCD and LCM above would be
5386 normalized to @code{P_a/P_b} = @code{(4*y+z)/(y+3*z)}.
5389 @subsection Numerator and denominator
5392 @cindex @code{numer()}
5393 @cindex @code{denom()}
5394 @cindex @code{numer_denom()}
5396 The numerator and denominator of an expression can be obtained with
5401 ex ex::numer_denom();
5404 These functions will first normalize the expression as described above and
5405 then return the numerator, denominator, or both as a list, respectively.
5406 If you need both numerator and denominator, calling @code{numer_denom()} is
5407 faster than using @code{numer()} and @code{denom()} separately.
5410 @subsection Converting to a polynomial or rational expression
5411 @cindex @code{to_polynomial()}
5412 @cindex @code{to_rational()}
5414 Some of the methods described so far only work on polynomials or rational
5415 functions. GiNaC provides a way to extend the domain of these functions to
5416 general expressions by using the temporary replacement algorithm described
5417 above. You do this by calling
5420 ex ex::to_polynomial(exmap & m);
5421 ex ex::to_polynomial(lst & l);
5425 ex ex::to_rational(exmap & m);
5426 ex ex::to_rational(lst & l);
5429 on the expression to be converted. The supplied @code{exmap} or @code{lst}
5430 will be filled with the generated temporary symbols and their replacement
5431 expressions in a format that can be used directly for the @code{subs()}
5432 method. It can also already contain a list of replacements from an earlier
5433 application of @code{.to_polynomial()} or @code{.to_rational()}, so it's
5434 possible to use it on multiple expressions and get consistent results.
5436 The difference between @code{.to_polynomial()} and @code{.to_rational()}
5437 is probably best illustrated with an example:
5441 symbol x("x"), y("y");
5442 ex a = 2*x/sin(x) - y/(3*sin(x));
5446 ex p = a.to_polynomial(lp);
5447 cout << " = " << p << "\n with " << lp << endl;
5448 // = symbol3*symbol2*y+2*symbol2*x
5449 // with @{symbol2==sin(x)^(-1),symbol3==-1/3@}
5452 ex r = a.to_rational(lr);
5453 cout << " = " << r << "\n with " << lr << endl;
5454 // = -1/3*symbol4^(-1)*y+2*symbol4^(-1)*x
5455 // with @{symbol4==sin(x)@}
5459 The following more useful example will print @samp{sin(x)-cos(x)}:
5464 ex a = pow(sin(x), 2) - pow(cos(x), 2);
5465 ex b = sin(x) + cos(x);
5468 divide(a.to_polynomial(m), b.to_polynomial(m), q);
5469 cout << q.subs(m) << endl;
5474 @node Symbolic differentiation, Series expansion, Rational expressions, Methods and functions
5475 @c node-name, next, previous, up
5476 @section Symbolic differentiation
5477 @cindex differentiation
5478 @cindex @code{diff()}
5480 @cindex product rule
5482 GiNaC's objects know how to differentiate themselves. Thus, a
5483 polynomial (class @code{add}) knows that its derivative is the sum of
5484 the derivatives of all the monomials:
5488 symbol x("x"), y("y"), z("z");
5489 ex P = pow(x, 5) + pow(x, 2) + y;
5491 cout << P.diff(x,2) << endl;
5493 cout << P.diff(y) << endl; // 1
5495 cout << P.diff(z) << endl; // 0
5500 If a second integer parameter @var{n} is given, the @code{diff} method
5501 returns the @var{n}th derivative.
5503 If @emph{every} object and every function is told what its derivative
5504 is, all derivatives of composed objects can be calculated using the
5505 chain rule and the product rule. Consider, for instance the expression
5506 @code{1/cosh(x)}. Since the derivative of @code{cosh(x)} is
5507 @code{sinh(x)} and the derivative of @code{pow(x,-1)} is
5508 @code{-pow(x,-2)}, GiNaC can readily compute the composition. It turns
5509 out that the composition is the generating function for Euler Numbers,
5510 i.e. the so called @var{n}th Euler number is the coefficient of
5511 @code{x^n/n!} in the expansion of @code{1/cosh(x)}. We may use this
5512 identity to code a function that generates Euler numbers in just three
5515 @cindex Euler numbers
5517 #include <ginac/ginac.h>
5518 using namespace GiNaC;
5520 ex EulerNumber(unsigned n)
5523 const ex generator = pow(cosh(x),-1);
5524 return generator.diff(x,n).subs(x==0);
5529 for (unsigned i=0; i<11; i+=2)
5530 std::cout << EulerNumber(i) << std::endl;
5535 When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
5536 @code{-61}, @code{1385}, @code{-50521}. We increment the loop variable
5537 @code{i} by two since all odd Euler numbers vanish anyways.
5540 @node Series expansion, Symmetrization, Symbolic differentiation, Methods and functions
5541 @c node-name, next, previous, up
5542 @section Series expansion
5543 @cindex @code{series()}
5544 @cindex Taylor expansion
5545 @cindex Laurent expansion
5546 @cindex @code{pseries} (class)
5547 @cindex @code{Order()}
5549 Expressions know how to expand themselves as a Taylor series or (more
5550 generally) a Laurent series. As in most conventional Computer Algebra
5551 Systems, no distinction is made between those two. There is a class of
5552 its own for storing such series (@code{class pseries}) and a built-in
5553 function (called @code{Order}) for storing the order term of the series.
5554 As a consequence, if you want to work with series, i.e. multiply two
5555 series, you need to call the method @code{ex::series} again to convert
5556 it to a series object with the usual structure (expansion plus order
5557 term). A sample application from special relativity could read:
5560 #include <ginac/ginac.h>
5561 using namespace std;
5562 using namespace GiNaC;
5566 symbol v("v"), c("c");
5568 ex gamma = 1/sqrt(1 - pow(v/c,2));
5569 ex mass_nonrel = gamma.series(v==0, 10);
5571 cout << "the relativistic mass increase with v is " << endl
5572 << mass_nonrel << endl;
5574 cout << "the inverse square of this series is " << endl
5575 << pow(mass_nonrel,-2).series(v==0, 10) << endl;
5579 Only calling the series method makes the last output simplify to
5580 @math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
5581 series raised to the power @math{-2}.
5583 @cindex Machin's formula
5584 As another instructive application, let us calculate the numerical
5585 value of Archimedes' constant
5589 (for which there already exists the built-in constant @code{Pi})
5590 using John Machin's amazing formula
5592 $\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
5595 @math{Pi==16*atan(1/5)-4*atan(1/239)}.
5597 This equation (and similar ones) were used for over 200 years for
5598 computing digits of pi (see @cite{Pi Unleashed}). We may expand the
5599 arcus tangent around @code{0} and insert the fractions @code{1/5} and
5600 @code{1/239}. However, as we have seen, a series in GiNaC carries an
5601 order term with it and the question arises what the system is supposed
5602 to do when the fractions are plugged into that order term. The solution
5603 is to use the function @code{series_to_poly()} to simply strip the order
5607 #include <ginac/ginac.h>
5608 using namespace GiNaC;
5610 ex machin_pi(int degr)
5613 ex pi_expansion = series_to_poly(atan(x).series(x,degr));
5614 ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
5615 -4*pi_expansion.subs(x==numeric(1,239));
5621 using std::cout; // just for fun, another way of...
5622 using std::endl; // ...dealing with this namespace std.
5624 for (int i=2; i<12; i+=2) @{
5625 pi_frac = machin_pi(i);
5626 cout << i << ":\t" << pi_frac << endl
5627 << "\t" << pi_frac.evalf() << endl;
5633 Note how we just called @code{.series(x,degr)} instead of
5634 @code{.series(x==0,degr)}. This is a simple shortcut for @code{ex}'s
5635 method @code{series()}: if the first argument is a symbol the expression
5636 is expanded in that symbol around point @code{0}. When you run this
5637 program, it will type out:
5641 3.1832635983263598326
5642 4: 5359397032/1706489875
5643 3.1405970293260603143
5644 6: 38279241713339684/12184551018734375
5645 3.141621029325034425
5646 8: 76528487109180192540976/24359780855939418203125
5647 3.141591772182177295
5648 10: 327853873402258685803048818236/104359128170408663038552734375
5649 3.1415926824043995174
5653 @node Symmetrization, Built-in functions, Series expansion, Methods and functions
5654 @c node-name, next, previous, up
5655 @section Symmetrization
5656 @cindex @code{symmetrize()}
5657 @cindex @code{antisymmetrize()}
5658 @cindex @code{symmetrize_cyclic()}
5663 ex ex::symmetrize(const lst & l);
5664 ex ex::antisymmetrize(const lst & l);
5665 ex ex::symmetrize_cyclic(const lst & l);
5668 symmetrize an expression by returning the sum over all symmetric,
5669 antisymmetric or cyclic permutations of the specified list of objects,
5670 weighted by the number of permutations.
5672 The three additional methods
5675 ex ex::symmetrize();
5676 ex ex::antisymmetrize();
5677 ex ex::symmetrize_cyclic();
5680 symmetrize or antisymmetrize an expression over its free indices.
5682 Symmetrization is most useful with indexed expressions but can be used with
5683 almost any kind of object (anything that is @code{subs()}able):
5687 idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
5688 symbol A("A"), B("B"), a("a"), b("b"), c("c");
5690 cout << indexed(A, i, j).symmetrize() << endl;
5691 // -> 1/2*A.j.i+1/2*A.i.j
5692 cout << indexed(A, i, j, k).antisymmetrize(lst(i, j)) << endl;
5693 // -> -1/2*A.j.i.k+1/2*A.i.j.k
5694 cout << lst(a, b, c).symmetrize_cyclic(lst(a, b, c)) << endl;
5695 // -> 1/3*@{a,b,c@}+1/3*@{b,c,a@}+1/3*@{c,a,b@}
5699 @node Built-in functions, Multiple polylogarithms, Symmetrization, Methods and functions
5700 @c node-name, next, previous, up
5701 @section Predefined mathematical functions
5703 @subsection Overview
5705 GiNaC contains the following predefined mathematical functions:
5708 @multitable @columnfractions .30 .70
5709 @item @strong{Name} @tab @strong{Function}
5712 @cindex @code{abs()}
5713 @item @code{step(x)}
5715 @cindex @code{step()}
5716 @item @code{csgn(x)}
5718 @cindex @code{conjugate()}
5719 @item @code{conjugate(x)}
5720 @tab complex conjugation
5721 @cindex @code{real_part()}
5722 @item @code{real_part(x)}
5724 @cindex @code{imag_part()}
5725 @item @code{imag_part(x)}
5727 @item @code{sqrt(x)}
5728 @tab square root (not a GiNaC function, rather an alias for @code{pow(x, numeric(1, 2))})
5729 @cindex @code{sqrt()}
5732 @cindex @code{sin()}
5735 @cindex @code{cos()}
5738 @cindex @code{tan()}
5739 @item @code{asin(x)}
5741 @cindex @code{asin()}
5742 @item @code{acos(x)}
5744 @cindex @code{acos()}
5745 @item @code{atan(x)}
5746 @tab inverse tangent
5747 @cindex @code{atan()}
5748 @item @code{atan2(y, x)}
5749 @tab inverse tangent with two arguments
5750 @item @code{sinh(x)}
5751 @tab hyperbolic sine
5752 @cindex @code{sinh()}
5753 @item @code{cosh(x)}
5754 @tab hyperbolic cosine
5755 @cindex @code{cosh()}
5756 @item @code{tanh(x)}
5757 @tab hyperbolic tangent
5758 @cindex @code{tanh()}
5759 @item @code{asinh(x)}
5760 @tab inverse hyperbolic sine
5761 @cindex @code{asinh()}
5762 @item @code{acosh(x)}
5763 @tab inverse hyperbolic cosine
5764 @cindex @code{acosh()}
5765 @item @code{atanh(x)}
5766 @tab inverse hyperbolic tangent
5767 @cindex @code{atanh()}
5769 @tab exponential function
5770 @cindex @code{exp()}
5772 @tab natural logarithm
5773 @cindex @code{log()}
5776 @cindex @code{Li2()}
5777 @item @code{Li(m, x)}
5778 @tab classical polylogarithm as well as multiple polylogarithm
5780 @item @code{G(a, y)}
5781 @tab multiple polylogarithm
5783 @item @code{G(a, s, y)}
5784 @tab multiple polylogarithm with explicit signs for the imaginary parts
5786 @item @code{S(n, p, x)}
5787 @tab Nielsen's generalized polylogarithm
5789 @item @code{H(m, x)}
5790 @tab harmonic polylogarithm
5792 @item @code{zeta(m)}
5793 @tab Riemann's zeta function as well as multiple zeta value
5794 @cindex @code{zeta()}
5795 @item @code{zeta(m, s)}
5796 @tab alternating Euler sum
5797 @cindex @code{zeta()}
5798 @item @code{zetaderiv(n, x)}
5799 @tab derivatives of Riemann's zeta function
5800 @item @code{tgamma(x)}
5802 @cindex @code{tgamma()}
5803 @cindex gamma function
5804 @item @code{lgamma(x)}
5805 @tab logarithm of gamma function
5806 @cindex @code{lgamma()}
5807 @item @code{beta(x, y)}
5808 @tab beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
5809 @cindex @code{beta()}
5811 @tab psi (digamma) function
5812 @cindex @code{psi()}
5813 @item @code{psi(n, x)}
5814 @tab derivatives of psi function (polygamma functions)
5815 @item @code{factorial(n)}
5816 @tab factorial function @math{n!}
5817 @cindex @code{factorial()}
5818 @item @code{binomial(n, k)}
5819 @tab binomial coefficients
5820 @cindex @code{binomial()}
5821 @item @code{Order(x)}
5822 @tab order term function in truncated power series
5823 @cindex @code{Order()}
5828 For functions that have a branch cut in the complex plane GiNaC follows
5829 the conventions for C++ as defined in the ANSI standard as far as
5830 possible. In particular: the natural logarithm (@code{log}) and the
5831 square root (@code{sqrt}) both have their branch cuts running along the
5832 negative real axis where the points on the axis itself belong to the
5833 upper part (i.e. continuous with quadrant II). The inverse
5834 trigonometric and hyperbolic functions are not defined for complex
5835 arguments by the C++ standard, however. In GiNaC we follow the
5836 conventions used by CLN, which in turn follow the carefully designed
5837 definitions in the Common Lisp standard. It should be noted that this
5838 convention is identical to the one used by the C99 standard and by most
5839 serious CAS. It is to be expected that future revisions of the C++
5840 standard incorporate these functions in the complex domain in a manner
5841 compatible with C99.
5843 @node Multiple polylogarithms, Complex expressions, Built-in functions, Methods and functions
5844 @c node-name, next, previous, up
5845 @subsection Multiple polylogarithms
5847 @cindex polylogarithm
5848 @cindex Nielsen's generalized polylogarithm
5849 @cindex harmonic polylogarithm
5850 @cindex multiple zeta value
5851 @cindex alternating Euler sum
5852 @cindex multiple polylogarithm
5854 The multiple polylogarithm is the most generic member of a family of functions,
5855 to which others like the harmonic polylogarithm, Nielsen's generalized
5856 polylogarithm and the multiple zeta value belong.
5857 Everyone of these functions can also be written as a multiple polylogarithm with specific
5858 parameters. This whole family of functions is therefore often referred to simply as
5859 multiple polylogarithms, containing @code{Li}, @code{G}, @code{H}, @code{S} and @code{zeta}.
5860 The multiple polylogarithm itself comes in two variants: @code{Li} and @code{G}. While
5861 @code{Li} and @code{G} in principle represent the same function, the different
5862 notations are more natural to the series representation or the integral
5863 representation, respectively.
5865 To facilitate the discussion of these functions we distinguish between indices and
5866 arguments as parameters. In the table above indices are printed as @code{m}, @code{s},
5867 @code{n} or @code{p}, whereas arguments are printed as @code{x}, @code{a} and @code{y}.
5869 To define a @code{Li}, @code{H} or @code{zeta} with a depth greater than one, you have to
5870 pass a GiNaC @code{lst} for the indices @code{m} and @code{s}, and in the case of @code{Li}
5871 for the argument @code{x} as well. The parameter @code{a} of @code{G} must always be a @code{lst} containing
5872 the arguments in expanded form. If @code{G} is used with a third parameter @code{s}, @code{s} must
5873 have the same length as @code{a}. It contains then the signs of the imaginary parts of the arguments. If
5874 @code{s} is not given, the signs default to +1.
5875 Note that @code{Li} and @code{zeta} are polymorphic in this respect. They can stand in for
5876 the classical polylogarithm and Riemann's zeta function (if depth is one), as well as for
5877 the multiple polylogarithm and the multiple zeta value, respectively. Note also, that
5878 GiNaC doesn't check whether the @code{lst}s for two parameters do have the same length.
5879 It is up to the user to ensure this, otherwise evaluating will result in undefined behavior.
5881 The functions print in LaTeX format as
5883 ${\rm Li\;\!}_{m_1,m_2,\ldots,m_k}(x_1,x_2,\ldots,x_k)$,
5889 ${\rm H\;\!}_{m_1,m_2,\ldots,m_k}(x)$ and
5892 $\zeta(m_1,m_2,\ldots,m_k)$.
5894 If @code{zeta} is an alternating zeta sum, i.e. @code{zeta(m,s)}, the indices with negative sign
5895 are printed with a line above, e.g.
5897 $\zeta(5,\overline{2})$.
5899 The order of indices and arguments in the GiNaC @code{lst}s and in the output is the same.
5901 Definitions and analytical as well as numerical properties of multiple polylogarithms
5902 are too numerous to be covered here. Instead, the user is referred to the publications listed at the
5903 end of this section. The implementation in GiNaC adheres to the definitions and conventions therein,
5904 except for a few differences which will be explicitly stated in the following.
5906 One difference is about the order of the indices and arguments. For GiNaC we adopt the convention
5907 that the indices and arguments are understood to be in the same order as in which they appear in
5908 the series representation. This means
5910 ${\rm Li\;\!}_{m_1,m_2,m_3}(x,1,1) = {\rm H\;\!}_{m_1,m_2,m_3}(x)$ and
5913 ${\rm Li\;\!}_{2,1}(1,1) = \zeta(2,1) = \zeta(3)$, but
5916 $\zeta(1,2)$ evaluates to infinity.
5918 So in comparison to the referenced publications the order of indices and arguments for @code{Li}
5921 The functions only evaluate if the indices are integers greater than zero, except for the indices
5922 @code{s} in @code{zeta} and @code{G} as well as @code{m} in @code{H}. Since @code{s}
5923 will be interpreted as the sequence of signs for the corresponding indices
5924 @code{m} or the sign of the imaginary part for the
5925 corresponding arguments @code{a}, it must contain 1 or -1, e.g.
5926 @code{zeta(lst(3,4), lst(-1,1))} means
5928 $\zeta(\overline{3},4)$
5931 @code{G(lst(a,b), lst(-1,1), c)} means
5933 $G(a-0\epsilon,b+0\epsilon;c)$.
5935 The definition of @code{H} allows indices to be 0, 1 or -1 (in expanded notation) or equally to
5936 be any integer (in compact notation). With GiNaC expanded and compact notation can be mixed,
5937 e.g. @code{lst(0,0,-1,0,1,0,0)}, @code{lst(0,0,-1,2,0,0)} and @code{lst(-3,2,0,0)} are equivalent as
5938 indices. The anonymous evaluator @code{eval()} tries to reduce the functions, if possible, to
5939 the least-generic multiple polylogarithm. If all arguments are unit, it returns @code{zeta}.
5940 Arguments equal to zero get considered, too. Riemann's zeta function @code{zeta} (with depth one)
5941 evaluates also for negative integers and positive even integers. For example:
5944 > Li(@{3,1@},@{x,1@});
5947 -zeta(@{3,2@},@{-1,-1@})
5952 It is easy to tell for a given function into which other function it can be rewritten, may
5953 it be a less-generic or a more-generic one, except for harmonic polylogarithms @code{H}
5954 with negative indices or trailing zeros (the example above gives a hint). Signs can
5955 quickly be messed up, for example. Therefore GiNaC offers a C++ function
5956 @code{convert_H_to_Li()} to deal with the upgrade of a @code{H} to a multiple polylogarithm
5957 @code{Li} (@code{eval()} already cares for the possible downgrade):
5960 > convert_H_to_Li(@{0,-2,-1,3@},x);
5961 Li(@{3,1,3@},@{-x,1,-1@})
5962 > convert_H_to_Li(@{2,-1,0@},x);
5963 -Li(@{2,1@},@{x,-1@})*log(x)+2*Li(@{3,1@},@{x,-1@})+Li(@{2,2@},@{x,-1@})
5966 Every function can be numerically evaluated for
5967 arbitrary real or complex arguments. The precision is arbitrary and can be set through the
5968 global variable @code{Digits}:
5973 > evalf(zeta(@{3,1,3,1@}));
5974 0.005229569563530960100930652283899231589890420784634635522547448972148869544...
5977 Note that the convention for arguments on the branch cut in GiNaC as stated above is
5978 different from the one Remiddi and Vermaseren have chosen for the harmonic polylogarithm.
5980 If a function evaluates to infinity, no exceptions are raised, but the function is returned
5985 In long expressions this helps a lot with debugging, because you can easily spot
5986 the divergencies. But on the other hand, you have to make sure for yourself, that no illegal
5987 cancellations of divergencies happen.
5989 Useful publications:
5991 @cite{Nested Sums, Expansion of Transcendental Functions and Multi-Scale Multi-Loop Integrals},
5992 S.Moch, P.Uwer, S.Weinzierl, hep-ph/0110083
5994 @cite{Harmonic Polylogarithms},
5995 E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
5997 @cite{Special Values of Multiple Polylogarithms},
5998 J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941
6000 @cite{Numerical Evaluation of Multiple Polylogarithms},
6001 J.Vollinga, S.Weinzierl, hep-ph/0410259
6003 @node Complex expressions, Solving linear systems of equations, Multiple polylogarithms, Methods and functions
6004 @c node-name, next, previous, up
6005 @section Complex expressions
6007 @cindex @code{conjugate()}
6009 For dealing with complex expressions there are the methods
6017 that return respectively the complex conjugate, the real part and the
6018 imaginary part of an expression. Complex conjugation works as expected
6019 for all built-in functinos and objects. Taking real and imaginary
6020 parts has not yet been implemented for all built-in functions. In cases where
6021 it is not known how to conjugate or take a real/imaginary part one
6022 of the functions @code{conjugate}, @code{real_part} or @code{imag_part}
6023 is returned. For instance, in case of a complex symbol @code{x}
6024 (symbols are complex by default), one could not simplify
6025 @code{conjugate(x)}. In the case of strings of gamma matrices,
6026 the @code{conjugate} method takes the Dirac conjugate.
6031 varidx a(symbol("a"), 4), b(symbol("b"), 4);
6035 cout << (3*I*x*y + sin(2*Pi*I*y)).conjugate() << endl;
6036 // -> -3*I*conjugate(x)*y+sin(-2*I*Pi*y)
6037 cout << (dirac_gamma(a)*dirac_gamma(b)*dirac_gamma5()).conjugate() << endl;
6038 // -> -gamma5*gamma~b*gamma~a
6042 If you declare your own GiNaC functions, then they will conjugate themselves
6043 by conjugating their arguments. This is the default strategy. If you want to
6044 change this behavior, you have to supply a specialized conjugation method
6045 for your function (see @ref{Symbolic functions} and the GiNaC source-code
6046 for @code{abs} as an example). Also, specialized methods can be provided
6047 to take real and imaginary parts of user-defined functions.
6049 @node Solving linear systems of equations, Input/output, Complex expressions, Methods and functions
6050 @c node-name, next, previous, up
6051 @section Solving linear systems of equations
6052 @cindex @code{lsolve()}
6054 The function @code{lsolve()} provides a convenient wrapper around some
6055 matrix operations that comes in handy when a system of linear equations
6059 ex lsolve(const ex & eqns, const ex & symbols,
6060 unsigned options = solve_algo::automatic);
6063 Here, @code{eqns} is a @code{lst} of equalities (i.e. class
6064 @code{relational}) while @code{symbols} is a @code{lst} of
6065 indeterminates. (@xref{The class hierarchy}, for an exposition of class
6068 It returns the @code{lst} of solutions as an expression. As an example,
6069 let us solve the two equations @code{a*x+b*y==3} and @code{x-y==b}:
6073 symbol a("a"), b("b"), x("x"), y("y");
6075 eqns = a*x+b*y==3, x-y==b;
6077 cout << lsolve(eqns, vars) << endl;
6078 // -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
6081 When the linear equations @code{eqns} are underdetermined, the solution
6082 will contain one or more tautological entries like @code{x==x},
6083 depending on the rank of the system. When they are overdetermined, the
6084 solution will be an empty @code{lst}. Note the third optional parameter
6085 to @code{lsolve()}: it accepts the same parameters as
6086 @code{matrix::solve()}. This is because @code{lsolve} is just a wrapper
6090 @node Input/output, Extending GiNaC, Solving linear systems of equations, Methods and functions
6091 @c node-name, next, previous, up
6092 @section Input and output of expressions
6095 @subsection Expression output
6097 @cindex output of expressions
6099 Expressions can simply be written to any stream:
6104 ex e = 4.5*I+pow(x,2)*3/2;
6105 cout << e << endl; // prints '4.5*I+3/2*x^2'
6109 The default output format is identical to the @command{ginsh} input syntax and
6110 to that used by most computer algebra systems, but not directly pastable
6111 into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
6112 is printed as @samp{x^2}).
6114 It is possible to print expressions in a number of different formats with
6115 a set of stream manipulators;
6118 std::ostream & dflt(std::ostream & os);
6119 std::ostream & latex(std::ostream & os);
6120 std::ostream & tree(std::ostream & os);
6121 std::ostream & csrc(std::ostream & os);
6122 std::ostream & csrc_float(std::ostream & os);
6123 std::ostream & csrc_double(std::ostream & os);
6124 std::ostream & csrc_cl_N(std::ostream & os);
6125 std::ostream & index_dimensions(std::ostream & os);
6126 std::ostream & no_index_dimensions(std::ostream & os);
6129 The @code{tree}, @code{latex} and @code{csrc} formats are also available in
6130 @command{ginsh} via the @code{print()}, @code{print_latex()} and
6131 @code{print_csrc()} functions, respectively.
6134 All manipulators affect the stream state permanently. To reset the output
6135 format to the default, use the @code{dflt} manipulator:
6139 cout << latex; // all output to cout will be in LaTeX format from
6141 cout << e << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
6142 cout << sin(x/2) << endl; // prints '\sin(\frac@{1@}@{2@} x)'
6143 cout << dflt; // revert to default output format
6144 cout << e << endl; // prints '4.5*I+3/2*x^2'
6148 If you don't want to affect the format of the stream you're working with,
6149 you can output to a temporary @code{ostringstream} like this:
6154 s << latex << e; // format of cout remains unchanged
6155 cout << s.str() << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
6159 @anchor{csrc printing}
6161 @cindex @code{csrc_float}
6162 @cindex @code{csrc_double}
6163 @cindex @code{csrc_cl_N}
6164 The @code{csrc} (an alias for @code{csrc_double}), @code{csrc_float},
6165 @code{csrc_double} and @code{csrc_cl_N} manipulators set the output to a
6166 format that can be directly used in a C or C++ program. The three possible
6167 formats select the data types used for numbers (@code{csrc_cl_N} uses the
6168 classes provided by the CLN library):
6172 cout << "f = " << csrc_float << e << ";\n";
6173 cout << "d = " << csrc_double << e << ";\n";
6174 cout << "n = " << csrc_cl_N << e << ";\n";
6178 The above example will produce (note the @code{x^2} being converted to
6182 f = (3.0/2.0)*(x*x)+std::complex<float>(0.0,4.5000000e+00);
6183 d = (3.0/2.0)*(x*x)+std::complex<double>(0.0,4.5000000000000000e+00);
6184 n = cln::cl_RA("3/2")*(x*x)+cln::complex(cln::cl_I("0"),cln::cl_F("4.5_17"));
6188 The @code{tree} manipulator allows dumping the internal structure of an
6189 expression for debugging purposes:
6200 add, hash=0x0, flags=0x3, nops=2
6201 power, hash=0x0, flags=0x3, nops=2
6202 x (symbol), serial=0, hash=0xc8d5bcdd, flags=0xf
6203 2 (numeric), hash=0x6526b0fa, flags=0xf
6204 3/2 (numeric), hash=0xf9828fbd, flags=0xf
6207 4.5L0i (numeric), hash=0xa40a97e0, flags=0xf
6211 @cindex @code{latex}
6212 The @code{latex} output format is for LaTeX parsing in mathematical mode.
6213 It is rather similar to the default format but provides some braces needed
6214 by LaTeX for delimiting boxes and also converts some common objects to
6215 conventional LaTeX names. It is possible to give symbols a special name for
6216 LaTeX output by supplying it as a second argument to the @code{symbol}
6219 For example, the code snippet
6223 symbol x("x", "\\circ");
6224 ex e = lgamma(x).series(x==0,3);
6225 cout << latex << e << endl;
6232 @{(-\ln(\circ))@}+@{(-\gamma_E)@} \circ+@{(\frac@{1@}@{12@} \pi^@{2@})@} \circ^@{2@}
6233 +\mathcal@{O@}(\circ^@{3@})
6236 @cindex @code{index_dimensions}
6237 @cindex @code{no_index_dimensions}
6238 Index dimensions are normally hidden in the output. To make them visible, use
6239 the @code{index_dimensions} manipulator. The dimensions will be written in
6240 square brackets behind each index value in the default and LaTeX output
6245 symbol x("x"), y("y");
6246 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
6247 ex e = indexed(x, mu) * indexed(y, nu);
6250 // prints 'x~mu*y~nu'
6251 cout << index_dimensions << e << endl;
6252 // prints 'x~mu[4]*y~nu[4]'
6253 cout << no_index_dimensions << e << endl;
6254 // prints 'x~mu*y~nu'
6259 @cindex Tree traversal
6260 If you need any fancy special output format, e.g. for interfacing GiNaC
6261 with other algebra systems or for producing code for different
6262 programming languages, you can always traverse the expression tree yourself:
6265 static void my_print(const ex & e)
6267 if (is_a<function>(e))
6268 cout << ex_to<function>(e).get_name();
6270 cout << ex_to<basic>(e).class_name();
6272 size_t n = e.nops();
6274 for (size_t i=0; i<n; i++) @{
6286 my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
6294 add(power(numeric(3),symbol(x)),mul(sin(mul(power(constant(Pi),numeric(-1)),
6295 symbol(y))),numeric(-2)))
6298 If you need an output format that makes it possible to accurately
6299 reconstruct an expression by feeding the output to a suitable parser or
6300 object factory, you should consider storing the expression in an
6301 @code{archive} object and reading the object properties from there.
6302 See the section on archiving for more information.
6305 @subsection Expression input
6306 @cindex input of expressions
6308 GiNaC provides no way to directly read an expression from a stream because
6309 you will usually want the user to be able to enter something like @samp{2*x+sin(y)}
6310 and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
6311 @code{y} you defined in your program and there is no way to specify the
6312 desired symbols to the @code{>>} stream input operator.
6314 Instead, GiNaC lets you construct an expression from a string, specifying the
6315 list of symbols to be used:
6319 symbol x("x"), y("y");
6320 ex e("2*x+sin(y)", lst(x, y));
6324 The input syntax is the same as that used by @command{ginsh} and the stream
6325 output operator @code{<<}. The symbols in the string are matched by name to
6326 the symbols in the list and if GiNaC encounters a symbol not specified in
6327 the list it will throw an exception.
6329 With this constructor, it's also easy to implement interactive GiNaC programs:
6334 #include <stdexcept>
6335 #include <ginac/ginac.h>
6336 using namespace std;
6337 using namespace GiNaC;
6344 cout << "Enter an expression containing 'x': ";
6349 cout << "The derivative of " << e << " with respect to x is ";
6350 cout << e.diff(x) << ".\n";
6351 @} catch (exception &p) @{
6352 cerr << p.what() << endl;
6357 @subsection Compiling expressions to C function pointers
6358 @cindex compiling expressions
6360 Numerical evaluation of algebraic expressions is seamlessly integrated into
6361 GiNaC by help of the CLN library. While CLN allows for very fast arbitrary
6362 precision numerics, which is more than sufficient for most users, sometimes only
6363 the speed of built-in floating point numbers is fast enough, e.g. for Monte
6364 Carlo integration. The only viable option then is the following: print the
6365 expression in C syntax format, manually add necessary C code, compile that
6366 program and run is as a separate application. This is not only cumbersome and
6367 involves a lot of manual intervention, but it also separates the algebraic and
6368 the numerical evaluation into different execution stages.
6370 GiNaC offers a couple of functions that help to avoid these inconveniences and
6371 problems. The functions automatically perform the printing of a GiNaC expression
6372 and the subsequent compiling of its associated C code. The created object code
6373 is then dynamically linked to the currently running program. A function pointer
6374 to the C function that performs the numerical evaluation is returned and can be
6375 used instantly. This all happens automatically, no user intervention is needed.
6377 The following example demonstrates the use of @code{compile_ex}:
6382 ex myexpr = sin(x) / x;
6385 compile_ex(myexpr, x, fp);
6387 cout << fp(3.2) << endl;
6391 The function @code{compile_ex} is called with the expression to be compiled and
6392 its only free variable @code{x}. Upon successful completion the third parameter
6393 contains a valid function pointer to the corresponding C code module. If called
6394 like in the last line only built-in double precision numerics is involved.
6399 The function pointer has to be defined in advance. GiNaC offers three function
6400 pointer types at the moment:
6403 typedef double (*FUNCP_1P) (double);
6404 typedef double (*FUNCP_2P) (double, double);
6405 typedef void (*FUNCP_CUBA) (const int*, const double[], const int*, double[]);
6408 @cindex CUBA library
6409 @cindex Monte Carlo integration
6410 @code{FUNCP_2P} allows for two variables in the expression. @code{FUNCP_CUBA} is
6411 the correct type to be used with the CUBA library
6412 (@uref{http://www.feynarts/cuba}) for numerical integrations. The details for the
6413 parameters of @code{FUNCP_CUBA} are explained in the CUBA manual.
6416 For every function pointer type there is a matching @code{compile_ex} available:
6419 void compile_ex(const ex& expr, const symbol& sym, FUNCP_1P& fp,
6420 const std::string filename = "");
6421 void compile_ex(const ex& expr, const symbol& sym1, const symbol& sym2,
6422 FUNCP_2P& fp, const std::string filename = "");
6423 void compile_ex(const lst& exprs, const lst& syms, FUNCP_CUBA& fp,
6424 const std::string filename = "");
6427 When the last parameter @code{filename} is not supplied, @code{compile_ex} will
6428 choose a unique random name for the intermediate source and object files it
6429 produces. On program termination these files will be deleted. If one wishes to
6430 keep the C code and the object files, one can supply the @code{filename}
6431 parameter. The intermediate files will use that filename and will not be
6435 @code{link_ex} is a function that allows to dynamically link an existing object
6436 file and to make it available via a function pointer. This is useful if you
6437 have already used @code{compile_ex} on an expression and want to avoid the
6438 compilation step to be performed over and over again when you restart your
6439 program. The precondition for this is of course, that you have chosen a
6440 filename when you did call @code{compile_ex}. For every above mentioned
6441 function pointer type there exists a corresponding @code{link_ex} function:
6444 void link_ex(const std::string filename, FUNCP_1P& fp);
6445 void link_ex(const std::string filename, FUNCP_2P& fp);
6446 void link_ex(const std::string filename, FUNCP_CUBA& fp);
6449 The complete filename (including the suffix @code{.so}) of the object file has
6456 void unlink_ex(const std::string filename);
6459 is supplied for the rare cases when one wishes to close the dynamically linked
6460 object files directly and have the intermediate files (only if filename has not
6461 been given) deleted. Normally one doesn't need this function, because all the
6462 clean-up will be done automatically upon (regular) program termination.
6464 All the described functions will throw an exception in case they cannot perform
6465 correctly, like for example when writing the file or starting the compiler
6466 fails. Since internally the same printing methods as described in section
6467 @ref{csrc printing} are used, only functions and objects that are available in
6468 standard C will compile successfully (that excludes polylogarithms for example
6469 at the moment). Another precondition for success is, of course, that it must be
6470 possible to evaluate the expression numerically. No free variables despite the
6471 ones supplied to @code{compile_ex} should appear in the expression.
6473 @cindex ginac-excompiler
6474 @code{compile_ex} uses the shell script @code{ginac-excompiler} to start the C
6475 compiler and produce the object files. This shell script comes with GiNaC and
6476 will be installed together with GiNaC in the configured @code{$PREFIX/bin}
6479 @subsection Archiving
6480 @cindex @code{archive} (class)
6483 GiNaC allows creating @dfn{archives} of expressions which can be stored
6484 to or retrieved from files. To create an archive, you declare an object
6485 of class @code{archive} and archive expressions in it, giving each
6486 expression a unique name:
6490 using namespace std;
6491 #include <ginac/ginac.h>
6492 using namespace GiNaC;
6496 symbol x("x"), y("y"), z("z");
6498 ex foo = sin(x + 2*y) + 3*z + 41;
6502 a.archive_ex(foo, "foo");
6503 a.archive_ex(bar, "the second one");
6507 The archive can then be written to a file:
6511 ofstream out("foobar.gar");
6517 The file @file{foobar.gar} contains all information that is needed to
6518 reconstruct the expressions @code{foo} and @code{bar}.
6520 @cindex @command{viewgar}
6521 The tool @command{viewgar} that comes with GiNaC can be used to view
6522 the contents of GiNaC archive files:
6525 $ viewgar foobar.gar
6526 foo = 41+sin(x+2*y)+3*z
6527 the second one = 42+sin(x+2*y)+3*z
6530 The point of writing archive files is of course that they can later be
6536 ifstream in("foobar.gar");
6541 And the stored expressions can be retrieved by their name:
6548 ex ex1 = a2.unarchive_ex(syms, "foo");
6549 ex ex2 = a2.unarchive_ex(syms, "the second one");
6551 cout << ex1 << endl; // prints "41+sin(x+2*y)+3*z"
6552 cout << ex2 << endl; // prints "42+sin(x+2*y)+3*z"
6553 cout << ex1.subs(x == 2) << endl; // prints "41+sin(2+2*y)+3*z"
6557 Note that you have to supply a list of the symbols which are to be inserted
6558 in the expressions. Symbols in archives are stored by their name only and
6559 if you don't specify which symbols you have, unarchiving the expression will
6560 create new symbols with that name. E.g. if you hadn't included @code{x} in
6561 the @code{syms} list above, the @code{ex1.subs(x == 2)} statement would
6562 have had no effect because the @code{x} in @code{ex1} would have been a
6563 different symbol than the @code{x} which was defined at the beginning of
6564 the program, although both would appear as @samp{x} when printed.
6566 You can also use the information stored in an @code{archive} object to
6567 output expressions in a format suitable for exact reconstruction. The
6568 @code{archive} and @code{archive_node} classes have a couple of member
6569 functions that let you access the stored properties:
6572 static void my_print2(const archive_node & n)
6575 n.find_string("class", class_name);
6576 cout << class_name << "(";
6578 archive_node::propinfovector p;
6579 n.get_properties(p);
6581 size_t num = p.size();
6582 for (size_t i=0; i<num; i++) @{
6583 const string &name = p[i].name;
6584 if (name == "class")
6586 cout << name << "=";
6588 unsigned count = p[i].count;
6592 for (unsigned j=0; j<count; j++) @{
6593 switch (p[i].type) @{
6594 case archive_node::PTYPE_BOOL: @{
6596 n.find_bool(name, x, j);
6597 cout << (x ? "true" : "false");
6600 case archive_node::PTYPE_UNSIGNED: @{
6602 n.find_unsigned(name, x, j);
6606 case archive_node::PTYPE_STRING: @{
6608 n.find_string(name, x, j);
6609 cout << '\"' << x << '\"';
6612 case archive_node::PTYPE_NODE: @{
6613 const archive_node &x = n.find_ex_node(name, j);
6635 ex e = pow(2, x) - y;
6637 my_print2(ar.get_top_node(0)); cout << endl;
6645 add(rest=@{power(basis=numeric(number="2"),exponent=symbol(name="x")),
6646 symbol(name="y")@},coeff=@{numeric(number="1"),numeric(number="-1")@},
6647 overall_coeff=numeric(number="0"))
6650 Be warned, however, that the set of properties and their meaning for each
6651 class may change between GiNaC versions.
6654 @node Extending GiNaC, What does not belong into GiNaC, Input/output, Top
6655 @c node-name, next, previous, up
6656 @chapter Extending GiNaC
6658 By reading so far you should have gotten a fairly good understanding of
6659 GiNaC's design patterns. From here on you should start reading the
6660 sources. All we can do now is issue some recommendations how to tackle
6661 GiNaC's many loose ends in order to fulfill everybody's dreams. If you
6662 develop some useful extension please don't hesitate to contact the GiNaC
6663 authors---they will happily incorporate them into future versions.
6666 * What does not belong into GiNaC:: What to avoid.
6667 * Symbolic functions:: Implementing symbolic functions.
6668 * Printing:: Adding new output formats.
6669 * Structures:: Defining new algebraic classes (the easy way).
6670 * Adding classes:: Defining new algebraic classes (the hard way).
6674 @node What does not belong into GiNaC, Symbolic functions, Extending GiNaC, Extending GiNaC
6675 @c node-name, next, previous, up
6676 @section What doesn't belong into GiNaC
6678 @cindex @command{ginsh}
6679 First of all, GiNaC's name must be read literally. It is designed to be
6680 a library for use within C++. The tiny @command{ginsh} accompanying
6681 GiNaC makes this even more clear: it doesn't even attempt to provide a
6682 language. There are no loops or conditional expressions in
6683 @command{ginsh}, it is merely a window into the library for the
6684 programmer to test stuff (or to show off). Still, the design of a
6685 complete CAS with a language of its own, graphical capabilities and all
6686 this on top of GiNaC is possible and is without doubt a nice project for
6689 There are many built-in functions in GiNaC that do not know how to
6690 evaluate themselves numerically to a precision declared at runtime
6691 (using @code{Digits}). Some may be evaluated at certain points, but not
6692 generally. This ought to be fixed. However, doing numerical
6693 computations with GiNaC's quite abstract classes is doomed to be
6694 inefficient. For this purpose, the underlying foundation classes
6695 provided by CLN are much better suited.
6698 @node Symbolic functions, Printing, What does not belong into GiNaC, Extending GiNaC
6699 @c node-name, next, previous, up
6700 @section Symbolic functions
6702 The easiest and most instructive way to start extending GiNaC is probably to
6703 create your own symbolic functions. These are implemented with the help of
6704 two preprocessor macros:
6706 @cindex @code{DECLARE_FUNCTION}
6707 @cindex @code{REGISTER_FUNCTION}
6709 DECLARE_FUNCTION_<n>P(<name>)
6710 REGISTER_FUNCTION(<name>, <options>)
6713 The @code{DECLARE_FUNCTION} macro will usually appear in a header file. It
6714 declares a C++ function with the given @samp{name} that takes exactly @samp{n}
6715 parameters of type @code{ex} and returns a newly constructed GiNaC
6716 @code{function} object that represents your function.
6718 The @code{REGISTER_FUNCTION} macro implements the function. It must be passed
6719 the same @samp{name} as the respective @code{DECLARE_FUNCTION} macro, and a
6720 set of options that associate the symbolic function with C++ functions you
6721 provide to implement the various methods such as evaluation, derivative,
6722 series expansion etc. They also describe additional attributes the function
6723 might have, such as symmetry and commutation properties, and a name for
6724 LaTeX output. Multiple options are separated by the member access operator
6725 @samp{.} and can be given in an arbitrary order.
6727 (By the way: in case you are worrying about all the macros above we can
6728 assure you that functions are GiNaC's most macro-intense classes. We have
6729 done our best to avoid macros where we can.)
6731 @subsection A minimal example
6733 Here is an example for the implementation of a function with two arguments
6734 that is not further evaluated:
6737 DECLARE_FUNCTION_2P(myfcn)
6739 REGISTER_FUNCTION(myfcn, dummy())
6742 Any code that has seen the @code{DECLARE_FUNCTION} line can use @code{myfcn()}
6743 in algebraic expressions:
6749 ex e = 2*myfcn(42, 1+3*x) - x;
6751 // prints '2*myfcn(42,1+3*x)-x'
6756 The @code{dummy()} option in the @code{REGISTER_FUNCTION} line signifies
6757 "no options". A function with no options specified merely acts as a kind of
6758 container for its arguments. It is a pure "dummy" function with no associated
6759 logic (which is, however, sometimes perfectly sufficient).
6761 Let's now have a look at the implementation of GiNaC's cosine function for an
6762 example of how to make an "intelligent" function.
6764 @subsection The cosine function
6766 The GiNaC header file @file{inifcns.h} contains the line
6769 DECLARE_FUNCTION_1P(cos)
6772 which declares to all programs using GiNaC that there is a function @samp{cos}
6773 that takes one @code{ex} as an argument. This is all they need to know to use
6774 this function in expressions.
6776 The implementation of the cosine function is in @file{inifcns_trans.cpp}. Here
6777 is its @code{REGISTER_FUNCTION} line:
6780 REGISTER_FUNCTION(cos, eval_func(cos_eval).
6781 evalf_func(cos_evalf).
6782 derivative_func(cos_deriv).
6783 latex_name("\\cos"));
6786 There are four options defined for the cosine function. One of them
6787 (@code{latex_name}) gives the function a proper name for LaTeX output; the
6788 other three indicate the C++ functions in which the "brains" of the cosine
6789 function are defined.
6791 @cindex @code{hold()}
6793 The @code{eval_func()} option specifies the C++ function that implements
6794 the @code{eval()} method, GiNaC's anonymous evaluator. This function takes
6795 the same number of arguments as the associated symbolic function (one in this
6796 case) and returns the (possibly transformed or in some way simplified)
6797 symbolically evaluated function (@xref{Automatic evaluation}, for a description
6798 of the automatic evaluation process). If no (further) evaluation is to take
6799 place, the @code{eval_func()} function must return the original function
6800 with @code{.hold()}, to avoid a potential infinite recursion. If your
6801 symbolic functions produce a segmentation fault or stack overflow when
6802 using them in expressions, you are probably missing a @code{.hold()}
6805 The @code{eval_func()} function for the cosine looks something like this
6806 (actually, it doesn't look like this at all, but it should give you an idea
6810 static ex cos_eval(const ex & x)
6812 if ("x is a multiple of 2*Pi")
6814 else if ("x is a multiple of Pi")
6816 else if ("x is a multiple of Pi/2")
6820 else if ("x has the form 'acos(y)'")
6822 else if ("x has the form 'asin(y)'")
6827 return cos(x).hold();
6831 This function is called every time the cosine is used in a symbolic expression:
6837 // this calls cos_eval(Pi), and inserts its return value into
6838 // the actual expression
6845 In this way, @code{cos(4*Pi)} automatically becomes @math{1},
6846 @code{cos(asin(a+b))} becomes @code{sqrt(1-(a+b)^2)}, etc. If no reasonable
6847 symbolic transformation can be done, the unmodified function is returned
6848 with @code{.hold()}.
6850 GiNaC doesn't automatically transform @code{cos(2)} to @samp{-0.416146...}.
6851 The user has to call @code{evalf()} for that. This is implemented in a
6855 static ex cos_evalf(const ex & x)
6857 if (is_a<numeric>(x))
6858 return cos(ex_to<numeric>(x));
6860 return cos(x).hold();
6864 Since we are lazy we defer the problem of numeric evaluation to somebody else,
6865 in this case the @code{cos()} function for @code{numeric} objects, which in
6866 turn hands it over to the @code{cos()} function in CLN. The @code{.hold()}
6867 isn't really needed here, but reminds us that the corresponding @code{eval()}
6868 function would require it in this place.
6870 Differentiation will surely turn up and so we need to tell @code{cos}
6871 what its first derivative is (higher derivatives, @code{.diff(x,3)} for
6872 instance, are then handled automatically by @code{basic::diff} and
6876 static ex cos_deriv(const ex & x, unsigned diff_param)
6882 @cindex product rule
6883 The second parameter is obligatory but uninteresting at this point. It
6884 specifies which parameter to differentiate in a partial derivative in
6885 case the function has more than one parameter, and its main application
6886 is for correct handling of the chain rule.
6888 An implementation of the series expansion is not needed for @code{cos()} as
6889 it doesn't have any poles and GiNaC can do Taylor expansion by itself (as
6890 long as it knows what the derivative of @code{cos()} is). @code{tan()}, on
6891 the other hand, does have poles and may need to do Laurent expansion:
6894 static ex tan_series(const ex & x, const relational & rel,
6895 int order, unsigned options)
6897 // Find the actual expansion point
6898 const ex x_pt = x.subs(rel);
6900 if ("x_pt is not an odd multiple of Pi/2")
6901 throw do_taylor(); // tell function::series() to do Taylor expansion
6903 // On a pole, expand sin()/cos()
6904 return (sin(x)/cos(x)).series(rel, order+2, options);
6908 The @code{series()} implementation of a function @emph{must} return a
6909 @code{pseries} object, otherwise your code will crash.
6911 @subsection Function options
6913 GiNaC functions understand several more options which are always
6914 specified as @code{.option(params)}. None of them are required, but you
6915 need to specify at least one option to @code{REGISTER_FUNCTION()}. There
6916 is a do-nothing option called @code{dummy()} which you can use to define
6917 functions without any special options.
6920 eval_func(<C++ function>)
6921 evalf_func(<C++ function>)
6922 derivative_func(<C++ function>)
6923 series_func(<C++ function>)
6924 conjugate_func(<C++ function>)
6927 These specify the C++ functions that implement symbolic evaluation,
6928 numeric evaluation, partial derivatives, and series expansion, respectively.
6929 They correspond to the GiNaC methods @code{eval()}, @code{evalf()},
6930 @code{diff()} and @code{series()}.
6932 The @code{eval_func()} function needs to use @code{.hold()} if no further
6933 automatic evaluation is desired or possible.
6935 If no @code{series_func()} is given, GiNaC defaults to simple Taylor
6936 expansion, which is correct if there are no poles involved. If the function
6937 has poles in the complex plane, the @code{series_func()} needs to check
6938 whether the expansion point is on a pole and fall back to Taylor expansion
6939 if it isn't. Otherwise, the pole usually needs to be regularized by some
6940 suitable transformation.
6943 latex_name(const string & n)
6946 specifies the LaTeX code that represents the name of the function in LaTeX
6947 output. The default is to put the function name in an @code{\mbox@{@}}.
6950 do_not_evalf_params()
6953 This tells @code{evalf()} to not recursively evaluate the parameters of the
6954 function before calling the @code{evalf_func()}.
6957 set_return_type(unsigned return_type, unsigned return_type_tinfo)
6960 This allows you to explicitly specify the commutation properties of the
6961 function (@xref{Non-commutative objects}, for an explanation of
6962 (non)commutativity in GiNaC). For example, you can use
6963 @code{set_return_type(return_types::noncommutative, TINFO_matrix)} to make
6964 GiNaC treat your function like a matrix. By default, functions inherit the
6965 commutation properties of their first argument.
6968 set_symmetry(const symmetry & s)
6971 specifies the symmetry properties of the function with respect to its
6972 arguments. @xref{Indexed objects}, for an explanation of symmetry
6973 specifications. GiNaC will automatically rearrange the arguments of
6974 symmetric functions into a canonical order.
6976 Sometimes you may want to have finer control over how functions are
6977 displayed in the output. For example, the @code{abs()} function prints
6978 itself as @samp{abs(x)} in the default output format, but as @samp{|x|}
6979 in LaTeX mode, and @code{fabs(x)} in C source output. This is achieved
6983 print_func<C>(<C++ function>)
6986 option which is explained in the next section.
6988 @subsection Functions with a variable number of arguments
6990 The @code{DECLARE_FUNCTION} and @code{REGISTER_FUNCTION} macros define
6991 functions with a fixed number of arguments. Sometimes, though, you may need
6992 to have a function that accepts a variable number of expressions. One way to
6993 accomplish this is to pass variable-length lists as arguments. The
6994 @code{Li()} function uses this method for multiple polylogarithms.
6996 It is also possible to define functions that accept a different number of
6997 parameters under the same function name, such as the @code{psi()} function
6998 which can be called either as @code{psi(z)} (the digamma function) or as
6999 @code{psi(n, z)} (polygamma functions). These are actually two different
7000 functions in GiNaC that, however, have the same name. Defining such
7001 functions is not possible with the macros but requires manually fiddling
7002 with GiNaC internals. If you are interested, please consult the GiNaC source
7003 code for the @code{psi()} function (@file{inifcns.h} and
7004 @file{inifcns_gamma.cpp}).
7007 @node Printing, Structures, Symbolic functions, Extending GiNaC
7008 @c node-name, next, previous, up
7009 @section GiNaC's expression output system
7011 GiNaC allows the output of expressions in a variety of different formats
7012 (@pxref{Input/output}). This section will explain how expression output
7013 is implemented internally, and how to define your own output formats or
7014 change the output format of built-in algebraic objects. You will also want
7015 to read this section if you plan to write your own algebraic classes or
7018 @cindex @code{print_context} (class)
7019 @cindex @code{print_dflt} (class)
7020 @cindex @code{print_latex} (class)
7021 @cindex @code{print_tree} (class)
7022 @cindex @code{print_csrc} (class)
7023 All the different output formats are represented by a hierarchy of classes
7024 rooted in the @code{print_context} class, defined in the @file{print.h}
7029 the default output format
7031 output in LaTeX mathematical mode
7033 a dump of the internal expression structure (for debugging)
7035 the base class for C source output
7036 @item print_csrc_float
7037 C source output using the @code{float} type
7038 @item print_csrc_double
7039 C source output using the @code{double} type
7040 @item print_csrc_cl_N
7041 C source output using CLN types
7044 The @code{print_context} base class provides two public data members:
7056 @code{s} is a reference to the stream to output to, while @code{options}
7057 holds flags and modifiers. Currently, there is only one flag defined:
7058 @code{print_options::print_index_dimensions} instructs the @code{idx} class
7059 to print the index dimension which is normally hidden.
7061 When you write something like @code{std::cout << e}, where @code{e} is
7062 an object of class @code{ex}, GiNaC will construct an appropriate
7063 @code{print_context} object (of a class depending on the selected output
7064 format), fill in the @code{s} and @code{options} members, and call
7066 @cindex @code{print()}
7068 void ex::print(const print_context & c, unsigned level = 0) const;
7071 which in turn forwards the call to the @code{print()} method of the
7072 top-level algebraic object contained in the expression.
7074 Unlike other methods, GiNaC classes don't usually override their
7075 @code{print()} method to implement expression output. Instead, the default
7076 implementation @code{basic::print(c, level)} performs a run-time double
7077 dispatch to a function selected by the dynamic type of the object and the
7078 passed @code{print_context}. To this end, GiNaC maintains a separate method
7079 table for each class, similar to the virtual function table used for ordinary
7080 (single) virtual function dispatch.
7082 The method table contains one slot for each possible @code{print_context}
7083 type, indexed by the (internally assigned) serial number of the type. Slots
7084 may be empty, in which case GiNaC will retry the method lookup with the
7085 @code{print_context} object's parent class, possibly repeating the process
7086 until it reaches the @code{print_context} base class. If there's still no
7087 method defined, the method table of the algebraic object's parent class
7088 is consulted, and so on, until a matching method is found (eventually it
7089 will reach the combination @code{basic/print_context}, which prints the
7090 object's class name enclosed in square brackets).
7092 You can think of the print methods of all the different classes and output
7093 formats as being arranged in a two-dimensional matrix with one axis listing
7094 the algebraic classes and the other axis listing the @code{print_context}
7097 Subclasses of @code{basic} can, of course, also overload @code{basic::print()}
7098 to implement printing, but then they won't get any of the benefits of the
7099 double dispatch mechanism (such as the ability for derived classes to
7100 inherit only certain print methods from its parent, or the replacement of
7101 methods at run-time).
7103 @subsection Print methods for classes
7105 The method table for a class is set up either in the definition of the class,
7106 by passing the appropriate @code{print_func<C>()} option to
7107 @code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT()} (@xref{Adding classes}, for
7108 an example), or at run-time using @code{set_print_func<T, C>()}. The latter
7109 can also be used to override existing methods dynamically.
7111 The argument to @code{print_func<C>()} and @code{set_print_func<T, C>()} can
7112 be a member function of the class (or one of its parent classes), a static
7113 member function, or an ordinary (global) C++ function. The @code{C} template
7114 parameter specifies the appropriate @code{print_context} type for which the
7115 method should be invoked, while, in the case of @code{set_print_func<>()}, the
7116 @code{T} parameter specifies the algebraic class (for @code{print_func<>()},
7117 the class is the one being implemented by
7118 @code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT}).
7120 For print methods that are member functions, their first argument must be of
7121 a type convertible to a @code{const C &}, and the second argument must be an
7124 For static members and global functions, the first argument must be of a type
7125 convertible to a @code{const T &}, the second argument must be of a type
7126 convertible to a @code{const C &}, and the third argument must be an
7127 @code{unsigned}. A global function will, of course, not have access to
7128 private and protected members of @code{T}.
7130 The @code{unsigned} argument of the print methods (and of @code{ex::print()}
7131 and @code{basic::print()}) is used for proper parenthesizing of the output
7132 (and by @code{print_tree} for proper indentation). It can be used for similar
7133 purposes if you write your own output formats.
7135 The explanations given above may seem complicated, but in practice it's
7136 really simple, as shown in the following example. Suppose that we want to
7137 display exponents in LaTeX output not as superscripts but with little
7138 upwards-pointing arrows. This can be achieved in the following way:
7141 void my_print_power_as_latex(const power & p,
7142 const print_latex & c,
7145 // get the precedence of the 'power' class
7146 unsigned power_prec = p.precedence();
7148 // if the parent operator has the same or a higher precedence
7149 // we need parentheses around the power
7150 if (level >= power_prec)
7153 // print the basis and exponent, each enclosed in braces, and
7154 // separated by an uparrow
7156 p.op(0).print(c, power_prec);
7157 c.s << "@}\\uparrow@{";
7158 p.op(1).print(c, power_prec);
7161 // don't forget the closing parenthesis
7162 if (level >= power_prec)
7168 // a sample expression
7169 symbol x("x"), y("y");
7170 ex e = -3*pow(x, 3)*pow(y, -2) + pow(x+y, 2) - 1;
7172 // switch to LaTeX mode
7175 // this prints "-1+@{(y+x)@}^@{2@}-3 \frac@{x^@{3@}@}@{y^@{2@}@}"
7178 // now we replace the method for the LaTeX output of powers with
7180 set_print_func<power, print_latex>(my_print_power_as_latex);
7182 // this prints "-1+@{@{(y+x)@}@}\uparrow@{2@}-3 \frac@{@{x@}\uparrow@{3@}@}@{@{y@}
7193 The first argument of @code{my_print_power_as_latex} could also have been
7194 a @code{const basic &}, the second one a @code{const print_context &}.
7197 The above code depends on @code{mul} objects converting their operands to
7198 @code{power} objects for the purpose of printing.
7201 The output of products including negative powers as fractions is also
7202 controlled by the @code{mul} class.
7205 The @code{power/print_latex} method provided by GiNaC prints square roots
7206 using @code{\sqrt}, but the above code doesn't.
7210 It's not possible to restore a method table entry to its previous or default
7211 value. Once you have called @code{set_print_func()}, you can only override
7212 it with another call to @code{set_print_func()}, but you can't easily go back
7213 to the default behavior again (you can, of course, dig around in the GiNaC
7214 sources, find the method that is installed at startup
7215 (@code{power::do_print_latex} in this case), and @code{set_print_func} that
7216 one; that is, after you circumvent the C++ member access control@dots{}).
7218 @subsection Print methods for functions
7220 Symbolic functions employ a print method dispatch mechanism similar to the
7221 one used for classes. The methods are specified with @code{print_func<C>()}
7222 function options. If you don't specify any special print methods, the function
7223 will be printed with its name (or LaTeX name, if supplied), followed by a
7224 comma-separated list of arguments enclosed in parentheses.
7226 For example, this is what GiNaC's @samp{abs()} function is defined like:
7229 static ex abs_eval(const ex & arg) @{ ... @}
7230 static ex abs_evalf(const ex & arg) @{ ... @}
7232 static void abs_print_latex(const ex & arg, const print_context & c)
7234 c.s << "@{|"; arg.print(c); c.s << "|@}";
7237 static void abs_print_csrc_float(const ex & arg, const print_context & c)
7239 c.s << "fabs("; arg.print(c); c.s << ")";
7242 REGISTER_FUNCTION(abs, eval_func(abs_eval).
7243 evalf_func(abs_evalf).
7244 print_func<print_latex>(abs_print_latex).
7245 print_func<print_csrc_float>(abs_print_csrc_float).
7246 print_func<print_csrc_double>(abs_print_csrc_float));
7249 This will display @samp{abs(x)} as @samp{|x|} in LaTeX mode and @code{fabs(x)}
7250 in non-CLN C source output, but as @code{abs(x)} in all other formats.
7252 There is currently no equivalent of @code{set_print_func()} for functions.
7254 @subsection Adding new output formats
7256 Creating a new output format involves subclassing @code{print_context},
7257 which is somewhat similar to adding a new algebraic class
7258 (@pxref{Adding classes}). There is a macro @code{GINAC_DECLARE_PRINT_CONTEXT}
7259 that needs to go into the class definition, and a corresponding macro
7260 @code{GINAC_IMPLEMENT_PRINT_CONTEXT} that has to appear at global scope.
7261 Every @code{print_context} class needs to provide a default constructor
7262 and a constructor from an @code{std::ostream} and an @code{unsigned}
7265 Here is an example for a user-defined @code{print_context} class:
7268 class print_myformat : public print_dflt
7270 GINAC_DECLARE_PRINT_CONTEXT(print_myformat, print_dflt)
7272 print_myformat(std::ostream & os, unsigned opt = 0)
7273 : print_dflt(os, opt) @{@}
7276 print_myformat::print_myformat() : print_dflt(std::cout) @{@}
7278 GINAC_IMPLEMENT_PRINT_CONTEXT(print_myformat, print_dflt)
7281 That's all there is to it. None of the actual expression output logic is
7282 implemented in this class. It merely serves as a selector for choosing
7283 a particular format. The algorithms for printing expressions in the new
7284 format are implemented as print methods, as described above.
7286 @code{print_myformat} is a subclass of @code{print_dflt}, so it behaves
7287 exactly like GiNaC's default output format:
7292 ex e = pow(x, 2) + 1;
7294 // this prints "1+x^2"
7297 // this also prints "1+x^2"
7298 e.print(print_myformat()); cout << endl;
7304 To fill @code{print_myformat} with life, we need to supply appropriate
7305 print methods with @code{set_print_func()}, like this:
7308 // This prints powers with '**' instead of '^'. See the LaTeX output
7309 // example above for explanations.
7310 void print_power_as_myformat(const power & p,
7311 const print_myformat & c,
7314 unsigned power_prec = p.precedence();
7315 if (level >= power_prec)
7317 p.op(0).print(c, power_prec);
7319 p.op(1).print(c, power_prec);
7320 if (level >= power_prec)
7326 // install a new print method for power objects
7327 set_print_func<power, print_myformat>(print_power_as_myformat);
7329 // now this prints "1+x**2"
7330 e.print(print_myformat()); cout << endl;
7332 // but the default format is still "1+x^2"
7338 @node Structures, Adding classes, Printing, Extending GiNaC
7339 @c node-name, next, previous, up
7342 If you are doing some very specialized things with GiNaC, or if you just
7343 need some more organized way to store data in your expressions instead of
7344 anonymous lists, you may want to implement your own algebraic classes.
7345 ('algebraic class' means any class directly or indirectly derived from
7346 @code{basic} that can be used in GiNaC expressions).
7348 GiNaC offers two ways of accomplishing this: either by using the
7349 @code{structure<T>} template class, or by rolling your own class from
7350 scratch. This section will discuss the @code{structure<T>} template which
7351 is easier to use but more limited, while the implementation of custom
7352 GiNaC classes is the topic of the next section. However, you may want to
7353 read both sections because many common concepts and member functions are
7354 shared by both concepts, and it will also allow you to decide which approach
7355 is most suited to your needs.
7357 The @code{structure<T>} template, defined in the GiNaC header file
7358 @file{structure.h}, wraps a type that you supply (usually a C++ @code{struct}
7359 or @code{class}) into a GiNaC object that can be used in expressions.
7361 @subsection Example: scalar products
7363 Let's suppose that we need a way to handle some kind of abstract scalar
7364 product of the form @samp{<x|y>} in expressions. Objects of the scalar
7365 product class have to store their left and right operands, which can in turn
7366 be arbitrary expressions. Here is a possible way to represent such a
7367 product in a C++ @code{struct}:
7371 using namespace std;
7373 #include <ginac/ginac.h>
7374 using namespace GiNaC;
7380 sprod_s(ex l, ex r) : left(l), right(r) @{@}
7384 The default constructor is required. Now, to make a GiNaC class out of this
7385 data structure, we need only one line:
7388 typedef structure<sprod_s> sprod;
7391 That's it. This line constructs an algebraic class @code{sprod} which
7392 contains objects of type @code{sprod_s}. We can now use @code{sprod} in
7393 expressions like any other GiNaC class:
7397 symbol a("a"), b("b");
7398 ex e = sprod(sprod_s(a, b));
7402 Note the difference between @code{sprod} which is the algebraic class, and
7403 @code{sprod_s} which is the unadorned C++ structure containing the @code{left}
7404 and @code{right} data members. As shown above, an @code{sprod} can be
7405 constructed from an @code{sprod_s} object.
7407 If you find the nested @code{sprod(sprod_s())} constructor too unwieldy,
7408 you could define a little wrapper function like this:
7411 inline ex make_sprod(ex left, ex right)
7413 return sprod(sprod_s(left, right));
7417 The @code{sprod_s} object contained in @code{sprod} can be accessed with
7418 the GiNaC @code{ex_to<>()} function followed by the @code{->} operator or
7419 @code{get_struct()}:
7423 cout << ex_to<sprod>(e)->left << endl;
7425 cout << ex_to<sprod>(e).get_struct().right << endl;
7430 You only have read access to the members of @code{sprod_s}.
7432 The type definition of @code{sprod} is enough to write your own algorithms
7433 that deal with scalar products, for example:
7438 if (is_a<sprod>(p)) @{
7439 const sprod_s & sp = ex_to<sprod>(p).get_struct();
7440 return make_sprod(sp.right, sp.left);
7451 @subsection Structure output
7453 While the @code{sprod} type is useable it still leaves something to be
7454 desired, most notably proper output:
7459 // -> [structure object]
7463 By default, any structure types you define will be printed as
7464 @samp{[structure object]}. To override this you can either specialize the
7465 template's @code{print()} member function, or specify print methods with
7466 @code{set_print_func<>()}, as described in @ref{Printing}. Unfortunately,
7467 it's not possible to supply class options like @code{print_func<>()} to
7468 structures, so for a self-contained structure type you need to resort to
7469 overriding the @code{print()} function, which is also what we will do here.
7471 The member functions of GiNaC classes are described in more detail in the
7472 next section, but it shouldn't be hard to figure out what's going on here:
7475 void sprod::print(const print_context & c, unsigned level) const
7477 // tree debug output handled by superclass
7478 if (is_a<print_tree>(c))
7479 inherited::print(c, level);
7481 // get the contained sprod_s object
7482 const sprod_s & sp = get_struct();
7484 // print_context::s is a reference to an ostream
7485 c.s << "<" << sp.left << "|" << sp.right << ">";
7489 Now we can print expressions containing scalar products:
7495 cout << swap_sprod(e) << endl;
7500 @subsection Comparing structures
7502 The @code{sprod} class defined so far still has one important drawback: all
7503 scalar products are treated as being equal because GiNaC doesn't know how to
7504 compare objects of type @code{sprod_s}. This can lead to some confusing
7505 and undesired behavior:
7509 cout << make_sprod(a, b) - make_sprod(a*a, b*b) << endl;
7511 cout << make_sprod(a, b) + make_sprod(a*a, b*b) << endl;
7512 // -> 2*<a|b> or 2*<a^2|b^2> (which one is undefined)
7516 To remedy this, we first need to define the operators @code{==} and @code{<}
7517 for objects of type @code{sprod_s}:
7520 inline bool operator==(const sprod_s & lhs, const sprod_s & rhs)
7522 return lhs.left.is_equal(rhs.left) && lhs.right.is_equal(rhs.right);
7525 inline bool operator<(const sprod_s & lhs, const sprod_s & rhs)
7527 return lhs.left.compare(rhs.left) < 0
7528 ? true : lhs.right.compare(rhs.right) < 0;
7532 The ordering established by the @code{<} operator doesn't have to make any
7533 algebraic sense, but it needs to be well defined. Note that we can't use
7534 expressions like @code{lhs.left == rhs.left} or @code{lhs.left < rhs.left}
7535 in the implementation of these operators because they would construct
7536 GiNaC @code{relational} objects which in the case of @code{<} do not
7537 establish a well defined ordering (for arbitrary expressions, GiNaC can't
7538 decide which one is algebraically 'less').
7540 Next, we need to change our definition of the @code{sprod} type to let
7541 GiNaC know that an ordering relation exists for the embedded objects:
7544 typedef structure<sprod_s, compare_std_less> sprod;
7547 @code{sprod} objects then behave as expected:
7551 cout << make_sprod(a, b) - make_sprod(a*a, b*b) << endl;
7552 // -> <a|b>-<a^2|b^2>
7553 cout << make_sprod(a, b) + make_sprod(a*a, b*b) << endl;
7554 // -> <a|b>+<a^2|b^2>
7555 cout << make_sprod(a, b) - make_sprod(a, b) << endl;
7557 cout << make_sprod(a, b) + make_sprod(a, b) << endl;
7562 The @code{compare_std_less} policy parameter tells GiNaC to use the
7563 @code{std::less} and @code{std::equal_to} functors to compare objects of
7564 type @code{sprod_s}. By default, these functors forward their work to the
7565 standard @code{<} and @code{==} operators, which we have overloaded.
7566 Alternatively, we could have specialized @code{std::less} and
7567 @code{std::equal_to} for class @code{sprod_s}.
7569 GiNaC provides two other comparison policies for @code{structure<T>}
7570 objects: the default @code{compare_all_equal}, and @code{compare_bitwise}
7571 which does a bit-wise comparison of the contained @code{T} objects.
7572 This should be used with extreme care because it only works reliably with
7573 built-in integral types, and it also compares any padding (filler bytes of
7574 undefined value) that the @code{T} class might have.
7576 @subsection Subexpressions
7578 Our scalar product class has two subexpressions: the left and right
7579 operands. It might be a good idea to make them accessible via the standard
7580 @code{nops()} and @code{op()} methods:
7583 size_t sprod::nops() const
7588 ex sprod::op(size_t i) const
7592 return get_struct().left;
7594 return get_struct().right;
7596 throw std::range_error("sprod::op(): no such operand");
7601 Implementing @code{nops()} and @code{op()} for container types such as
7602 @code{sprod} has two other nice side effects:
7606 @code{has()} works as expected
7608 GiNaC generates better hash keys for the objects (the default implementation
7609 of @code{calchash()} takes subexpressions into account)
7612 @cindex @code{let_op()}
7613 There is a non-const variant of @code{op()} called @code{let_op()} that
7614 allows replacing subexpressions:
7617 ex & sprod::let_op(size_t i)
7619 // every non-const member function must call this
7620 ensure_if_modifiable();
7624 return get_struct().left;
7626 return get_struct().right;
7628 throw std::range_error("sprod::let_op(): no such operand");
7633 Once we have provided @code{let_op()} we also get @code{subs()} and
7634 @code{map()} for free. In fact, every container class that returns a non-null
7635 @code{nops()} value must either implement @code{let_op()} or provide custom
7636 implementations of @code{subs()} and @code{map()}.
7638 In turn, the availability of @code{map()} enables the recursive behavior of a
7639 couple of other default method implementations, in particular @code{evalf()},
7640 @code{evalm()}, @code{normal()}, @code{diff()} and @code{expand()}. Although
7641 we probably want to provide our own version of @code{expand()} for scalar
7642 products that turns expressions like @samp{<a+b|c>} into @samp{<a|c>+<b|c>}.
7643 This is left as an exercise for the reader.
7645 The @code{structure<T>} template defines many more member functions that
7646 you can override by specialization to customize the behavior of your
7647 structures. You are referred to the next section for a description of
7648 some of these (especially @code{eval()}). There is, however, one topic
7649 that shall be addressed here, as it demonstrates one peculiarity of the
7650 @code{structure<T>} template: archiving.
7652 @subsection Archiving structures
7654 If you don't know how the archiving of GiNaC objects is implemented, you
7655 should first read the next section and then come back here. You're back?
7658 To implement archiving for structures it is not enough to provide
7659 specializations for the @code{archive()} member function and the
7660 unarchiving constructor (the @code{unarchive()} function has a default
7661 implementation). You also need to provide a unique name (as a string literal)
7662 for each structure type you define. This is because in GiNaC archives,
7663 the class of an object is stored as a string, the class name.
7665 By default, this class name (as returned by the @code{class_name()} member
7666 function) is @samp{structure} for all structure classes. This works as long
7667 as you have only defined one structure type, but if you use two or more you
7668 need to provide a different name for each by specializing the
7669 @code{get_class_name()} member function. Here is a sample implementation
7670 for enabling archiving of the scalar product type defined above:
7673 const char *sprod::get_class_name() @{ return "sprod"; @}
7675 void sprod::archive(archive_node & n) const
7677 inherited::archive(n);
7678 n.add_ex("left", get_struct().left);
7679 n.add_ex("right", get_struct().right);
7682 sprod::structure(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst)
7684 n.find_ex("left", get_struct().left, sym_lst);
7685 n.find_ex("right", get_struct().right, sym_lst);
7689 Note that the unarchiving constructor is @code{sprod::structure} and not
7690 @code{sprod::sprod}, and that we don't need to supply an
7691 @code{sprod::unarchive()} function.
7694 @node Adding classes, A comparison with other CAS, Structures, Extending GiNaC
7695 @c node-name, next, previous, up
7696 @section Adding classes
7698 The @code{structure<T>} template provides an way to extend GiNaC with custom
7699 algebraic classes that is easy to use but has its limitations, the most
7700 severe of which being that you can't add any new member functions to
7701 structures. To be able to do this, you need to write a new class definition
7704 This section will explain how to implement new algebraic classes in GiNaC by
7705 giving the example of a simple 'string' class. After reading this section
7706 you will know how to properly declare a GiNaC class and what the minimum
7707 required member functions are that you have to implement. We only cover the
7708 implementation of a 'leaf' class here (i.e. one that doesn't contain
7709 subexpressions). Creating a container class like, for example, a class
7710 representing tensor products is more involved but this section should give
7711 you enough information so you can consult the source to GiNaC's predefined
7712 classes if you want to implement something more complicated.
7714 @subsection GiNaC's run-time type information system
7716 @cindex hierarchy of classes
7718 All algebraic classes (that is, all classes that can appear in expressions)
7719 in GiNaC are direct or indirect subclasses of the class @code{basic}. So a
7720 @code{basic *} (which is essentially what an @code{ex} is) represents a
7721 generic pointer to an algebraic class. Occasionally it is necessary to find
7722 out what the class of an object pointed to by a @code{basic *} really is.
7723 Also, for the unarchiving of expressions it must be possible to find the
7724 @code{unarchive()} function of a class given the class name (as a string). A
7725 system that provides this kind of information is called a run-time type
7726 information (RTTI) system. The C++ language provides such a thing (see the
7727 standard header file @file{<typeinfo>}) but for efficiency reasons GiNaC
7728 implements its own, simpler RTTI.
7730 The RTTI in GiNaC is based on two mechanisms:
7735 The @code{basic} class declares a member variable @code{tinfo_key} which
7736 holds an unsigned integer that identifies the object's class. These numbers
7737 are defined in the @file{tinfos.h} header file for the built-in GiNaC
7738 classes. They all start with @code{TINFO_}.
7741 By means of some clever tricks with static members, GiNaC maintains a list
7742 of information for all classes derived from @code{basic}. The information
7743 available includes the class names, the @code{tinfo_key}s, and pointers
7744 to the unarchiving functions. This class registry is defined in the
7745 @file{registrar.h} header file.
7749 The disadvantage of this proprietary RTTI implementation is that there's
7750 a little more to do when implementing new classes (C++'s RTTI works more
7751 or less automatically) but don't worry, most of the work is simplified by
7754 @subsection A minimalistic example
7756 Now we will start implementing a new class @code{mystring} that allows
7757 placing character strings in algebraic expressions (this is not very useful,
7758 but it's just an example). This class will be a direct subclass of
7759 @code{basic}. You can use this sample implementation as a starting point
7760 for your own classes.
7762 The code snippets given here assume that you have included some header files
7768 #include <stdexcept>
7769 using namespace std;
7771 #include <ginac/ginac.h>
7772 using namespace GiNaC;
7775 The first thing we have to do is to define a @code{tinfo_key} for our new
7776 class. This can be any arbitrary unsigned number that is not already taken
7777 by one of the existing classes but it's better to come up with something
7778 that is unlikely to clash with keys that might be added in the future. The
7779 numbers in @file{tinfos.h} are modeled somewhat after the class hierarchy
7780 which is not a requirement but we are going to stick with this scheme:
7783 const unsigned TINFO_mystring = 0x42420001U;
7786 Now we can write down the class declaration. The class stores a C++
7787 @code{string} and the user shall be able to construct a @code{mystring}
7788 object from a C or C++ string:
7791 class mystring : public basic
7793 GINAC_DECLARE_REGISTERED_CLASS(mystring, basic)
7796 mystring(const string &s);
7797 mystring(const char *s);
7803 GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
7806 The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
7807 macros are defined in @file{registrar.h}. They take the name of the class
7808 and its direct superclass as arguments and insert all required declarations
7809 for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
7810 the first line after the opening brace of the class definition. The
7811 @code{GINAC_IMPLEMENT_REGISTERED_CLASS} may appear anywhere else in the
7812 source (at global scope, of course, not inside a function).
7814 @code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
7815 declarations of the default constructor and a couple of other functions that
7816 are required. It also defines a type @code{inherited} which refers to the
7817 superclass so you don't have to modify your code every time you shuffle around
7818 the class hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} registers the
7819 class with the GiNaC RTTI (there is also a
7820 @code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT} which allows specifying additional
7821 options for the class, and which we will be using instead in a few minutes).
7823 Now there are seven member functions we have to implement to get a working
7829 @code{mystring()}, the default constructor.
7832 @code{void archive(archive_node &n)}, the archiving function. This stores all
7833 information needed to reconstruct an object of this class inside an
7834 @code{archive_node}.
7837 @code{mystring(const archive_node &n, lst &sym_lst)}, the unarchiving
7838 constructor. This constructs an instance of the class from the information
7839 found in an @code{archive_node}.
7842 @code{ex unarchive(const archive_node &n, lst &sym_lst)}, the static
7843 unarchiving function. It constructs a new instance by calling the unarchiving
7847 @cindex @code{compare_same_type()}
7848 @code{int compare_same_type(const basic &other)}, which is used internally
7849 by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
7850 -1, depending on the relative order of this object and the @code{other}
7851 object. If it returns 0, the objects are considered equal.
7852 @strong{Please notice:} This has nothing to do with the (numeric) ordering
7853 relationship expressed by @code{<}, @code{>=} etc (which cannot be defined
7854 for non-numeric classes). For example, @code{numeric(1).compare_same_type(numeric(2))}
7855 may return +1 even though 1 is clearly smaller than 2. Every GiNaC class
7856 must provide a @code{compare_same_type()} function, even those representing
7857 objects for which no reasonable algebraic ordering relationship can be
7861 And, of course, @code{mystring(const string &s)} and @code{mystring(const char *s)}
7862 which are the two constructors we declared.
7866 Let's proceed step-by-step. The default constructor looks like this:
7869 mystring::mystring() : inherited(TINFO_mystring) @{@}
7872 The golden rule is that in all constructors you have to set the
7873 @code{tinfo_key} member to the @code{TINFO_*} value of your class. Otherwise
7874 it will be set by the constructor of the superclass and all hell will break
7875 loose in the RTTI. For your convenience, the @code{basic} class provides
7876 a constructor that takes a @code{tinfo_key} value, which we are using here
7877 (remember that in our case @code{inherited == basic}). If the superclass
7878 didn't have such a constructor, we would have to set the @code{tinfo_key}
7879 to the right value manually.
7881 In the default constructor you should set all other member variables to
7882 reasonable default values (we don't need that here since our @code{str}
7883 member gets set to an empty string automatically).
7885 Next are the three functions for archiving. You have to implement them even
7886 if you don't plan to use archives, but the minimum required implementation
7887 is really simple. First, the archiving function:
7890 void mystring::archive(archive_node &n) const
7892 inherited::archive(n);
7893 n.add_string("string", str);
7897 The only thing that is really required is calling the @code{archive()}
7898 function of the superclass. Optionally, you can store all information you
7899 deem necessary for representing the object into the passed
7900 @code{archive_node}. We are just storing our string here. For more
7901 information on how the archiving works, consult the @file{archive.h} header
7904 The unarchiving constructor is basically the inverse of the archiving
7908 mystring::mystring(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
7910 n.find_string("string", str);
7914 If you don't need archiving, just leave this function empty (but you must
7915 invoke the unarchiving constructor of the superclass). Note that we don't
7916 have to set the @code{tinfo_key} here because it is done automatically
7917 by the unarchiving constructor of the @code{basic} class.
7919 Finally, the unarchiving function:
7922 ex mystring::unarchive(const archive_node &n, lst &sym_lst)
7924 return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
7928 You don't have to understand how exactly this works. Just copy these
7929 four lines into your code literally (replacing the class name, of
7930 course). It calls the unarchiving constructor of the class and unless
7931 you are doing something very special (like matching @code{archive_node}s
7932 to global objects) you don't need a different implementation. For those
7933 who are interested: setting the @code{dynallocated} flag puts the object
7934 under the control of GiNaC's garbage collection. It will get deleted
7935 automatically once it is no longer referenced.
7937 Our @code{compare_same_type()} function uses a provided function to compare
7941 int mystring::compare_same_type(const basic &other) const
7943 const mystring &o = static_cast<const mystring &>(other);
7944 int cmpval = str.compare(o.str);
7947 else if (cmpval < 0)
7954 Although this function takes a @code{basic &}, it will always be a reference
7955 to an object of exactly the same class (objects of different classes are not
7956 comparable), so the cast is safe. If this function returns 0, the two objects
7957 are considered equal (in the sense that @math{A-B=0}), so you should compare
7958 all relevant member variables.
7960 Now the only thing missing is our two new constructors:
7963 mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s) @{@}
7964 mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s) @{@}
7967 No surprises here. We set the @code{str} member from the argument and
7968 remember to pass the right @code{tinfo_key} to the @code{basic} constructor.
7970 That's it! We now have a minimal working GiNaC class that can store
7971 strings in algebraic expressions. Let's confirm that the RTTI works:
7974 ex e = mystring("Hello, world!");
7975 cout << is_a<mystring>(e) << endl;
7978 cout << ex_to<basic>(e).class_name() << endl;
7982 Obviously it does. Let's see what the expression @code{e} looks like:
7986 // -> [mystring object]
7989 Hm, not exactly what we expect, but of course the @code{mystring} class
7990 doesn't yet know how to print itself. This can be done either by implementing
7991 the @code{print()} member function, or, preferably, by specifying a
7992 @code{print_func<>()} class option. Let's say that we want to print the string
7993 surrounded by double quotes:
7996 class mystring : public basic
8000 void do_print(const print_context &c, unsigned level = 0) const;
8004 void mystring::do_print(const print_context &c, unsigned level) const
8006 // print_context::s is a reference to an ostream
8007 c.s << '\"' << str << '\"';
8011 The @code{level} argument is only required for container classes to
8012 correctly parenthesize the output.
8014 Now we need to tell GiNaC that @code{mystring} objects should use the
8015 @code{do_print()} member function for printing themselves. For this, we
8019 GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
8025 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mystring, basic,
8026 print_func<print_context>(&mystring::do_print))
8029 Let's try again to print the expression:
8033 // -> "Hello, world!"
8036 Much better. If we wanted to have @code{mystring} objects displayed in a
8037 different way depending on the output format (default, LaTeX, etc.), we
8038 would have supplied multiple @code{print_func<>()} options with different
8039 template parameters (@code{print_dflt}, @code{print_latex}, etc.),
8040 separated by dots. This is similar to the way options are specified for
8041 symbolic functions. @xref{Printing}, for a more in-depth description of the
8042 way expression output is implemented in GiNaC.
8044 The @code{mystring} class can be used in arbitrary expressions:
8047 e += mystring("GiNaC rulez");
8049 // -> "GiNaC rulez"+"Hello, world!"
8052 (GiNaC's automatic term reordering is in effect here), or even
8055 e = pow(mystring("One string"), 2*sin(Pi-mystring("Another string")));
8057 // -> "One string"^(2*sin(-"Another string"+Pi))
8060 Whether this makes sense is debatable but remember that this is only an
8061 example. At least it allows you to implement your own symbolic algorithms
8064 Note that GiNaC's algebraic rules remain unchanged:
8067 e = mystring("Wow") * mystring("Wow");
8071 e = pow(mystring("First")-mystring("Second"), 2);
8072 cout << e.expand() << endl;
8073 // -> -2*"First"*"Second"+"First"^2+"Second"^2
8076 There's no way to, for example, make GiNaC's @code{add} class perform string
8077 concatenation. You would have to implement this yourself.
8079 @subsection Automatic evaluation
8082 @cindex @code{eval()}
8083 @cindex @code{hold()}
8084 When dealing with objects that are just a little more complicated than the
8085 simple string objects we have implemented, chances are that you will want to
8086 have some automatic simplifications or canonicalizations performed on them.
8087 This is done in the evaluation member function @code{eval()}. Let's say that
8088 we wanted all strings automatically converted to lowercase with
8089 non-alphabetic characters stripped, and empty strings removed:
8092 class mystring : public basic
8096 ex eval(int level = 0) const;
8100 ex mystring::eval(int level) const
8103 for (int i=0; i<str.length(); i++) @{
8105 if (c >= 'A' && c <= 'Z')
8106 new_str += tolower(c);
8107 else if (c >= 'a' && c <= 'z')
8111 if (new_str.length() == 0)
8114 return mystring(new_str).hold();
8118 The @code{level} argument is used to limit the recursion depth of the
8119 evaluation. We don't have any subexpressions in the @code{mystring}
8120 class so we are not concerned with this. If we had, we would call the
8121 @code{eval()} functions of the subexpressions with @code{level - 1} as
8122 the argument if @code{level != 1}. The @code{hold()} member function
8123 sets a flag in the object that prevents further evaluation. Otherwise
8124 we might end up in an endless loop. When you want to return the object
8125 unmodified, use @code{return this->hold();}.
8127 Let's confirm that it works:
8130 ex e = mystring("Hello, world!") + mystring("!?#");
8134 e = mystring("Wow!") + mystring("WOW") + mystring(" W ** o ** W");
8139 @subsection Optional member functions
8141 We have implemented only a small set of member functions to make the class
8142 work in the GiNaC framework. There are two functions that are not strictly
8143 required but will make operations with objects of the class more efficient:
8145 @cindex @code{calchash()}
8146 @cindex @code{is_equal_same_type()}
8148 unsigned calchash() const;
8149 bool is_equal_same_type(const basic &other) const;
8152 The @code{calchash()} method returns an @code{unsigned} hash value for the
8153 object which will allow GiNaC to compare and canonicalize expressions much
8154 more efficiently. You should consult the implementation of some of the built-in
8155 GiNaC classes for examples of hash functions. The default implementation of
8156 @code{calchash()} calculates a hash value out of the @code{tinfo_key} of the
8157 class and all subexpressions that are accessible via @code{op()}.
8159 @code{is_equal_same_type()} works like @code{compare_same_type()} but only
8160 tests for equality without establishing an ordering relation, which is often
8161 faster. The default implementation of @code{is_equal_same_type()} just calls
8162 @code{compare_same_type()} and tests its result for zero.
8164 @subsection Other member functions
8166 For a real algebraic class, there are probably some more functions that you
8167 might want to provide:
8170 bool info(unsigned inf) const;
8171 ex evalf(int level = 0) const;
8172 ex series(const relational & r, int order, unsigned options = 0) const;
8173 ex derivative(const symbol & s) const;
8176 If your class stores sub-expressions (see the scalar product example in the
8177 previous section) you will probably want to override
8179 @cindex @code{let_op()}
8182 ex op(size_t i) const;
8183 ex & let_op(size_t i);
8184 ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
8185 ex map(map_function & f) const;
8188 @code{let_op()} is a variant of @code{op()} that allows write access. The
8189 default implementations of @code{subs()} and @code{map()} use it, so you have
8190 to implement either @code{let_op()}, or @code{subs()} and @code{map()}.
8192 You can, of course, also add your own new member functions. Remember
8193 that the RTTI may be used to get information about what kinds of objects
8194 you are dealing with (the position in the class hierarchy) and that you
8195 can always extract the bare object from an @code{ex} by stripping the
8196 @code{ex} off using the @code{ex_to<mystring>(e)} function when that
8197 should become a need.
8199 That's it. May the source be with you!
8202 @node A comparison with other CAS, Advantages, Adding classes, Top
8203 @c node-name, next, previous, up
8204 @chapter A Comparison With Other CAS
8207 This chapter will give you some information on how GiNaC compares to
8208 other, traditional Computer Algebra Systems, like @emph{Maple},
8209 @emph{Mathematica} or @emph{Reduce}, where it has advantages and
8210 disadvantages over these systems.
8213 * Advantages:: Strengths of the GiNaC approach.
8214 * Disadvantages:: Weaknesses of the GiNaC approach.
8215 * Why C++?:: Attractiveness of C++.
8218 @node Advantages, Disadvantages, A comparison with other CAS, A comparison with other CAS
8219 @c node-name, next, previous, up
8222 GiNaC has several advantages over traditional Computer
8223 Algebra Systems, like
8228 familiar language: all common CAS implement their own proprietary
8229 grammar which you have to learn first (and maybe learn again when your
8230 vendor decides to `enhance' it). With GiNaC you can write your program
8231 in common C++, which is standardized.
8235 structured data types: you can build up structured data types using
8236 @code{struct}s or @code{class}es together with STL features instead of
8237 using unnamed lists of lists of lists.
8240 strongly typed: in CAS, you usually have only one kind of variables
8241 which can hold contents of an arbitrary type. This 4GL like feature is
8242 nice for novice programmers, but dangerous.
8245 development tools: powerful development tools exist for C++, like fancy
8246 editors (e.g. with automatic indentation and syntax highlighting),
8247 debuggers, visualization tools, documentation generators@dots{}
8250 modularization: C++ programs can easily be split into modules by
8251 separating interface and implementation.
8254 price: GiNaC is distributed under the GNU Public License which means
8255 that it is free and available with source code. And there are excellent
8256 C++-compilers for free, too.
8259 extendable: you can add your own classes to GiNaC, thus extending it on
8260 a very low level. Compare this to a traditional CAS that you can
8261 usually only extend on a high level by writing in the language defined
8262 by the parser. In particular, it turns out to be almost impossible to
8263 fix bugs in a traditional system.
8266 multiple interfaces: Though real GiNaC programs have to be written in
8267 some editor, then be compiled, linked and executed, there are more ways
8268 to work with the GiNaC engine. Many people want to play with
8269 expressions interactively, as in traditional CASs. Currently, two such
8270 windows into GiNaC have been implemented and many more are possible: the
8271 tiny @command{ginsh} that is part of the distribution exposes GiNaC's
8272 types to a command line and second, as a more consistent approach, an
8273 interactive interface to the Cint C++ interpreter has been put together
8274 (called GiNaC-cint) that allows an interactive scripting interface
8275 consistent with the C++ language. It is available from the usual GiNaC
8279 seamless integration: it is somewhere between difficult and impossible
8280 to call CAS functions from within a program written in C++ or any other
8281 programming language and vice versa. With GiNaC, your symbolic routines
8282 are part of your program. You can easily call third party libraries,
8283 e.g. for numerical evaluation or graphical interaction. All other
8284 approaches are much more cumbersome: they range from simply ignoring the
8285 problem (i.e. @emph{Maple}) to providing a method for `embedding' the
8286 system (i.e. @emph{Yacas}).
8289 efficiency: often large parts of a program do not need symbolic
8290 calculations at all. Why use large integers for loop variables or
8291 arbitrary precision arithmetics where @code{int} and @code{double} are
8292 sufficient? For pure symbolic applications, GiNaC is comparable in
8293 speed with other CAS.
8298 @node Disadvantages, Why C++?, Advantages, A comparison with other CAS
8299 @c node-name, next, previous, up
8300 @section Disadvantages
8302 Of course it also has some disadvantages:
8307 advanced features: GiNaC cannot compete with a program like
8308 @emph{Reduce} which exists for more than 30 years now or @emph{Maple}
8309 which grows since 1981 by the work of dozens of programmers, with
8310 respect to mathematical features. Integration, factorization,
8311 non-trivial simplifications, limits etc. are missing in GiNaC (and are
8312 not planned for the near future).
8315 portability: While the GiNaC library itself is designed to avoid any
8316 platform dependent features (it should compile on any ANSI compliant C++
8317 compiler), the currently used version of the CLN library (fast large
8318 integer and arbitrary precision arithmetics) can only by compiled
8319 without hassle on systems with the C++ compiler from the GNU Compiler
8320 Collection (GCC).@footnote{This is because CLN uses PROVIDE/REQUIRE like
8321 macros to let the compiler gather all static initializations, which
8322 works for GNU C++ only. Feel free to contact the authors in case you
8323 really believe that you need to use a different compiler. We have
8324 occasionally used other compilers and may be able to give you advice.}
8325 GiNaC uses recent language features like explicit constructors, mutable
8326 members, RTTI, @code{dynamic_cast}s and STL, so ANSI compliance is meant
8327 literally. Recent GCC versions starting at 2.95.3, although itself not
8328 yet ANSI compliant, support all needed features.
8333 @node Why C++?, Internal structures, Disadvantages, A comparison with other CAS
8334 @c node-name, next, previous, up
8337 Why did we choose to implement GiNaC in C++ instead of Java or any other
8338 language? C++ is not perfect: type checking is not strict (casting is
8339 possible), separation between interface and implementation is not
8340 complete, object oriented design is not enforced. The main reason is
8341 the often scolded feature of operator overloading in C++. While it may
8342 be true that operating on classes with a @code{+} operator is rarely
8343 meaningful, it is perfectly suited for algebraic expressions. Writing
8344 @math{3x+5y} as @code{3*x+5*y} instead of
8345 @code{x.times(3).plus(y.times(5))} looks much more natural.
8346 Furthermore, the main developers are more familiar with C++ than with
8347 any other programming language.
8350 @node Internal structures, Expressions are reference counted, Why C++? , Top
8351 @c node-name, next, previous, up
8352 @appendix Internal structures
8355 * Expressions are reference counted::
8356 * Internal representation of products and sums::
8359 @node Expressions are reference counted, Internal representation of products and sums, Internal structures, Internal structures
8360 @c node-name, next, previous, up
8361 @appendixsection Expressions are reference counted
8363 @cindex reference counting
8364 @cindex copy-on-write
8365 @cindex garbage collection
8366 In GiNaC, there is an @emph{intrusive reference-counting} mechanism at work
8367 where the counter belongs to the algebraic objects derived from class
8368 @code{basic} but is maintained by the smart pointer class @code{ptr}, of
8369 which @code{ex} contains an instance. If you understood that, you can safely
8370 skip the rest of this passage.
8372 Expressions are extremely light-weight since internally they work like
8373 handles to the actual representation. They really hold nothing more
8374 than a pointer to some other object. What this means in practice is
8375 that whenever you create two @code{ex} and set the second equal to the
8376 first no copying process is involved. Instead, the copying takes place
8377 as soon as you try to change the second. Consider the simple sequence
8382 #include <ginac/ginac.h>
8383 using namespace std;
8384 using namespace GiNaC;
8388 symbol x("x"), y("y"), z("z");
8391 e1 = sin(x + 2*y) + 3*z + 41;
8392 e2 = e1; // e2 points to same object as e1
8393 cout << e2 << endl; // prints sin(x+2*y)+3*z+41
8394 e2 += 1; // e2 is copied into a new object
8395 cout << e2 << endl; // prints sin(x+2*y)+3*z+42
8399 The line @code{e2 = e1;} creates a second expression pointing to the
8400 object held already by @code{e1}. The time involved for this operation
8401 is therefore constant, no matter how large @code{e1} was. Actual
8402 copying, however, must take place in the line @code{e2 += 1;} because
8403 @code{e1} and @code{e2} are not handles for the same object any more.
8404 This concept is called @dfn{copy-on-write semantics}. It increases
8405 performance considerably whenever one object occurs multiple times and
8406 represents a simple garbage collection scheme because when an @code{ex}
8407 runs out of scope its destructor checks whether other expressions handle
8408 the object it points to too and deletes the object from memory if that
8409 turns out not to be the case. A slightly less trivial example of
8410 differentiation using the chain-rule should make clear how powerful this
8415 symbol x("x"), y("y");
8419 ex e3 = diff(sin(e2), x); // first derivative of sin(e2) by x
8420 cout << e1 << endl // prints x+3*y
8421 << e2 << endl // prints (x+3*y)^3
8422 << e3 << endl; // prints 3*(x+3*y)^2*cos((x+3*y)^3)
8426 Here, @code{e1} will actually be referenced three times while @code{e2}
8427 will be referenced two times. When the power of an expression is built,
8428 that expression needs not be copied. Likewise, since the derivative of
8429 a power of an expression can be easily expressed in terms of that
8430 expression, no copying of @code{e1} is involved when @code{e3} is
8431 constructed. So, when @code{e3} is constructed it will print as
8432 @code{3*(x+3*y)^2*cos((x+3*y)^3)} but the argument of @code{cos()} only
8433 holds a reference to @code{e2} and the factor in front is just
8436 As a user of GiNaC, you cannot see this mechanism of copy-on-write
8437 semantics. When you insert an expression into a second expression, the
8438 result behaves exactly as if the contents of the first expression were
8439 inserted. But it may be useful to remember that this is not what
8440 happens. Knowing this will enable you to write much more efficient
8441 code. If you still have an uncertain feeling with copy-on-write
8442 semantics, we recommend you have a look at the
8443 @uref{http://www.parashift.com/c++-faq-lite/, C++-FAQ lite} by
8444 Marshall Cline. Chapter 16 covers this issue and presents an
8445 implementation which is pretty close to the one in GiNaC.
8448 @node Internal representation of products and sums, Package tools, Expressions are reference counted, Internal structures
8449 @c node-name, next, previous, up
8450 @appendixsection Internal representation of products and sums
8452 @cindex representation
8455 @cindex @code{power}
8456 Although it should be completely transparent for the user of
8457 GiNaC a short discussion of this topic helps to understand the sources
8458 and also explain performance to a large degree. Consider the
8459 unexpanded symbolic expression
8461 $2d^3 \left( 4a + 5b - 3 \right)$
8464 @math{2*d^3*(4*a+5*b-3)}
8466 which could naively be represented by a tree of linear containers for
8467 addition and multiplication, one container for exponentiation with base
8468 and exponent and some atomic leaves of symbols and numbers in this
8473 @cindex pair-wise representation
8474 However, doing so results in a rather deeply nested tree which will
8475 quickly become inefficient to manipulate. We can improve on this by
8476 representing the sum as a sequence of terms, each one being a pair of a
8477 purely numeric multiplicative coefficient and its rest. In the same
8478 spirit we can store the multiplication as a sequence of terms, each
8479 having a numeric exponent and a possibly complicated base, the tree
8480 becomes much more flat:
8484 The number @code{3} above the symbol @code{d} shows that @code{mul}
8485 objects are treated similarly where the coefficients are interpreted as
8486 @emph{exponents} now. Addition of sums of terms or multiplication of
8487 products with numerical exponents can be coded to be very efficient with
8488 such a pair-wise representation. Internally, this handling is performed
8489 by most CAS in this way. It typically speeds up manipulations by an
8490 order of magnitude. The overall multiplicative factor @code{2} and the
8491 additive term @code{-3} look somewhat out of place in this
8492 representation, however, since they are still carrying a trivial
8493 exponent and multiplicative factor @code{1} respectively. Within GiNaC,
8494 this is avoided by adding a field that carries an overall numeric
8495 coefficient. This results in the realistic picture of internal
8498 $2d^3 \left( 4a + 5b - 3 \right)$:
8501 @math{2*d^3*(4*a+5*b-3)}:
8507 This also allows for a better handling of numeric radicals, since
8508 @code{sqrt(2)} can now be carried along calculations. Now it should be
8509 clear, why both classes @code{add} and @code{mul} are derived from the
8510 same abstract class: the data representation is the same, only the
8511 semantics differs. In the class hierarchy, methods for polynomial
8512 expansion and the like are reimplemented for @code{add} and @code{mul},
8513 but the data structure is inherited from @code{expairseq}.
8516 @node Package tools, ginac-config, Internal representation of products and sums, Top
8517 @c node-name, next, previous, up
8518 @appendix Package tools
8520 If you are creating a software package that uses the GiNaC library,
8521 setting the correct command line options for the compiler and linker
8522 can be difficult. GiNaC includes two tools to make this process easier.
8525 * ginac-config:: A shell script to detect compiler and linker flags.
8526 * AM_PATH_GINAC:: Macro for GNU automake.
8530 @node ginac-config, AM_PATH_GINAC, Package tools, Package tools
8531 @c node-name, next, previous, up
8532 @section @command{ginac-config}
8533 @cindex ginac-config
8535 @command{ginac-config} is a shell script that you can use to determine
8536 the compiler and linker command line options required to compile and
8537 link a program with the GiNaC library.
8539 @command{ginac-config} takes the following flags:
8543 Prints out the version of GiNaC installed.
8545 Prints '-I' flags pointing to the installed header files.
8547 Prints out the linker flags necessary to link a program against GiNaC.
8548 @item --prefix[=@var{PREFIX}]
8549 If @var{PREFIX} is specified, overrides the configured value of @env{$prefix}.
8550 (And of exec-prefix, unless @code{--exec-prefix} is also specified)
8551 Otherwise, prints out the configured value of @env{$prefix}.
8552 @item --exec-prefix[=@var{PREFIX}]
8553 If @var{PREFIX} is specified, overrides the configured value of @env{$exec_prefix}.
8554 Otherwise, prints out the configured value of @env{$exec_prefix}.
8557 Typically, @command{ginac-config} will be used within a configure
8558 script, as described below. It, however, can also be used directly from
8559 the command line using backquotes to compile a simple program. For
8563 c++ -o simple `ginac-config --cppflags` simple.cpp `ginac-config --libs`
8566 This command line might expand to (for example):
8569 cc -o simple -I/usr/local/include simple.cpp -L/usr/local/lib \
8570 -lginac -lcln -lstdc++
8573 Not only is the form using @command{ginac-config} easier to type, it will
8574 work on any system, no matter how GiNaC was configured.
8577 @node AM_PATH_GINAC, Configure script options, ginac-config, Package tools
8578 @c node-name, next, previous, up
8579 @section @samp{AM_PATH_GINAC}
8580 @cindex AM_PATH_GINAC
8582 For packages configured using GNU automake, GiNaC also provides
8583 a macro to automate the process of checking for GiNaC.
8586 AM_PATH_GINAC([@var{MINIMUM-VERSION}, [@var{ACTION-IF-FOUND}
8587 [, @var{ACTION-IF-NOT-FOUND}]]])
8595 Determines the location of GiNaC using @command{ginac-config}, which is
8596 either found in the user's path, or from the environment variable
8597 @env{GINACLIB_CONFIG}.
8600 Tests the installed libraries to make sure that their version
8601 is later than @var{MINIMUM-VERSION}. (A default version will be used
8605 If the required version was found, sets the @env{GINACLIB_CPPFLAGS} variable
8606 to the output of @command{ginac-config --cppflags} and the @env{GINACLIB_LIBS}
8607 variable to the output of @command{ginac-config --libs}, and calls
8608 @samp{AC_SUBST()} for these variables so they can be used in generated
8609 makefiles, and then executes @var{ACTION-IF-FOUND}.
8612 If the required version was not found, sets @env{GINACLIB_CPPFLAGS} and
8613 @env{GINACLIB_LIBS} to empty strings, and executes @var{ACTION-IF-NOT-FOUND}.
8617 This macro is in file @file{ginac.m4} which is installed in
8618 @file{$datadir/aclocal}. Note that if automake was installed with a
8619 different @samp{--prefix} than GiNaC, you will either have to manually
8620 move @file{ginac.m4} to automake's @file{$datadir/aclocal}, or give
8621 aclocal the @samp{-I} option when running it.
8624 * Configure script options:: Configuring a package that uses AM_PATH_GINAC.
8625 * Example package:: Example of a package using AM_PATH_GINAC.
8629 @node Configure script options, Example package, AM_PATH_GINAC, AM_PATH_GINAC
8630 @c node-name, next, previous, up
8631 @subsection Configuring a package that uses @samp{AM_PATH_GINAC}
8633 Simply make sure that @command{ginac-config} is in your path, and run
8634 the configure script.
8641 The directory where the GiNaC libraries are installed needs
8642 to be found by your system's dynamic linker.
8644 This is generally done by
8647 editing @file{/etc/ld.so.conf} and running @command{ldconfig}
8653 setting the environment variable @env{LD_LIBRARY_PATH},
8656 or, as a last resort,
8659 giving a @samp{-R} or @samp{-rpath} flag (depending on your linker) when
8660 running configure, for instance:
8663 LDFLAGS=-R/home/cbauer/lib ./configure
8668 You can also specify a @command{ginac-config} not in your path by
8669 setting the @env{GINACLIB_CONFIG} environment variable to the
8670 name of the executable
8673 If you move the GiNaC package from its installed location,
8674 you will either need to modify @command{ginac-config} script
8675 manually to point to the new location or rebuild GiNaC.
8686 --with-ginac-prefix=@var{PREFIX}
8687 --with-ginac-exec-prefix=@var{PREFIX}
8690 are provided to override the prefix and exec-prefix that were stored
8691 in the @command{ginac-config} shell script by GiNaC's configure. You are
8692 generally better off configuring GiNaC with the right path to begin with.
8696 @node Example package, Bibliography, Configure script options, AM_PATH_GINAC
8697 @c node-name, next, previous, up
8698 @subsection Example of a package using @samp{AM_PATH_GINAC}
8700 The following shows how to build a simple package using automake
8701 and the @samp{AM_PATH_GINAC} macro. The program used here is @file{simple.cpp}:
8705 #include <ginac/ginac.h>
8709 GiNaC::symbol x("x");
8710 GiNaC::ex a = GiNaC::sin(x);
8711 std::cout << "Derivative of " << a
8712 << " is " << a.diff(x) << std::endl;
8717 You should first read the introductory portions of the automake
8718 Manual, if you are not already familiar with it.
8720 Two files are needed, @file{configure.in}, which is used to build the
8724 dnl Process this file with autoconf to produce a configure script.
8726 AM_INIT_AUTOMAKE(simple.cpp, 1.0.0)
8732 AM_PATH_GINAC(0.9.0, [
8733 LIBS="$LIBS $GINACLIB_LIBS"
8734 CPPFLAGS="$CPPFLAGS $GINACLIB_CPPFLAGS"
8735 ], AC_MSG_ERROR([need to have GiNaC installed]))
8740 The only command in this which is not standard for automake
8741 is the @samp{AM_PATH_GINAC} macro.
8743 That command does the following: If a GiNaC version greater or equal
8744 than 0.7.0 is found, then it adds @env{$GINACLIB_LIBS} to @env{$LIBS}
8745 and @env{$GINACLIB_CPPFLAGS} to @env{$CPPFLAGS}. Otherwise, it dies with
8746 the error message `need to have GiNaC installed'
8748 And the @file{Makefile.am}, which will be used to build the Makefile.
8751 ## Process this file with automake to produce Makefile.in
8752 bin_PROGRAMS = simple
8753 simple_SOURCES = simple.cpp
8756 This @file{Makefile.am}, says that we are building a single executable,
8757 from a single source file @file{simple.cpp}. Since every program
8758 we are building uses GiNaC we simply added the GiNaC options
8759 to @env{$LIBS} and @env{$CPPFLAGS}, but in other circumstances, we might
8760 want to specify them on a per-program basis: for instance by
8764 simple_LDADD = $(GINACLIB_LIBS)
8765 INCLUDES = $(GINACLIB_CPPFLAGS)
8768 to the @file{Makefile.am}.
8770 To try this example out, create a new directory and add the three
8773 Now execute the following commands:
8776 $ automake --add-missing
8781 You now have a package that can be built in the normal fashion
8790 @node Bibliography, Concept index, Example package, Top
8791 @c node-name, next, previous, up
8792 @appendix Bibliography
8797 @cite{ISO/IEC 14882:1998: Programming Languages: C++}
8800 @cite{CLN: A Class Library for Numbers}, @email{haible@@ilog.fr, Bruno Haible}
8803 @cite{The C++ Programming Language}, Bjarne Stroustrup, 3rd Edition, ISBN 0-201-88954-4, Addison Wesley
8806 @cite{C++ FAQs}, Marshall Cline, ISBN 0-201-58958-3, 1995, Addison Wesley
8809 @cite{Algorithms for Computer Algebra}, Keith O. Geddes, Stephen R. Czapor,
8810 and George Labahn, ISBN 0-7923-9259-0, 1992, Kluwer Academic Publishers, Norwell, Massachusetts
8813 @cite{Computer Algebra: Systems and Algorithms for Algebraic Computation},
8814 James H. Davenport, Yvon Siret and Evelyne Tournier, ISBN 0-12-204230-1, 1988,
8815 Academic Press, London
8818 @cite{Computer Algebra Systems - A Practical Guide},
8819 Michael J. Wester (editor), ISBN 0-471-98353-5, 1999, Wiley, Chichester
8822 @cite{The Art of Computer Programming, Vol 2: Seminumerical Algorithms},
8823 Donald E. Knuth, ISBN 0-201-89684-2, 1998, Addison Wesley
8826 @cite{Pi Unleashed}, J@"org Arndt and Christoph Haenel,
8827 ISBN 3-540-66572-2, 2001, Springer, Heidelberg
8830 @cite{The Role of gamma5 in Dimensional Regularization}, Dirk Kreimer, hep-ph/9401354
8835 @node Concept index, , Bibliography, Top
8836 @c node-name, next, previous, up
8837 @unnumbered Concept index