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-2002 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 The GiNaC Group:
51 @author Christian Bauer, Alexander Frink, Richard Kreckel
54 @vskip 0pt plus 1filll
55 Copyright @copyright{} 1999-2002 Johannes Gutenberg University Mainz, Germany
57 Permission is granted to make and distribute verbatim copies of
58 this manual provided the copyright notice and this permission notice
59 are preserved on all copies.
61 Permission is granted to copy and distribute modified versions of this
62 manual under the conditions for verbatim copying, provided that the entire
63 resulting derived work is distributed under the terms of a permission
64 notice identical to this one.
73 @node Top, Introduction, (dir), (dir)
74 @c node-name, next, previous, up
77 This is a tutorial that documents GiNaC @value{VERSION}, an open
78 framework for symbolic computation within the C++ programming language.
81 * Introduction:: GiNaC's purpose.
82 * A Tour of GiNaC:: A quick tour of the library.
83 * Installation:: How to install the package.
84 * Basic Concepts:: Description of fundamental classes.
85 * Methods and Functions:: Algorithms for symbolic manipulations.
86 * Extending GiNaC:: How to extend the library.
87 * A Comparison With Other CAS:: Compares GiNaC to traditional CAS.
88 * Internal Structures:: Description of some internal structures.
89 * Package Tools:: Configuring packages to work with GiNaC.
95 @node Introduction, A Tour of GiNaC, Top, Top
96 @c node-name, next, previous, up
98 @cindex history of GiNaC
100 The motivation behind GiNaC derives from the observation that most
101 present day computer algebra systems (CAS) are linguistically and
102 semantically impoverished. Although they are quite powerful tools for
103 learning math and solving particular problems they lack modern
104 linguistic structures that allow for the creation of large-scale
105 projects. GiNaC is an attempt to overcome this situation by extending a
106 well established and standardized computer language (C++) by some
107 fundamental symbolic capabilities, thus allowing for integrated systems
108 that embed symbolic manipulations together with more established areas
109 of computer science (like computation-intense numeric applications,
110 graphical interfaces, etc.) under one roof.
112 The particular problem that led to the writing of the GiNaC framework is
113 still a very active field of research, namely the calculation of higher
114 order corrections to elementary particle interactions. There,
115 theoretical physicists are interested in matching present day theories
116 against experiments taking place at particle accelerators. The
117 computations involved are so complex they call for a combined symbolical
118 and numerical approach. This turned out to be quite difficult to
119 accomplish with the present day CAS we have worked with so far and so we
120 tried to fill the gap by writing GiNaC. But of course its applications
121 are in no way restricted to theoretical physics.
123 This tutorial is intended for the novice user who is new to GiNaC but
124 already has some background in C++ programming. However, since a
125 hand-made documentation like this one is difficult to keep in sync with
126 the development, the actual documentation is inside the sources in the
127 form of comments. That documentation may be parsed by one of the many
128 Javadoc-like documentation systems. If you fail at generating it you
129 may access it from @uref{http://www.ginac.de/reference/, the GiNaC home
130 page}. It is an invaluable resource not only for the advanced user who
131 wishes to extend the system (or chase bugs) but for everybody who wants
132 to comprehend the inner workings of GiNaC. This little tutorial on the
133 other hand only covers the basic things that are unlikely to change in
137 The GiNaC framework for symbolic computation within the C++ programming
138 language is Copyright @copyright{} 1999-2002 Johannes Gutenberg
139 University Mainz, Germany.
141 This program is free software; you can redistribute it and/or
142 modify it under the terms of the GNU General Public License as
143 published by the Free Software Foundation; either version 2 of the
144 License, or (at your option) any later version.
146 This program is distributed in the hope that it will be useful, but
147 WITHOUT ANY WARRANTY; without even the implied warranty of
148 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
149 General Public License for more details.
151 You should have received a copy of the GNU General Public License
152 along with this program; see the file COPYING. If not, write to the
153 Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
157 @node A Tour of GiNaC, How to use it from within C++, Introduction, Top
158 @c node-name, next, previous, up
159 @chapter A Tour of GiNaC
161 This quick tour of GiNaC wants to arise your interest in the
162 subsequent chapters by showing off a bit. Please excuse us if it
163 leaves many open questions.
166 * How to use it from within C++:: Two simple examples.
167 * What it can do for you:: A Tour of GiNaC's features.
171 @node How to use it from within C++, What it can do for you, A Tour of GiNaC, A Tour of GiNaC
172 @c node-name, next, previous, up
173 @section How to use it from within C++
175 The GiNaC open framework for symbolic computation within the C++ programming
176 language does not try to define a language of its own as conventional
177 CAS do. Instead, it extends the capabilities of C++ by symbolic
178 manipulations. Here is how to generate and print a simple (and rather
179 pointless) bivariate polynomial with some large coefficients:
183 #include <ginac/ginac.h>
185 using namespace GiNaC;
189 symbol x("x"), y("y");
192 for (int i=0; i<3; ++i)
193 poly += factorial(i+16)*pow(x,i)*pow(y,2-i);
195 cout << poly << endl;
200 Assuming the file is called @file{hello.cc}, on our system we can compile
201 and run it like this:
204 $ c++ hello.cc -o hello -lcln -lginac
206 355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
209 (@xref{Package Tools}, for tools that help you when creating a software
210 package that uses GiNaC.)
212 @cindex Hermite polynomial
213 Next, there is a more meaningful C++ program that calls a function which
214 generates Hermite polynomials in a specified free variable.
218 #include <ginac/ginac.h>
220 using namespace GiNaC;
222 ex HermitePoly(const symbol & x, int n)
224 ex HKer=exp(-pow(x, 2));
225 // uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2)
226 return normal(pow(-1, n) * diff(HKer, x, n) / HKer);
233 for (int i=0; i<6; ++i)
234 cout << "H_" << i << "(z) == " << HermitePoly(z,i) << endl;
240 When run, this will type out
246 H_3(z) == -12*z+8*z^3
247 H_4(z) == -48*z^2+16*z^4+12
248 H_5(z) == 120*z-160*z^3+32*z^5
251 This method of generating the coefficients is of course far from optimal
252 for production purposes.
254 In order to show some more examples of what GiNaC can do we will now use
255 the @command{ginsh}, a simple GiNaC interactive shell that provides a
256 convenient window into GiNaC's capabilities.
259 @node What it can do for you, Installation, How to use it from within C++, A Tour of GiNaC
260 @c node-name, next, previous, up
261 @section What it can do for you
263 @cindex @command{ginsh}
264 After invoking @command{ginsh} one can test and experiment with GiNaC's
265 features much like in other Computer Algebra Systems except that it does
266 not provide programming constructs like loops or conditionals. For a
267 concise description of the @command{ginsh} syntax we refer to its
268 accompanied man page. Suffice to say that assignments and comparisons in
269 @command{ginsh} are written as they are in C, i.e. @code{=} assigns and
272 It can manipulate arbitrary precision integers in a very fast way.
273 Rational numbers are automatically converted to fractions of coprime
278 369988485035126972924700782451696644186473100389722973815184405301748249
280 123329495011708990974900260817232214728824366796574324605061468433916083
287 Exact numbers are always retained as exact numbers and only evaluated as
288 floating point numbers if requested. For instance, with numeric
289 radicals is dealt pretty much as with symbols. Products of sums of them
293 > expand((1+a^(1/5)-a^(2/5))^3);
294 1+3*a+3*a^(1/5)-5*a^(3/5)-a^(6/5)
295 > expand((1+3^(1/5)-3^(2/5))^3);
297 > evalf((1+3^(1/5)-3^(2/5))^3);
298 0.33408977534118624228
301 The function @code{evalf} that was used above converts any number in
302 GiNaC's expressions into floating point numbers. This can be done to
303 arbitrary predefined accuracy:
307 0.14285714285714285714
311 0.1428571428571428571428571428571428571428571428571428571428571428571428
312 5714285714285714285714285714285714285
315 Exact numbers other than rationals that can be manipulated in GiNaC
316 include predefined constants like Archimedes' @code{Pi}. They can both
317 be used in symbolic manipulations (as an exact number) as well as in
318 numeric expressions (as an inexact number):
324 9.869604401089358619+x
328 11.869604401089358619
331 Built-in functions evaluate immediately to exact numbers if
332 this is possible. Conversions that can be safely performed are done
333 immediately; conversions that are not generally valid are not done:
344 (Note that converting the last input to @code{x} would allow one to
345 conclude that @code{42*Pi} is equal to @code{0}.)
347 Linear equation systems can be solved along with basic linear
348 algebra manipulations over symbolic expressions. In C++ GiNaC offers
349 a matrix class for this purpose but we can see what it can do using
350 @command{ginsh}'s bracket notation to type them in:
353 > lsolve(a+x*y==z,x);
355 > lsolve(@{3*x+5*y == 7, -2*x+10*y == -5@}, @{x, y@});
357 > M = [ [1, 3], [-3, 2] ];
361 > charpoly(M,lambda);
363 > A = [ [1, 1], [2, -1] ];
366 [[1,1],[2,-1]]+2*[[1,3],[-3,2]]
369 > B = [ [0, 0, a], [b, 1, -b], [-1/a, 0, 0] ];
370 > evalm(B^(2^12345));
371 [[1,0,0],[0,1,0],[0,0,1]]
374 Multivariate polynomials and rational functions may be expanded,
375 collected and normalized (i.e. converted to a ratio of two coprime
379 > a = x^4 + 2*x^2*y^2 + 4*x^3*y + 12*x*y^3 - 3*y^4;
380 12*x*y^3+2*x^2*y^2+4*x^3*y-3*y^4+x^4
381 > b = x^2 + 4*x*y - y^2;
384 8*x^5*y+17*x^4*y^2+43*x^2*y^4-24*x*y^5+16*x^3*y^3+3*y^6+x^6
386 4*x^3*y-y^2-3*y^4+(12*y^3+4*y)*x+x^4+x^2*(1+2*y^2)
388 12*x*y^3-3*y^4+(-1+2*x^2)*y^2+(4*x+4*x^3)*y+x^2+x^4
393 You can differentiate functions and expand them as Taylor or Laurent
394 series in a very natural syntax (the second argument of @code{series} is
395 a relation defining the evaluation point, the third specifies the
398 @cindex Zeta function
402 > series(sin(x),x==0,4);
404 > series(1/tan(x),x==0,4);
405 x^(-1)-1/3*x+Order(x^2)
406 > series(tgamma(x),x==0,3);
407 x^(-1)-Euler+(1/12*Pi^2+1/2*Euler^2)*x+
408 (-1/3*zeta(3)-1/12*Pi^2*Euler-1/6*Euler^3)*x^2+Order(x^3)
410 x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
411 -(0.90747907608088628905)*x^2+Order(x^3)
412 > series(tgamma(2*sin(x)-2),x==Pi/2,6);
413 -(x-1/2*Pi)^(-2)+(-1/12*Pi^2-1/2*Euler^2-1/240)*(x-1/2*Pi)^2
414 -Euler-1/12+Order((x-1/2*Pi)^3)
417 Here we have made use of the @command{ginsh}-command @code{%} to pop the
418 previously evaluated element from @command{ginsh}'s internal stack.
420 If you ever wanted to convert units in C or C++ and found this is
421 cumbersome, here is the solution. Symbolic types can always be used as
422 tags for different types of objects. Converting from wrong units to the
423 metric system is now easy:
431 140613.91592783185568*kg*m^(-2)
435 @node Installation, Prerequisites, What it can do for you, Top
436 @c node-name, next, previous, up
437 @chapter Installation
440 GiNaC's installation follows the spirit of most GNU software. It is
441 easily installed on your system by three steps: configuration, build,
445 * Prerequisites:: Packages upon which GiNaC depends.
446 * Configuration:: How to configure GiNaC.
447 * Building GiNaC:: How to compile GiNaC.
448 * Installing GiNaC:: How to install GiNaC on your system.
452 @node Prerequisites, Configuration, Installation, Installation
453 @c node-name, next, previous, up
454 @section Prerequisites
456 In order to install GiNaC on your system, some prerequisites need to be
457 met. First of all, you need to have a C++-compiler adhering to the
458 ANSI-standard @cite{ISO/IEC 14882:1998(E)}. We used GCC for development
459 so if you have a different compiler you are on your own. For the
460 configuration to succeed you need a Posix compliant shell installed in
461 @file{/bin/sh}, GNU @command{bash} is fine. Perl is needed by the built
462 process as well, since some of the source files are automatically
463 generated by Perl scripts. Last but not least, Bruno Haible's library
464 CLN is extensively used and needs to be installed on your system.
465 Please get it either from @uref{ftp://ftp.santafe.edu/pub/gnu/}, from
466 @uref{ftp://ftpthep.physik.uni-mainz.de/pub/gnu/, GiNaC's FTP site} or
467 from @uref{ftp://ftp.ilog.fr/pub/Users/haible/gnu/, Bruno Haible's FTP
468 site} (it is covered by GPL) and install it prior to trying to install
469 GiNaC. The configure script checks if it can find it and if it cannot
470 it will refuse to continue.
473 @node Configuration, Building GiNaC, Prerequisites, Installation
474 @c node-name, next, previous, up
475 @section Configuration
476 @cindex configuration
479 To configure GiNaC means to prepare the source distribution for
480 building. It is done via a shell script called @command{configure} that
481 is shipped with the sources and was originally generated by GNU
482 Autoconf. Since a configure script generated by GNU Autoconf never
483 prompts, all customization must be done either via command line
484 parameters or environment variables. It accepts a list of parameters,
485 the complete set of which can be listed by calling it with the
486 @option{--help} option. The most important ones will be shortly
487 described in what follows:
492 @option{--disable-shared}: When given, this option switches off the
493 build of a shared library, i.e. a @file{.so} file. This may be convenient
494 when developing because it considerably speeds up compilation.
497 @option{--prefix=@var{PREFIX}}: The directory where the compiled library
498 and headers are installed. It defaults to @file{/usr/local} which means
499 that the library is installed in the directory @file{/usr/local/lib},
500 the header files in @file{/usr/local/include/ginac} and the documentation
501 (like this one) into @file{/usr/local/share/doc/GiNaC}.
504 @option{--libdir=@var{LIBDIR}}: Use this option in case you want to have
505 the library installed in some other directory than
506 @file{@var{PREFIX}/lib/}.
509 @option{--includedir=@var{INCLUDEDIR}}: Use this option in case you want
510 to have the header files installed in some other directory than
511 @file{@var{PREFIX}/include/ginac/}. For instance, if you specify
512 @option{--includedir=/usr/include} you will end up with the header files
513 sitting in the directory @file{/usr/include/ginac/}. Note that the
514 subdirectory @file{ginac} is enforced by this process in order to
515 keep the header files separated from others. This avoids some
516 clashes and allows for an easier deinstallation of GiNaC. This ought
517 to be considered A Good Thing (tm).
520 @option{--datadir=@var{DATADIR}}: This option may be given in case you
521 want to have the documentation installed in some other directory than
522 @file{@var{PREFIX}/share/doc/GiNaC/}.
526 In addition, you may specify some environment variables. @env{CXX}
527 holds the path and the name of the C++ compiler in case you want to
528 override the default in your path. (The @command{configure} script
529 searches your path for @command{c++}, @command{g++}, @command{gcc},
530 @command{CC}, @command{cxx} and @command{cc++} in that order.) It may
531 be very useful to define some compiler flags with the @env{CXXFLAGS}
532 environment variable, like optimization, debugging information and
533 warning levels. If omitted, it defaults to @option{-g
534 -O2}.@footnote{The @command{configure} script is itself generated from
535 the file @file{configure.ac}. It is only distributed in packaged
536 releases of GiNaC. If you got the naked sources, e.g. from CVS, you
537 must generate @command{configure} along with the various
538 @file{Makefile.in} by using the @command{autogen.sh} script. This will
539 require a fair amount of support from your local toolchain, though.}
541 The whole process is illustrated in the following two
542 examples. (Substitute @command{setenv @var{VARIABLE} @var{value}} for
543 @command{export @var{VARIABLE}=@var{value}} if the Berkeley C shell is
546 Here is a simple configuration for a site-wide GiNaC library assuming
547 everything is in default paths:
550 $ export CXXFLAGS="-Wall -O2"
554 And here is a configuration for a private static GiNaC library with
555 several components sitting in custom places (site-wide GCC and private
556 CLN). The compiler is persuaded to be picky and full assertions and
557 debugging information are switched on:
560 $ export CXX=/usr/local/gnu/bin/c++
561 $ export CPPFLAGS="$(CPPFLAGS) -I$(HOME)/include"
562 $ export CXXFLAGS="$(CXXFLAGS) -DDO_GINAC_ASSERT -ggdb -Wall -pedantic"
563 $ export LDFLAGS="$(LDFLAGS) -L$(HOME)/lib"
564 $ ./configure --disable-shared --prefix=$(HOME)
568 @node Building GiNaC, Installing GiNaC, Configuration, Installation
569 @c node-name, next, previous, up
570 @section Building GiNaC
571 @cindex building GiNaC
573 After proper configuration you should just build the whole
578 at the command prompt and go for a cup of coffee. The exact time it
579 takes to compile GiNaC depends not only on the speed of your machines
580 but also on other parameters, for instance what value for @env{CXXFLAGS}
581 you entered. Optimization may be very time-consuming.
583 Just to make sure GiNaC works properly you may run a collection of
584 regression tests by typing
590 This will compile some sample programs, run them and check the output
591 for correctness. The regression tests fall in three categories. First,
592 the so called @emph{exams} are performed, simple tests where some
593 predefined input is evaluated (like a pupils' exam). Second, the
594 @emph{checks} test the coherence of results among each other with
595 possible random input. Third, some @emph{timings} are performed, which
596 benchmark some predefined problems with different sizes and display the
597 CPU time used in seconds. Each individual test should return a message
598 @samp{passed}. This is mostly intended to be a QA-check if something
599 was broken during development, not a sanity check of your system. Some
600 of the tests in sections @emph{checks} and @emph{timings} may require
601 insane amounts of memory and CPU time. Feel free to kill them if your
602 machine catches fire. Another quite important intent is to allow people
603 to fiddle around with optimization.
605 Generally, the top-level Makefile runs recursively to the
606 subdirectories. It is therefore safe to go into any subdirectory
607 (@code{doc/}, @code{ginsh/}, @dots{}) and simply type @code{make}
608 @var{target} there in case something went wrong.
611 @node Installing GiNaC, Basic Concepts, Building GiNaC, Installation
612 @c node-name, next, previous, up
613 @section Installing GiNaC
616 To install GiNaC on your system, simply type
622 As described in the section about configuration the files will be
623 installed in the following directories (the directories will be created
624 if they don't already exist):
629 @file{libginac.a} will go into @file{@var{PREFIX}/lib/} (or
630 @file{@var{LIBDIR}}) which defaults to @file{/usr/local/lib/}.
631 So will @file{libginac.so} unless the configure script was
632 given the option @option{--disable-shared}. The proper symlinks
633 will be established as well.
636 All the header files will be installed into @file{@var{PREFIX}/include/ginac/}
637 (or @file{@var{INCLUDEDIR}/ginac/}, if specified).
640 All documentation (HTML and Postscript) will be stuffed into
641 @file{@var{PREFIX}/share/doc/GiNaC/} (or
642 @file{@var{DATADIR}/doc/GiNaC/}, if @var{DATADIR} was specified).
646 For the sake of completeness we will list some other useful make
647 targets: @command{make clean} deletes all files generated by
648 @command{make}, i.e. all the object files. In addition @command{make
649 distclean} removes all files generated by the configuration and
650 @command{make maintainer-clean} goes one step further and deletes files
651 that may require special tools to rebuild (like the @command{libtool}
652 for instance). Finally @command{make uninstall} removes the installed
653 library, header files and documentation@footnote{Uninstallation does not
654 work after you have called @command{make distclean} since the
655 @file{Makefile} is itself generated by the configuration from
656 @file{Makefile.in} and hence deleted by @command{make distclean}. There
657 are two obvious ways out of this dilemma. First, you can run the
658 configuration again with the same @var{PREFIX} thus creating a
659 @file{Makefile} with a working @samp{uninstall} target. Second, you can
660 do it by hand since you now know where all the files went during
664 @node Basic Concepts, Expressions, Installing GiNaC, Top
665 @c node-name, next, previous, up
666 @chapter Basic Concepts
668 This chapter will describe the different fundamental objects that can be
669 handled by GiNaC. But before doing so, it is worthwhile introducing you
670 to the more commonly used class of expressions, representing a flexible
671 meta-class for storing all mathematical objects.
674 * Expressions:: The fundamental GiNaC class.
675 * The Class Hierarchy:: Overview of GiNaC's classes.
676 * Error handling:: How the library reports errors.
677 * Symbols:: Symbolic objects.
678 * Numbers:: Numerical objects.
679 * Constants:: Pre-defined constants.
680 * Fundamental containers:: The power, add and mul classes.
681 * Lists:: Lists of expressions.
682 * Mathematical functions:: Mathematical functions.
683 * Relations:: Equality, Inequality and all that.
684 * Matrices:: Matrices.
685 * Indexed objects:: Handling indexed quantities.
686 * Non-commutative objects:: Algebras with non-commutative products.
690 @node Expressions, The Class Hierarchy, Basic Concepts, Basic Concepts
691 @c node-name, next, previous, up
693 @cindex expression (class @code{ex})
696 The most common class of objects a user deals with is the expression
697 @code{ex}, representing a mathematical object like a variable, number,
698 function, sum, product, etc@dots{} Expressions may be put together to form
699 new expressions, passed as arguments to functions, and so on. Here is a
700 little collection of valid expressions:
703 ex MyEx1 = 5; // simple number
704 ex MyEx2 = x + 2*y; // polynomial in x and y
705 ex MyEx3 = (x + 1)/(x - 1); // rational expression
706 ex MyEx4 = sin(x + 2*y) + 3*z + 41; // containing a function
707 ex MyEx5 = MyEx4 + 1; // similar to above
710 Expressions are handles to other more fundamental objects, that often
711 contain other expressions thus creating a tree of expressions
712 (@xref{Internal Structures}, for particular examples). Most methods on
713 @code{ex} therefore run top-down through such an expression tree. For
714 example, the method @code{has()} scans recursively for occurrences of
715 something inside an expression. Thus, if you have declared @code{MyEx4}
716 as in the example above @code{MyEx4.has(y)} will find @code{y} inside
717 the argument of @code{sin} and hence return @code{true}.
719 The next sections will outline the general picture of GiNaC's class
720 hierarchy and describe the classes of objects that are handled by
724 @node The Class Hierarchy, Error handling, Expressions, Basic Concepts
725 @c node-name, next, previous, up
726 @section The Class Hierarchy
728 GiNaC's class hierarchy consists of several classes representing
729 mathematical objects, all of which (except for @code{ex} and some
730 helpers) are internally derived from one abstract base class called
731 @code{basic}. You do not have to deal with objects of class
732 @code{basic}, instead you'll be dealing with symbols, numbers,
733 containers of expressions and so on.
737 To get an idea about what kinds of symbolic composits may be built we
738 have a look at the most important classes in the class hierarchy and
739 some of the relations among the classes:
741 @image{classhierarchy}
743 The abstract classes shown here (the ones without drop-shadow) are of no
744 interest for the user. They are used internally in order to avoid code
745 duplication if two or more classes derived from them share certain
746 features. An example is @code{expairseq}, a container for a sequence of
747 pairs each consisting of one expression and a number (@code{numeric}).
748 What @emph{is} visible to the user are the derived classes @code{add}
749 and @code{mul}, representing sums and products. @xref{Internal
750 Structures}, where these two classes are described in more detail. The
751 following table shortly summarizes what kinds of mathematical objects
752 are stored in the different classes:
755 @multitable @columnfractions .22 .78
756 @item @code{symbol} @tab Algebraic symbols @math{a}, @math{x}, @math{y}@dots{}
757 @item @code{constant} @tab Constants like
764 @item @code{numeric} @tab All kinds of numbers, @math{42}, @math{7/3*I}, @math{3.14159}@dots{}
765 @item @code{add} @tab Sums like @math{x+y} or @math{a-(2*b)+3}
766 @item @code{mul} @tab Products like @math{x*y} or @math{2*a^2*(x+y+z)/b}
767 @item @code{ncmul} @tab Products of non-commutative objects
768 @item @code{power} @tab Exponentials such as @math{x^2}, @math{a^b},
773 @code{sqrt(}@math{2}@code{)}
776 @item @code{pseries} @tab Power Series, e.g. @math{x-1/6*x^3+1/120*x^5+O(x^7)}
777 @item @code{function} @tab A symbolic function like @math{sin(2*x)}
778 @item @code{lst} @tab Lists of expressions @{@math{x}, @math{2*y}, @math{3+z}@}
779 @item @code{matrix} @tab @math{m}x@math{n} matrices of expressions
780 @item @code{relational} @tab A relation like the identity @math{x}@code{==}@math{y}
781 @item @code{indexed} @tab Indexed object like @math{A_ij}
782 @item @code{tensor} @tab Special tensor like the delta and metric tensors
783 @item @code{idx} @tab Index of an indexed object
784 @item @code{varidx} @tab Index with variance
785 @item @code{spinidx} @tab Index with variance and dot (used in Weyl-van-der-Waerden spinor formalism)
786 @item @code{wildcard} @tab Wildcard for pattern matching
791 @node Error handling, Symbols, The Class Hierarchy, Basic Concepts
792 @c node-name, next, previous, up
793 @section Error handling
795 @cindex @code{pole_error} (class)
797 GiNaC reports run-time errors by throwing C++ exceptions. All exceptions
798 generated by GiNaC are subclassed from the standard @code{exception} class
799 defined in the @file{<stdexcept>} header. In addition to the predefined
800 @code{logic_error}, @code{domain_error}, @code{out_of_range},
801 @code{invalid_argument}, @code{runtime_error}, @code{range_error} and
802 @code{overflow_error} types, GiNaC also defines a @code{pole_error}
803 exception that gets thrown when trying to evaluate a mathematical function
806 The @code{pole_error} class has a member function
809 int pole_error::degree(void) const;
812 that returns the order of the singularity (or 0 when the pole is
813 logarithmic or the order is undefined).
815 When using GiNaC it is useful to arrange for exceptions to be catched in
816 the main program even if you don't want to do any special error handling.
817 Otherwise whenever an error occurs in GiNaC, it will be delegated to the
818 default exception handler of your C++ compiler's run-time system which
819 usually only aborts the program without giving any information what went
822 Here is an example for a @code{main()} function that catches and prints
823 exceptions generated by GiNaC:
828 #include <ginac/ginac.h>
830 using namespace GiNaC;
838 @} catch (exception &p) @{
839 cerr << p.what() << endl;
847 @node Symbols, Numbers, Error handling, Basic Concepts
848 @c node-name, next, previous, up
850 @cindex @code{symbol} (class)
851 @cindex hierarchy of classes
854 Symbols are for symbolic manipulation what atoms are for chemistry. You
855 can declare objects of class @code{symbol} as any other object simply by
856 saying @code{symbol x,y;}. There is, however, a catch in here having to
857 do with the fact that C++ is a compiled language. The information about
858 the symbol's name is thrown away by the compiler but at a later stage
859 you may want to print expressions holding your symbols. In order to
860 avoid confusion GiNaC's symbols are able to know their own name. This
861 is accomplished by declaring its name for output at construction time in
862 the fashion @code{symbol x("x");}. If you declare a symbol using the
863 default constructor (i.e. without string argument) the system will deal
864 out a unique name. That name may not be suitable for printing but for
865 internal routines when no output is desired it is often enough. We'll
866 come across examples of such symbols later in this tutorial.
868 This implies that the strings passed to symbols at construction time may
869 not be used for comparing two of them. It is perfectly legitimate to
870 write @code{symbol x("x"),y("x");} but it is likely to lead into
871 trouble. Here, @code{x} and @code{y} are different symbols and
872 statements like @code{x-y} will not be simplified to zero although the
873 output @code{x-x} looks funny. Such output may also occur when there
874 are two different symbols in two scopes, for instance when you call a
875 function that declares a symbol with a name already existent in a symbol
876 in the calling function. Again, comparing them (using @code{operator==}
877 for instance) will always reveal their difference. Watch out, please.
879 @cindex @code{subs()}
880 Although symbols can be assigned expressions for internal reasons, you
881 should not do it (and we are not going to tell you how it is done). If
882 you want to replace a symbol with something else in an expression, you
883 can use the expression's @code{.subs()} method (@pxref{Substituting Expressions}).
886 @node Numbers, Constants, Symbols, Basic Concepts
887 @c node-name, next, previous, up
889 @cindex @code{numeric} (class)
895 For storing numerical things, GiNaC uses Bruno Haible's library CLN.
896 The classes therein serve as foundation classes for GiNaC. CLN stands
897 for Class Library for Numbers or alternatively for Common Lisp Numbers.
898 In order to find out more about CLN's internals the reader is refered to
899 the documentation of that library. @inforef{Introduction, , cln}, for
900 more information. Suffice to say that it is by itself build on top of
901 another library, the GNU Multiple Precision library GMP, which is an
902 extremely fast library for arbitrary long integers and rationals as well
903 as arbitrary precision floating point numbers. It is very commonly used
904 by several popular cryptographic applications. CLN extends GMP by
905 several useful things: First, it introduces the complex number field
906 over either reals (i.e. floating point numbers with arbitrary precision)
907 or rationals. Second, it automatically converts rationals to integers
908 if the denominator is unity and complex numbers to real numbers if the
909 imaginary part vanishes and also correctly treats algebraic functions.
910 Third it provides good implementations of state-of-the-art algorithms
911 for all trigonometric and hyperbolic functions as well as for
912 calculation of some useful constants.
914 The user can construct an object of class @code{numeric} in several
915 ways. The following example shows the four most important constructors.
916 It uses construction from C-integer, construction of fractions from two
917 integers, construction from C-float and construction from a string:
921 #include <ginac/ginac.h>
922 using namespace GiNaC;
926 numeric two = 2; // exact integer 2
927 numeric r(2,3); // exact fraction 2/3
928 numeric e(2.71828); // floating point number
929 numeric p = "3.14159265358979323846"; // constructor from string
930 // Trott's constant in scientific notation:
931 numeric trott("1.0841015122311136151E-2");
933 std::cout << two*p << std::endl; // floating point 6.283...
938 @cindex complex numbers
939 The imaginary unit in GiNaC is a predefined @code{numeric} object with the
944 numeric z1 = 2-3*I; // exact complex number 2-3i
945 numeric z2 = 5.9+1.6*I; // complex floating point number
949 It may be tempting to construct fractions by writing @code{numeric r(3/2)}.
950 This would, however, call C's built-in operator @code{/} for integers
951 first and result in a numeric holding a plain integer 1. @strong{Never
952 use the operator @code{/} on integers} unless you know exactly what you
953 are doing! Use the constructor from two integers instead, as shown in
954 the example above. Writing @code{numeric(1)/2} may look funny but works
957 @cindex @code{Digits}
959 We have seen now the distinction between exact numbers and floating
960 point numbers. Clearly, the user should never have to worry about
961 dynamically created exact numbers, since their `exactness' always
962 determines how they ought to be handled, i.e. how `long' they are. The
963 situation is different for floating point numbers. Their accuracy is
964 controlled by one @emph{global} variable, called @code{Digits}. (For
965 those readers who know about Maple: it behaves very much like Maple's
966 @code{Digits}). All objects of class numeric that are constructed from
967 then on will be stored with a precision matching that number of decimal
972 #include <ginac/ginac.h>
974 using namespace GiNaC;
978 numeric three(3.0), one(1.0);
979 numeric x = one/three;
981 cout << "in " << Digits << " digits:" << endl;
983 cout << Pi.evalf() << endl;
995 The above example prints the following output to screen:
999 0.33333333333333333334
1000 3.1415926535897932385
1002 0.33333333333333333333333333333333333333333333333333333333333333333334
1003 3.1415926535897932384626433832795028841971693993751058209749445923078
1007 Note that the last number is not necessarily rounded as you would
1008 naively expect it to be rounded in the decimal system. But note also,
1009 that in both cases you got a couple of extra digits. This is because
1010 numbers are internally stored by CLN as chunks of binary digits in order
1011 to match your machine's word size and to not waste precision. Thus, on
1012 architectures with differnt word size, the above output might even
1013 differ with regard to actually computed digits.
1015 It should be clear that objects of class @code{numeric} should be used
1016 for constructing numbers or for doing arithmetic with them. The objects
1017 one deals with most of the time are the polymorphic expressions @code{ex}.
1019 @subsection Tests on numbers
1021 Once you have declared some numbers, assigned them to expressions and
1022 done some arithmetic with them it is frequently desired to retrieve some
1023 kind of information from them like asking whether that number is
1024 integer, rational, real or complex. For those cases GiNaC provides
1025 several useful methods. (Internally, they fall back to invocations of
1026 certain CLN functions.)
1028 As an example, let's construct some rational number, multiply it with
1029 some multiple of its denominator and test what comes out:
1033 #include <ginac/ginac.h>
1034 using namespace std;
1035 using namespace GiNaC;
1037 // some very important constants:
1038 const numeric twentyone(21);
1039 const numeric ten(10);
1040 const numeric five(5);
1044 numeric answer = twentyone;
1047 cout << answer.is_integer() << endl; // false, it's 21/5
1049 cout << answer.is_integer() << endl; // true, it's 42 now!
1053 Note that the variable @code{answer} is constructed here as an integer
1054 by @code{numeric}'s copy constructor but in an intermediate step it
1055 holds a rational number represented as integer numerator and integer
1056 denominator. When multiplied by 10, the denominator becomes unity and
1057 the result is automatically converted to a pure integer again.
1058 Internally, the underlying CLN is responsible for this behavior and we
1059 refer the reader to CLN's documentation. Suffice to say that
1060 the same behavior applies to complex numbers as well as return values of
1061 certain functions. Complex numbers are automatically converted to real
1062 numbers if the imaginary part becomes zero. The full set of tests that
1063 can be applied is listed in the following table.
1066 @multitable @columnfractions .30 .70
1067 @item @strong{Method} @tab @strong{Returns true if the object is@dots{}}
1068 @item @code{.is_zero()}
1069 @tab @dots{}equal to zero
1070 @item @code{.is_positive()}
1071 @tab @dots{}not complex and greater than 0
1072 @item @code{.is_integer()}
1073 @tab @dots{}a (non-complex) integer
1074 @item @code{.is_pos_integer()}
1075 @tab @dots{}an integer and greater than 0
1076 @item @code{.is_nonneg_integer()}
1077 @tab @dots{}an integer and greater equal 0
1078 @item @code{.is_even()}
1079 @tab @dots{}an even integer
1080 @item @code{.is_odd()}
1081 @tab @dots{}an odd integer
1082 @item @code{.is_prime()}
1083 @tab @dots{}a prime integer (probabilistic primality test)
1084 @item @code{.is_rational()}
1085 @tab @dots{}an exact rational number (integers are rational, too)
1086 @item @code{.is_real()}
1087 @tab @dots{}a real integer, rational or float (i.e. is not complex)
1088 @item @code{.is_cinteger()}
1089 @tab @dots{}a (complex) integer (such as @math{2-3*I})
1090 @item @code{.is_crational()}
1091 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
1096 @node Constants, Fundamental containers, Numbers, Basic Concepts
1097 @c node-name, next, previous, up
1099 @cindex @code{constant} (class)
1102 @cindex @code{Catalan}
1103 @cindex @code{Euler}
1104 @cindex @code{evalf()}
1105 Constants behave pretty much like symbols except that they return some
1106 specific number when the method @code{.evalf()} is called.
1108 The predefined known constants are:
1111 @multitable @columnfractions .14 .30 .56
1112 @item @strong{Name} @tab @strong{Common Name} @tab @strong{Numerical Value (to 35 digits)}
1114 @tab Archimedes' constant
1115 @tab 3.14159265358979323846264338327950288
1116 @item @code{Catalan}
1117 @tab Catalan's constant
1118 @tab 0.91596559417721901505460351493238411
1120 @tab Euler's (or Euler-Mascheroni) constant
1121 @tab 0.57721566490153286060651209008240243
1126 @node Fundamental containers, Lists, Constants, Basic Concepts
1127 @c node-name, next, previous, up
1128 @section Fundamental containers: the @code{power}, @code{add} and @code{mul} classes
1132 @cindex @code{power}
1134 Simple polynomial expressions are written down in GiNaC pretty much like
1135 in other CAS or like expressions involving numerical variables in C.
1136 The necessary operators @code{+}, @code{-}, @code{*} and @code{/} have
1137 been overloaded to achieve this goal. When you run the following
1138 code snippet, the constructor for an object of type @code{mul} is
1139 automatically called to hold the product of @code{a} and @code{b} and
1140 then the constructor for an object of type @code{add} is called to hold
1141 the sum of that @code{mul} object and the number one:
1145 symbol a("a"), b("b");
1150 @cindex @code{pow()}
1151 For exponentiation, you have already seen the somewhat clumsy (though C-ish)
1152 statement @code{pow(x,2);} to represent @code{x} squared. This direct
1153 construction is necessary since we cannot safely overload the constructor
1154 @code{^} in C++ to construct a @code{power} object. If we did, it would
1155 have several counterintuitive and undesired effects:
1159 Due to C's operator precedence, @code{2*x^2} would be parsed as @code{(2*x)^2}.
1161 Due to the binding of the operator @code{^}, @code{x^a^b} would result in
1162 @code{(x^a)^b}. This would be confusing since most (though not all) other CAS
1163 interpret this as @code{x^(a^b)}.
1165 Also, expressions involving integer exponents are very frequently used,
1166 which makes it even more dangerous to overload @code{^} since it is then
1167 hard to distinguish between the semantics as exponentiation and the one
1168 for exclusive or. (It would be embarrassing to return @code{1} where one
1169 has requested @code{2^3}.)
1172 @cindex @command{ginsh}
1173 All effects are contrary to mathematical notation and differ from the
1174 way most other CAS handle exponentiation, therefore overloading @code{^}
1175 is ruled out for GiNaC's C++ part. The situation is different in
1176 @command{ginsh}, there the exponentiation-@code{^} exists. (Also note
1177 that the other frequently used exponentiation operator @code{**} does
1178 not exist at all in C++).
1180 To be somewhat more precise, objects of the three classes described
1181 here, are all containers for other expressions. An object of class
1182 @code{power} is best viewed as a container with two slots, one for the
1183 basis, one for the exponent. All valid GiNaC expressions can be
1184 inserted. However, basic transformations like simplifying
1185 @code{pow(pow(x,2),3)} to @code{x^6} automatically are only performed
1186 when this is mathematically possible. If we replace the outer exponent
1187 three in the example by some symbols @code{a}, the simplification is not
1188 safe and will not be performed, since @code{a} might be @code{1/2} and
1191 Objects of type @code{add} and @code{mul} are containers with an
1192 arbitrary number of slots for expressions to be inserted. Again, simple
1193 and safe simplifications are carried out like transforming
1194 @code{3*x+4-x} to @code{2*x+4}.
1196 The general rule is that when you construct such objects, GiNaC
1197 automatically creates them in canonical form, which might differ from
1198 the form you typed in your program. This allows for rapid comparison of
1199 expressions, since after all @code{a-a} is simply zero. Note, that the
1200 canonical form is not necessarily lexicographical ordering or in any way
1201 easily guessable. It is only guaranteed that constructing the same
1202 expression twice, either implicitly or explicitly, results in the same
1206 @node Lists, Mathematical functions, Fundamental containers, Basic Concepts
1207 @c node-name, next, previous, up
1208 @section Lists of expressions
1209 @cindex @code{lst} (class)
1211 @cindex @code{nops()}
1213 @cindex @code{append()}
1214 @cindex @code{prepend()}
1215 @cindex @code{remove_first()}
1216 @cindex @code{remove_last()}
1218 The GiNaC class @code{lst} serves for holding a @dfn{list} of arbitrary
1219 expressions. These are sometimes used to supply a variable number of
1220 arguments of the same type to GiNaC methods such as @code{subs()} and
1221 @code{to_rational()}, so you should have a basic understanding about them.
1223 Lists of up to 16 expressions can be directly constructed from single
1228 symbol x("x"), y("y");
1229 lst l(x, 2, y, x+y);
1230 // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y'
1234 Use the @code{nops()} method to determine the size (number of expressions) of
1235 a list and the @code{op()} method to access individual elements:
1239 cout << l.nops() << endl; // prints '4'
1240 cout << l.op(2) << " " << l.op(0) << endl; // prints 'y x'
1244 You can append or prepend an expression to a list with the @code{append()}
1245 and @code{prepend()} methods:
1249 l.append(4*x); // l is now @{x, 2, y, x+y, 4*x@}
1250 l.prepend(0); // l is now @{0, x, 2, y, x+y, 4*x@}
1254 Finally you can remove the first or last element of a list with
1255 @code{remove_first()} and @code{remove_last()}:
1259 l.remove_first(); // l is now @{x, 2, y, x+y, 4*x@}
1260 l.remove_last(); // l is now @{x, 2, y, x+y@}
1265 @node Mathematical functions, Relations, Lists, Basic Concepts
1266 @c node-name, next, previous, up
1267 @section Mathematical functions
1268 @cindex @code{function} (class)
1269 @cindex trigonometric function
1270 @cindex hyperbolic function
1272 There are quite a number of useful functions hard-wired into GiNaC. For
1273 instance, all trigonometric and hyperbolic functions are implemented
1274 (@xref{Built-in Functions}, for a complete list).
1276 These functions (better called @emph{pseudofunctions}) are all objects
1277 of class @code{function}. They accept one or more expressions as
1278 arguments and return one expression. If the arguments are not
1279 numerical, the evaluation of the function may be halted, as it does in
1280 the next example, showing how a function returns itself twice and
1281 finally an expression that may be really useful:
1283 @cindex Gamma function
1284 @cindex @code{subs()}
1287 symbol x("x"), y("y");
1289 cout << tgamma(foo) << endl;
1290 // -> tgamma(x+(1/2)*y)
1291 ex bar = foo.subs(y==1);
1292 cout << tgamma(bar) << endl;
1294 ex foobar = bar.subs(x==7);
1295 cout << tgamma(foobar) << endl;
1296 // -> (135135/128)*Pi^(1/2)
1300 Besides evaluation most of these functions allow differentiation, series
1301 expansion and so on. Read the next chapter in order to learn more about
1304 It must be noted that these pseudofunctions are created by inline
1305 functions, where the argument list is templated. This means that
1306 whenever you call @code{GiNaC::sin(1)} it is equivalent to
1307 @code{sin(ex(1))} and will therefore not result in a floating point
1308 number. Unless of course the function prototype is explicitly
1309 overridden -- which is the case for arguments of type @code{numeric}
1310 (not wrapped inside an @code{ex}). Hence, in order to obtain a floating
1311 point number of class @code{numeric} you should call
1312 @code{sin(numeric(1))}. This is almost the same as calling
1313 @code{sin(1).evalf()} except that the latter will return a numeric
1314 wrapped inside an @code{ex}.
1317 @node Relations, Matrices, Mathematical functions, Basic Concepts
1318 @c node-name, next, previous, up
1320 @cindex @code{relational} (class)
1322 Sometimes, a relation holding between two expressions must be stored
1323 somehow. The class @code{relational} is a convenient container for such
1324 purposes. A relation is by definition a container for two @code{ex} and
1325 a relation between them that signals equality, inequality and so on.
1326 They are created by simply using the C++ operators @code{==}, @code{!=},
1327 @code{<}, @code{<=}, @code{>} and @code{>=} between two expressions.
1329 @xref{Mathematical functions}, for examples where various applications
1330 of the @code{.subs()} method show how objects of class relational are
1331 used as arguments. There they provide an intuitive syntax for
1332 substitutions. They are also used as arguments to the @code{ex::series}
1333 method, where the left hand side of the relation specifies the variable
1334 to expand in and the right hand side the expansion point. They can also
1335 be used for creating systems of equations that are to be solved for
1336 unknown variables. But the most common usage of objects of this class
1337 is rather inconspicuous in statements of the form @code{if
1338 (expand(pow(a+b,2))==a*a+2*a*b+b*b) @{...@}}. Here, an implicit
1339 conversion from @code{relational} to @code{bool} takes place. Note,
1340 however, that @code{==} here does not perform any simplifications, hence
1341 @code{expand()} must be called explicitly.
1344 @node Matrices, Indexed objects, Relations, Basic Concepts
1345 @c node-name, next, previous, up
1347 @cindex @code{matrix} (class)
1349 A @dfn{matrix} is a two-dimensional array of expressions. The elements of a
1350 matrix with @math{m} rows and @math{n} columns are accessed with two
1351 @code{unsigned} indices, the first one in the range 0@dots{}@math{m-1}, the
1352 second one in the range 0@dots{}@math{n-1}.
1354 There are a couple of ways to construct matrices, with or without preset
1358 matrix::matrix(unsigned r, unsigned c);
1359 matrix::matrix(unsigned r, unsigned c, const lst & l);
1360 ex lst_to_matrix(const lst & l);
1361 ex diag_matrix(const lst & l);
1364 The first two functions are @code{matrix} constructors which create a matrix
1365 with @samp{r} rows and @samp{c} columns. The matrix elements can be
1366 initialized from a (flat) list of expressions @samp{l}. Otherwise they are
1367 all set to zero. The @code{lst_to_matrix()} function constructs a matrix
1368 from a list of lists, each list representing a matrix row. Finally,
1369 @code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
1370 elements. Note that the last two functions return expressions, not matrix
1373 Matrix elements can be accessed and set using the parenthesis (function call)
1377 const ex & matrix::operator()(unsigned r, unsigned c) const;
1378 ex & matrix::operator()(unsigned r, unsigned c);
1381 It is also possible to access the matrix elements in a linear fashion with
1382 the @code{op()} method. But C++-style subscripting with square brackets
1383 @samp{[]} is not available.
1385 Here are a couple of examples that all construct the same 2x2 diagonal
1390 symbol a("a"), b("b");
1398 e = matrix(2, 2, lst(a, 0, 0, b));
1400 e = lst_to_matrix(lst(lst(a, 0), lst(0, b)));
1402 e = diag_matrix(lst(a, b));
1409 @cindex @code{transpose()}
1410 @cindex @code{inverse()}
1411 There are three ways to do arithmetic with matrices. The first (and most
1412 efficient one) is to use the methods provided by the @code{matrix} class:
1415 matrix matrix::add(const matrix & other) const;
1416 matrix matrix::sub(const matrix & other) const;
1417 matrix matrix::mul(const matrix & other) const;
1418 matrix matrix::mul_scalar(const ex & other) const;
1419 matrix matrix::pow(const ex & expn) const;
1420 matrix matrix::transpose(void) const;
1421 matrix matrix::inverse(void) const;
1424 All of these methods return the result as a new matrix object. Here is an
1425 example that calculates @math{A*B-2*C} for three matrices @math{A}, @math{B}
1430 matrix A(2, 2, lst(1, 2, 3, 4));
1431 matrix B(2, 2, lst(-1, 0, 2, 1));
1432 matrix C(2, 2, lst(8, 4, 2, 1));
1434 matrix result = A.mul(B).sub(C.mul_scalar(2));
1435 cout << result << endl;
1436 // -> [[-13,-6],[1,2]]
1441 @cindex @code{evalm()}
1442 The second (and probably the most natural) way is to construct an expression
1443 containing matrices with the usual arithmetic operators and @code{pow()}.
1444 For efficiency reasons, expressions with sums, products and powers of
1445 matrices are not automatically evaluated in GiNaC. You have to call the
1449 ex ex::evalm() const;
1452 to obtain the result:
1459 // -> [[1,2],[3,4]]*[[-1,0],[2,1]]-2*[[8,4],[2,1]]
1460 cout << e.evalm() << endl;
1461 // -> [[-13,-6],[1,2]]
1466 The non-commutativity of the product @code{A*B} in this example is
1467 automatically recognized by GiNaC. There is no need to use a special
1468 operator here. @xref{Non-commutative objects}, for more information about
1469 dealing with non-commutative expressions.
1471 Finally, you can work with indexed matrices and call @code{simplify_indexed()}
1472 to perform the arithmetic:
1477 idx i(symbol("i"), 2), j(symbol("j"), 2), k(symbol("k"), 2);
1478 e = indexed(A, i, k) * indexed(B, k, j) - 2 * indexed(C, i, j);
1480 // -> -2*[[8,4],[2,1]].i.j+[[-1,0],[2,1]].k.j*[[1,2],[3,4]].i.k
1481 cout << e.simplify_indexed() << endl;
1482 // -> [[-13,-6],[1,2]].i.j
1486 Using indices is most useful when working with rectangular matrices and
1487 one-dimensional vectors because you don't have to worry about having to
1488 transpose matrices before multiplying them. @xref{Indexed objects}, for
1489 more information about using matrices with indices, and about indices in
1492 The @code{matrix} class provides a couple of additional methods for
1493 computing determinants, traces, and characteristic polynomials:
1496 ex matrix::determinant(unsigned algo = determinant_algo::automatic) const;
1497 ex matrix::trace(void) const;
1498 ex matrix::charpoly(const symbol & lambda) const;
1501 The @samp{algo} argument of @code{determinant()} allows to select between
1502 different algorithms for calculating the determinant. The possible values
1503 are defined in the @file{flags.h} header file. By default, GiNaC uses a
1504 heuristic to automatically select an algorithm that is likely to give the
1505 result most quickly.
1508 @node Indexed objects, Non-commutative objects, Matrices, Basic Concepts
1509 @c node-name, next, previous, up
1510 @section Indexed objects
1512 GiNaC allows you to handle expressions containing general indexed objects in
1513 arbitrary spaces. It is also able to canonicalize and simplify such
1514 expressions and perform symbolic dummy index summations. There are a number
1515 of predefined indexed objects provided, like delta and metric tensors.
1517 There are few restrictions placed on indexed objects and their indices and
1518 it is easy to construct nonsense expressions, but our intention is to
1519 provide a general framework that allows you to implement algorithms with
1520 indexed quantities, getting in the way as little as possible.
1522 @cindex @code{idx} (class)
1523 @cindex @code{indexed} (class)
1524 @subsection Indexed quantities and their indices
1526 Indexed expressions in GiNaC are constructed of two special types of objects,
1527 @dfn{index objects} and @dfn{indexed objects}.
1531 @cindex contravariant
1534 @item Index objects are of class @code{idx} or a subclass. Every index has
1535 a @dfn{value} and a @dfn{dimension} (which is the dimension of the space
1536 the index lives in) which can both be arbitrary expressions but are usually
1537 a number or a simple symbol. In addition, indices of class @code{varidx} have
1538 a @dfn{variance} (they can be co- or contravariant), and indices of class
1539 @code{spinidx} have a variance and can be @dfn{dotted} or @dfn{undotted}.
1541 @item Indexed objects are of class @code{indexed} or a subclass. They
1542 contain a @dfn{base expression} (which is the expression being indexed), and
1543 one or more indices.
1547 @strong{Note:} when printing expressions, covariant indices and indices
1548 without variance are denoted @samp{.i} while contravariant indices are
1549 denoted @samp{~i}. Dotted indices have a @samp{*} in front of the index
1550 value. In the following, we are going to use that notation in the text so
1551 instead of @math{A^i_jk} we will write @samp{A~i.j.k}. Index dimensions are
1552 not visible in the output.
1554 A simple example shall illustrate the concepts:
1558 #include <ginac/ginac.h>
1559 using namespace std;
1560 using namespace GiNaC;
1564 symbol i_sym("i"), j_sym("j");
1565 idx i(i_sym, 3), j(j_sym, 3);
1568 cout << indexed(A, i, j) << endl;
1573 The @code{idx} constructor takes two arguments, the index value and the
1574 index dimension. First we define two index objects, @code{i} and @code{j},
1575 both with the numeric dimension 3. The value of the index @code{i} is the
1576 symbol @code{i_sym} (which prints as @samp{i}) and the value of the index
1577 @code{j} is the symbol @code{j_sym} (which prints as @samp{j}). Next we
1578 construct an expression containing one indexed object, @samp{A.i.j}. It has
1579 the symbol @code{A} as its base expression and the two indices @code{i} and
1582 Note the difference between the indices @code{i} and @code{j} which are of
1583 class @code{idx}, and the index values which are the symbols @code{i_sym}
1584 and @code{j_sym}. The indices of indexed objects cannot directly be symbols
1585 or numbers but must be index objects. For example, the following is not
1586 correct and will raise an exception:
1589 symbol i("i"), j("j");
1590 e = indexed(A, i, j); // ERROR: indices must be of type idx
1593 You can have multiple indexed objects in an expression, index values can
1594 be numeric, and index dimensions symbolic:
1598 symbol B("B"), dim("dim");
1599 cout << 4 * indexed(A, i)
1600 + indexed(B, idx(j_sym, 4), idx(2, 3), idx(i_sym, dim)) << endl;
1605 @code{B} has a 4-dimensional symbolic index @samp{k}, a 3-dimensional numeric
1606 index of value 2, and a symbolic index @samp{i} with the symbolic dimension
1607 @samp{dim}. Note that GiNaC doesn't automatically notify you that the free
1608 indices of @samp{A} and @samp{B} in the sum don't match (you have to call
1609 @code{simplify_indexed()} for that, see below).
1611 In fact, base expressions, index values and index dimensions can be
1612 arbitrary expressions:
1616 cout << indexed(A+B, idx(2*i_sym+1, dim/2)) << endl;
1621 It's also possible to construct nonsense like @samp{Pi.sin(x)}. You will not
1622 get an error message from this but you will probably not be able to do
1623 anything useful with it.
1625 @cindex @code{get_value()}
1626 @cindex @code{get_dimension()}
1630 ex idx::get_value(void);
1631 ex idx::get_dimension(void);
1634 return the value and dimension of an @code{idx} object. If you have an index
1635 in an expression, such as returned by calling @code{.op()} on an indexed
1636 object, you can get a reference to the @code{idx} object with the function
1637 @code{ex_to<idx>()} on the expression.
1639 There are also the methods
1642 bool idx::is_numeric(void);
1643 bool idx::is_symbolic(void);
1644 bool idx::is_dim_numeric(void);
1645 bool idx::is_dim_symbolic(void);
1648 for checking whether the value and dimension are numeric or symbolic
1649 (non-numeric). Using the @code{info()} method of an index (see @ref{Information
1650 About Expressions}) returns information about the index value.
1652 @cindex @code{varidx} (class)
1653 If you need co- and contravariant indices, use the @code{varidx} class:
1657 symbol mu_sym("mu"), nu_sym("nu");
1658 varidx mu(mu_sym, 4), nu(nu_sym, 4); // default is contravariant ~mu, ~nu
1659 varidx mu_co(mu_sym, 4, true); // covariant index .mu
1661 cout << indexed(A, mu, nu) << endl;
1663 cout << indexed(A, mu_co, nu) << endl;
1665 cout << indexed(A, mu.toggle_variance(), nu) << endl;
1670 A @code{varidx} is an @code{idx} with an additional flag that marks it as
1671 co- or contravariant. The default is a contravariant (upper) index, but
1672 this can be overridden by supplying a third argument to the @code{varidx}
1673 constructor. The two methods
1676 bool varidx::is_covariant(void);
1677 bool varidx::is_contravariant(void);
1680 allow you to check the variance of a @code{varidx} object (use @code{ex_to<varidx>()}
1681 to get the object reference from an expression). There's also the very useful
1685 ex varidx::toggle_variance(void);
1688 which makes a new index with the same value and dimension but the opposite
1689 variance. By using it you only have to define the index once.
1691 @cindex @code{spinidx} (class)
1692 The @code{spinidx} class provides dotted and undotted variant indices, as
1693 used in the Weyl-van-der-Waerden spinor formalism:
1697 symbol K("K"), C_sym("C"), D_sym("D");
1698 spinidx C(C_sym, 2), D(D_sym); // default is 2-dimensional,
1699 // contravariant, undotted
1700 spinidx C_co(C_sym, 2, true); // covariant index
1701 spinidx D_dot(D_sym, 2, false, true); // contravariant, dotted
1702 spinidx D_co_dot(D_sym, 2, true, true); // covariant, dotted
1704 cout << indexed(K, C, D) << endl;
1706 cout << indexed(K, C_co, D_dot) << endl;
1708 cout << indexed(K, D_co_dot, D) << endl;
1713 A @code{spinidx} is a @code{varidx} with an additional flag that marks it as
1714 dotted or undotted. The default is undotted but this can be overridden by
1715 supplying a fourth argument to the @code{spinidx} constructor. The two
1719 bool spinidx::is_dotted(void);
1720 bool spinidx::is_undotted(void);
1723 allow you to check whether or not a @code{spinidx} object is dotted (use
1724 @code{ex_to<spinidx>()} to get the object reference from an expression).
1725 Finally, the two methods
1728 ex spinidx::toggle_dot(void);
1729 ex spinidx::toggle_variance_dot(void);
1732 create a new index with the same value and dimension but opposite dottedness
1733 and the same or opposite variance.
1735 @subsection Substituting indices
1737 @cindex @code{subs()}
1738 Sometimes you will want to substitute one symbolic index with another
1739 symbolic or numeric index, for example when calculating one specific element
1740 of a tensor expression. This is done with the @code{.subs()} method, as it
1741 is done for symbols (see @ref{Substituting Expressions}).
1743 You have two possibilities here. You can either substitute the whole index
1744 by another index or expression:
1748 ex e = indexed(A, mu_co);
1749 cout << e << " becomes " << e.subs(mu_co == nu) << endl;
1750 // -> A.mu becomes A~nu
1751 cout << e << " becomes " << e.subs(mu_co == varidx(0, 4)) << endl;
1752 // -> A.mu becomes A~0
1753 cout << e << " becomes " << e.subs(mu_co == 0) << endl;
1754 // -> A.mu becomes A.0
1758 The third example shows that trying to replace an index with something that
1759 is not an index will substitute the index value instead.
1761 Alternatively, you can substitute the @emph{symbol} of a symbolic index by
1766 ex e = indexed(A, mu_co);
1767 cout << e << " becomes " << e.subs(mu_sym == nu_sym) << endl;
1768 // -> A.mu becomes A.nu
1769 cout << e << " becomes " << e.subs(mu_sym == 0) << endl;
1770 // -> A.mu becomes A.0
1774 As you see, with the second method only the value of the index will get
1775 substituted. Its other properties, including its dimension, remain unchanged.
1776 If you want to change the dimension of an index you have to substitute the
1777 whole index by another one with the new dimension.
1779 Finally, substituting the base expression of an indexed object works as
1784 ex e = indexed(A, mu_co);
1785 cout << e << " becomes " << e.subs(A == A+B) << endl;
1786 // -> A.mu becomes (B+A).mu
1790 @subsection Symmetries
1791 @cindex @code{symmetry} (class)
1792 @cindex @code{sy_none()}
1793 @cindex @code{sy_symm()}
1794 @cindex @code{sy_anti()}
1795 @cindex @code{sy_cycl()}
1797 Indexed objects can have certain symmetry properties with respect to their
1798 indices. Symmetries are specified as a tree of objects of class @code{symmetry}
1799 that is constructed with the helper functions
1802 symmetry sy_none(...);
1803 symmetry sy_symm(...);
1804 symmetry sy_anti(...);
1805 symmetry sy_cycl(...);
1808 @code{sy_none()} stands for no symmetry, @code{sy_symm()} and @code{sy_anti()}
1809 specify fully symmetric or antisymmetric, respectively, and @code{sy_cycl()}
1810 represents a cyclic symmetry. Each of these functions accepts up to four
1811 arguments which can be either symmetry objects themselves or unsigned integer
1812 numbers that represent an index position (counting from 0). A symmetry
1813 specification that consists of only a single @code{sy_symm()}, @code{sy_anti()}
1814 or @code{sy_cycl()} with no arguments specifies the respective symmetry for
1817 Here are some examples of symmetry definitions:
1822 e = indexed(A, i, j);
1823 e = indexed(A, sy_none(), i, j); // equivalent
1824 e = indexed(A, sy_none(0, 1), i, j); // equivalent
1826 // Symmetric in all three indices:
1827 e = indexed(A, sy_symm(), i, j, k);
1828 e = indexed(A, sy_symm(0, 1, 2), i, j, k); // equivalent
1829 e = indexed(A, sy_symm(2, 0, 1), i, j, k); // same symmetry, but yields a
1830 // different canonical order
1832 // Symmetric in the first two indices only:
1833 e = indexed(A, sy_symm(0, 1), i, j, k);
1834 e = indexed(A, sy_none(sy_symm(0, 1), 2), i, j, k); // equivalent
1836 // Antisymmetric in the first and last index only (index ranges need not
1838 e = indexed(A, sy_anti(0, 2), i, j, k);
1839 e = indexed(A, sy_none(sy_anti(0, 2), 1), i, j, k); // equivalent
1841 // An example of a mixed symmetry: antisymmetric in the first two and
1842 // last two indices, symmetric when swapping the first and last index
1843 // pairs (like the Riemann curvature tensor):
1844 e = indexed(A, sy_symm(sy_anti(0, 1), sy_anti(2, 3)), i, j, k, l);
1846 // Cyclic symmetry in all three indices:
1847 e = indexed(A, sy_cycl(), i, j, k);
1848 e = indexed(A, sy_cycl(0, 1, 2), i, j, k); // equivalent
1850 // The following examples are invalid constructions that will throw
1851 // an exception at run time.
1853 // An index may not appear multiple times:
1854 e = indexed(A, sy_symm(0, 0, 1), i, j, k); // ERROR
1855 e = indexed(A, sy_none(sy_symm(0, 1), sy_anti(0, 2)), i, j, k); // ERROR
1857 // Every child of sy_symm(), sy_anti() and sy_cycl() must refer to the
1858 // same number of indices:
1859 e = indexed(A, sy_symm(sy_anti(0, 1), 2), i, j, k); // ERROR
1861 // And of course, you cannot specify indices which are not there:
1862 e = indexed(A, sy_symm(0, 1, 2, 3), i, j, k); // ERROR
1866 If you need to specify more than four indices, you have to use the
1867 @code{.add()} method of the @code{symmetry} class. For example, to specify
1868 full symmetry in the first six indices you would write
1869 @code{sy_symm(0, 1, 2, 3).add(4).add(5)}.
1871 If an indexed object has a symmetry, GiNaC will automatically bring the
1872 indices into a canonical order which allows for some immediate simplifications:
1876 cout << indexed(A, sy_symm(), i, j)
1877 + indexed(A, sy_symm(), j, i) << endl;
1879 cout << indexed(B, sy_anti(), i, j)
1880 + indexed(B, sy_anti(), j, i) << endl;
1882 cout << indexed(B, sy_anti(), i, j, k)
1883 - indexed(B, sy_anti(), j, k, i) << endl;
1888 @cindex @code{get_free_indices()}
1890 @subsection Dummy indices
1892 GiNaC treats certain symbolic index pairs as @dfn{dummy indices} meaning
1893 that a summation over the index range is implied. Symbolic indices which are
1894 not dummy indices are called @dfn{free indices}. Numeric indices are neither
1895 dummy nor free indices.
1897 To be recognized as a dummy index pair, the two indices must be of the same
1898 class and their value must be the same single symbol (an index like
1899 @samp{2*n+1} is never a dummy index). If the indices are of class
1900 @code{varidx} they must also be of opposite variance; if they are of class
1901 @code{spinidx} they must be both dotted or both undotted.
1903 The method @code{.get_free_indices()} returns a vector containing the free
1904 indices of an expression. It also checks that the free indices of the terms
1905 of a sum are consistent:
1909 symbol A("A"), B("B"), C("C");
1911 symbol i_sym("i"), j_sym("j"), k_sym("k"), l_sym("l");
1912 idx i(i_sym, 3), j(j_sym, 3), k(k_sym, 3), l(l_sym, 3);
1914 ex e = indexed(A, i, j) * indexed(B, j, k) + indexed(C, k, l, i, l);
1915 cout << exprseq(e.get_free_indices()) << endl;
1917 // 'j' and 'l' are dummy indices
1919 symbol mu_sym("mu"), nu_sym("nu"), rho_sym("rho"), sigma_sym("sigma");
1920 varidx mu(mu_sym, 4), nu(nu_sym, 4), rho(rho_sym, 4), sigma(sigma_sym, 4);
1922 e = indexed(A, mu, nu) * indexed(B, nu.toggle_variance(), rho)
1923 + indexed(C, mu, sigma, rho, sigma.toggle_variance());
1924 cout << exprseq(e.get_free_indices()) << endl;
1926 // 'nu' is a dummy index, but 'sigma' is not
1928 e = indexed(A, mu, mu);
1929 cout << exprseq(e.get_free_indices()) << endl;
1931 // 'mu' is not a dummy index because it appears twice with the same
1934 e = indexed(A, mu, nu) + 42;
1935 cout << exprseq(e.get_free_indices()) << endl; // ERROR
1936 // this will throw an exception:
1937 // "add::get_free_indices: inconsistent indices in sum"
1941 @cindex @code{simplify_indexed()}
1942 @subsection Simplifying indexed expressions
1944 In addition to the few automatic simplifications that GiNaC performs on
1945 indexed expressions (such as re-ordering the indices of symmetric tensors
1946 and calculating traces and convolutions of matrices and predefined tensors)
1950 ex ex::simplify_indexed(void);
1951 ex ex::simplify_indexed(const scalar_products & sp);
1954 that performs some more expensive operations:
1957 @item it checks the consistency of free indices in sums in the same way
1958 @code{get_free_indices()} does
1959 @item it tries to give dummy indices that appear in different terms of a sum
1960 the same name to allow simplifications like @math{a_i*b_i-a_j*b_j=0}
1961 @item it (symbolically) calculates all possible dummy index summations/contractions
1962 with the predefined tensors (this will be explained in more detail in the
1964 @item it detects contractions that vanish for symmetry reasons, for example
1965 the contraction of a symmetric and a totally antisymmetric tensor
1966 @item as a special case of dummy index summation, it can replace scalar products
1967 of two tensors with a user-defined value
1970 The last point is done with the help of the @code{scalar_products} class
1971 which is used to store scalar products with known values (this is not an
1972 arithmetic class, you just pass it to @code{simplify_indexed()}):
1976 symbol A("A"), B("B"), C("C"), i_sym("i");
1980 sp.add(A, B, 0); // A and B are orthogonal
1981 sp.add(A, C, 0); // A and C are orthogonal
1982 sp.add(A, A, 4); // A^2 = 4 (A has length 2)
1984 e = indexed(A + B, i) * indexed(A + C, i);
1986 // -> (B+A).i*(A+C).i
1988 cout << e.expand(expand_options::expand_indexed).simplify_indexed(sp)
1994 The @code{scalar_products} object @code{sp} acts as a storage for the
1995 scalar products added to it with the @code{.add()} method. This method
1996 takes three arguments: the two expressions of which the scalar product is
1997 taken, and the expression to replace it with. After @code{sp.add(A, B, 0)},
1998 @code{simplify_indexed()} will replace all scalar products of indexed
1999 objects that have the symbols @code{A} and @code{B} as base expressions
2000 with the single value 0. The number, type and dimension of the indices
2001 don't matter; @samp{A~mu~nu*B.mu.nu} would also be replaced by 0.
2003 @cindex @code{expand()}
2004 The example above also illustrates a feature of the @code{expand()} method:
2005 if passed the @code{expand_indexed} option it will distribute indices
2006 over sums, so @samp{(A+B).i} becomes @samp{A.i+B.i}.
2008 @cindex @code{tensor} (class)
2009 @subsection Predefined tensors
2011 Some frequently used special tensors such as the delta, epsilon and metric
2012 tensors are predefined in GiNaC. They have special properties when
2013 contracted with other tensor expressions and some of them have constant
2014 matrix representations (they will evaluate to a number when numeric
2015 indices are specified).
2017 @cindex @code{delta_tensor()}
2018 @subsubsection Delta tensor
2020 The delta tensor takes two indices, is symmetric and has the matrix
2021 representation @code{diag(1, 1, 1, ...)}. It is constructed by the function
2022 @code{delta_tensor()}:
2026 symbol A("A"), B("B");
2028 idx i(symbol("i"), 3), j(symbol("j"), 3),
2029 k(symbol("k"), 3), l(symbol("l"), 3);
2031 ex e = indexed(A, i, j) * indexed(B, k, l)
2032 * delta_tensor(i, k) * delta_tensor(j, l) << endl;
2033 cout << e.simplify_indexed() << endl;
2036 cout << delta_tensor(i, i) << endl;
2041 @cindex @code{metric_tensor()}
2042 @subsubsection General metric tensor
2044 The function @code{metric_tensor()} creates a general symmetric metric
2045 tensor with two indices that can be used to raise/lower tensor indices. The
2046 metric tensor is denoted as @samp{g} in the output and if its indices are of
2047 mixed variance it is automatically replaced by a delta tensor:
2053 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
2055 ex e = metric_tensor(mu, nu) * indexed(A, nu.toggle_variance(), rho);
2056 cout << e.simplify_indexed() << endl;
2059 e = delta_tensor(mu, nu.toggle_variance()) * metric_tensor(nu, rho);
2060 cout << e.simplify_indexed() << endl;
2063 e = metric_tensor(mu.toggle_variance(), nu.toggle_variance())
2064 * metric_tensor(nu, rho);
2065 cout << e.simplify_indexed() << endl;
2068 e = metric_tensor(nu.toggle_variance(), rho.toggle_variance())
2069 * metric_tensor(mu, nu) * (delta_tensor(mu.toggle_variance(), rho)
2070 + indexed(A, mu.toggle_variance(), rho));
2071 cout << e.simplify_indexed() << endl;
2076 @cindex @code{lorentz_g()}
2077 @subsubsection Minkowski metric tensor
2079 The Minkowski metric tensor is a special metric tensor with a constant
2080 matrix representation which is either @code{diag(1, -1, -1, ...)} (negative
2081 signature, the default) or @code{diag(-1, 1, 1, ...)} (positive signature).
2082 It is created with the function @code{lorentz_g()} (although it is output as
2087 varidx mu(symbol("mu"), 4);
2089 e = delta_tensor(varidx(0, 4), mu.toggle_variance())
2090 * lorentz_g(mu, varidx(0, 4)); // negative signature
2091 cout << e.simplify_indexed() << endl;
2094 e = delta_tensor(varidx(0, 4), mu.toggle_variance())
2095 * lorentz_g(mu, varidx(0, 4), true); // positive signature
2096 cout << e.simplify_indexed() << endl;
2101 @cindex @code{spinor_metric()}
2102 @subsubsection Spinor metric tensor
2104 The function @code{spinor_metric()} creates an antisymmetric tensor with
2105 two indices that is used to raise/lower indices of 2-component spinors.
2106 It is output as @samp{eps}:
2112 spinidx A(symbol("A")), B(symbol("B")), C(symbol("C"));
2113 ex A_co = A.toggle_variance(), B_co = B.toggle_variance();
2115 e = spinor_metric(A, B) * indexed(psi, B_co);
2116 cout << e.simplify_indexed() << endl;
2119 e = spinor_metric(A, B) * indexed(psi, A_co);
2120 cout << e.simplify_indexed() << endl;
2123 e = spinor_metric(A_co, B_co) * indexed(psi, B);
2124 cout << e.simplify_indexed() << endl;
2127 e = spinor_metric(A_co, B_co) * indexed(psi, A);
2128 cout << e.simplify_indexed() << endl;
2131 e = spinor_metric(A_co, B_co) * spinor_metric(A, B);
2132 cout << e.simplify_indexed() << endl;
2135 e = spinor_metric(A_co, B_co) * spinor_metric(B, C);
2136 cout << e.simplify_indexed() << endl;
2141 The matrix representation of the spinor metric is @code{[[0, 1], [-1, 0]]}.
2143 @cindex @code{epsilon_tensor()}
2144 @cindex @code{lorentz_eps()}
2145 @subsubsection Epsilon tensor
2147 The epsilon tensor is totally antisymmetric, its number of indices is equal
2148 to the dimension of the index space (the indices must all be of the same
2149 numeric dimension), and @samp{eps.1.2.3...} (resp. @samp{eps~0~1~2...}) is
2150 defined to be 1. Its behavior with indices that have a variance also
2151 depends on the signature of the metric. Epsilon tensors are output as
2154 There are three functions defined to create epsilon tensors in 2, 3 and 4
2158 ex epsilon_tensor(const ex & i1, const ex & i2);
2159 ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
2160 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig = false);
2163 The first two functions create an epsilon tensor in 2 or 3 Euclidean
2164 dimensions, the last function creates an epsilon tensor in a 4-dimensional
2165 Minkowski space (the last @code{bool} argument specifies whether the metric
2166 has negative or positive signature, as in the case of the Minkowski metric
2171 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4),
2172 sig(symbol("sig"), 4), lam(symbol("lam"), 4), bet(symbol("bet"), 4);
2173 e = lorentz_eps(mu, nu, rho, sig) *
2174 lorentz_eps(mu.toggle_variance(), nu.toggle_variance(), lam, bet);
2175 cout << simplify_indexed(e) << endl;
2176 // -> 2*eta~bet~rho*eta~sig~lam-2*eta~sig~bet*eta~rho~lam
2178 idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
2179 symbol A("A"), B("B");
2180 e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(B, k);
2181 cout << simplify_indexed(e) << endl;
2182 // -> -B.k*A.j*eps.i.k.j
2183 e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(A, k);
2184 cout << simplify_indexed(e) << endl;
2189 @subsection Linear algebra
2191 The @code{matrix} class can be used with indices to do some simple linear
2192 algebra (linear combinations and products of vectors and matrices, traces
2193 and scalar products):
2197 idx i(symbol("i"), 2), j(symbol("j"), 2);
2198 symbol x("x"), y("y");
2200 // A is a 2x2 matrix, X is a 2x1 vector
2201 matrix A(2, 2, lst(1, 2, 3, 4)), X(2, 1, lst(x, y));
2203 cout << indexed(A, i, i) << endl;
2206 ex e = indexed(A, i, j) * indexed(X, j);
2207 cout << e.simplify_indexed() << endl;
2208 // -> [[2*y+x],[4*y+3*x]].i
2210 e = indexed(A, i, j) * indexed(X, i) + indexed(X, j) * 2;
2211 cout << e.simplify_indexed() << endl;
2212 // -> [[3*y+3*x,6*y+2*x]].j
2216 You can of course obtain the same results with the @code{matrix::add()},
2217 @code{matrix::mul()} and @code{matrix::trace()} methods (@pxref{Matrices})
2218 but with indices you don't have to worry about transposing matrices.
2220 Matrix indices always start at 0 and their dimension must match the number
2221 of rows/columns of the matrix. Matrices with one row or one column are
2222 vectors and can have one or two indices (it doesn't matter whether it's a
2223 row or a column vector). Other matrices must have two indices.
2225 You should be careful when using indices with variance on matrices. GiNaC
2226 doesn't look at the variance and doesn't know that @samp{F~mu~nu} and
2227 @samp{F.mu.nu} are different matrices. In this case you should use only
2228 one form for @samp{F} and explicitly multiply it with a matrix representation
2229 of the metric tensor.
2232 @node Non-commutative objects, Methods and Functions, Indexed objects, Basic Concepts
2233 @c node-name, next, previous, up
2234 @section Non-commutative objects
2236 GiNaC is equipped to handle certain non-commutative algebras. Three classes of
2237 non-commutative objects are built-in which are mostly of use in high energy
2241 @item Clifford (Dirac) algebra (class @code{clifford})
2242 @item su(3) Lie algebra (class @code{color})
2243 @item Matrices (unindexed) (class @code{matrix})
2246 The @code{clifford} and @code{color} classes are subclasses of
2247 @code{indexed} because the elements of these algebras usually carry
2248 indices. The @code{matrix} class is described in more detail in
2251 Unlike most computer algebra systems, GiNaC does not primarily provide an
2252 operator (often denoted @samp{&*}) for representing inert products of
2253 arbitrary objects. Rather, non-commutativity in GiNaC is a property of the
2254 classes of objects involved, and non-commutative products are formed with
2255 the usual @samp{*} operator, as are ordinary products. GiNaC is capable of
2256 figuring out by itself which objects commute and will group the factors
2257 by their class. Consider this example:
2261 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
2262 idx a(symbol("a"), 8), b(symbol("b"), 8);
2263 ex e = -dirac_gamma(mu) * (2*color_T(a)) * 8 * color_T(b) * dirac_gamma(nu);
2265 // -> -16*(gamma~mu*gamma~nu)*(T.a*T.b)
2269 As can be seen, GiNaC pulls out the overall commutative factor @samp{-16} and
2270 groups the non-commutative factors (the gammas and the su(3) generators)
2271 together while preserving the order of factors within each class (because
2272 Clifford objects commute with color objects). The resulting expression is a
2273 @emph{commutative} product with two factors that are themselves non-commutative
2274 products (@samp{gamma~mu*gamma~nu} and @samp{T.a*T.b}). For clarification,
2275 parentheses are placed around the non-commutative products in the output.
2277 @cindex @code{ncmul} (class)
2278 Non-commutative products are internally represented by objects of the class
2279 @code{ncmul}, as opposed to commutative products which are handled by the
2280 @code{mul} class. You will normally not have to worry about this distinction,
2283 The advantage of this approach is that you never have to worry about using
2284 (or forgetting to use) a special operator when constructing non-commutative
2285 expressions. Also, non-commutative products in GiNaC are more intelligent
2286 than in other computer algebra systems; they can, for example, automatically
2287 canonicalize themselves according to rules specified in the implementation
2288 of the non-commutative classes. The drawback is that to work with other than
2289 the built-in algebras you have to implement new classes yourself. Symbols
2290 always commute and it's not possible to construct non-commutative products
2291 using symbols to represent the algebra elements or generators. User-defined
2292 functions can, however, be specified as being non-commutative.
2294 @cindex @code{return_type()}
2295 @cindex @code{return_type_tinfo()}
2296 Information about the commutativity of an object or expression can be
2297 obtained with the two member functions
2300 unsigned ex::return_type(void) const;
2301 unsigned ex::return_type_tinfo(void) const;
2304 The @code{return_type()} function returns one of three values (defined in
2305 the header file @file{flags.h}), corresponding to three categories of
2306 expressions in GiNaC:
2309 @item @code{return_types::commutative}: Commutes with everything. Most GiNaC
2310 classes are of this kind.
2311 @item @code{return_types::noncommutative}: Non-commutative, belonging to a
2312 certain class of non-commutative objects which can be determined with the
2313 @code{return_type_tinfo()} method. Expressions of this category commute
2314 with everything except @code{noncommutative} expressions of the same
2316 @item @code{return_types::noncommutative_composite}: Non-commutative, composed
2317 of non-commutative objects of different classes. Expressions of this
2318 category don't commute with any other @code{noncommutative} or
2319 @code{noncommutative_composite} expressions.
2322 The value returned by the @code{return_type_tinfo()} method is valid only
2323 when the return type of the expression is @code{noncommutative}. It is a
2324 value that is unique to the class of the object and usually one of the
2325 constants in @file{tinfos.h}, or derived therefrom.
2327 Here are a couple of examples:
2330 @multitable @columnfractions 0.33 0.33 0.34
2331 @item @strong{Expression} @tab @strong{@code{return_type()}} @tab @strong{@code{return_type_tinfo()}}
2332 @item @code{42} @tab @code{commutative} @tab -
2333 @item @code{2*x-y} @tab @code{commutative} @tab -
2334 @item @code{dirac_ONE()} @tab @code{noncommutative} @tab @code{TINFO_clifford}
2335 @item @code{dirac_gamma(mu)*dirac_gamma(nu)} @tab @code{noncommutative} @tab @code{TINFO_clifford}
2336 @item @code{2*color_T(a)} @tab @code{noncommutative} @tab @code{TINFO_color}
2337 @item @code{dirac_ONE()*color_T(a)} @tab @code{noncommutative_composite} @tab -
2341 Note: the @code{return_type_tinfo()} of Clifford objects is only equal to
2342 @code{TINFO_clifford} for objects with a representation label of zero.
2343 Other representation labels yield a different @code{return_type_tinfo()},
2344 but it's the same for any two objects with the same label. This is also true
2347 A last note: With the exception of matrices, positive integer powers of
2348 non-commutative objects are automatically expanded in GiNaC. For example,
2349 @code{pow(a*b, 2)} becomes @samp{a*b*a*b} if @samp{a} and @samp{b} are
2350 non-commutative expressions).
2353 @cindex @code{clifford} (class)
2354 @subsection Clifford algebra
2356 @cindex @code{dirac_gamma()}
2357 Clifford algebra elements (also called Dirac gamma matrices, although GiNaC
2358 doesn't treat them as matrices) are designated as @samp{gamma~mu} and satisfy
2359 @samp{gamma~mu*gamma~nu + gamma~nu*gamma~mu = 2*eta~mu~nu} where @samp{eta~mu~nu}
2360 is the Minkowski metric tensor. Dirac gammas are constructed by the function
2363 ex dirac_gamma(const ex & mu, unsigned char rl = 0);
2366 which takes two arguments: the index and a @dfn{representation label} in the
2367 range 0 to 255 which is used to distinguish elements of different Clifford
2368 algebras (this is also called a @dfn{spin line index}). Gammas with different
2369 labels commute with each other. The dimension of the index can be 4 or (in
2370 the framework of dimensional regularization) any symbolic value. Spinor
2371 indices on Dirac gammas are not supported in GiNaC.
2373 @cindex @code{dirac_ONE()}
2374 The unity element of a Clifford algebra is constructed by
2377 ex dirac_ONE(unsigned char rl = 0);
2380 @strong{Note:} You must always use @code{dirac_ONE()} when referring to
2381 multiples of the unity element, even though it's customary to omit it.
2382 E.g. instead of @code{dirac_gamma(mu)*(dirac_slash(q,4)+m)} you have to
2383 write @code{dirac_gamma(mu)*(dirac_slash(q,4)+m*dirac_ONE())}. Otherwise,
2384 GiNaC will complain and/or produce incorrect results.
2386 @cindex @code{dirac_gamma5()}
2387 There is a special element @samp{gamma5} that commutes with all other
2388 gammas, has a unit square, and in 4 dimensions equals
2389 @samp{gamma~0 gamma~1 gamma~2 gamma~3}, provided by
2392 ex dirac_gamma5(unsigned char rl = 0);
2395 @cindex @code{dirac_gammaL()}
2396 @cindex @code{dirac_gammaR()}
2397 The chiral projectors @samp{(1+/-gamma5)/2} are also available as proper
2398 objects, constructed by
2401 ex dirac_gammaL(unsigned char rl = 0);
2402 ex dirac_gammaR(unsigned char rl = 0);
2405 They observe the relations @samp{gammaL^2 = gammaL}, @samp{gammaR^2 = gammaR},
2406 and @samp{gammaL gammaR = gammaR gammaL = 0}.
2408 @cindex @code{dirac_slash()}
2409 Finally, the function
2412 ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
2415 creates a term that represents a contraction of @samp{e} with the Dirac
2416 Lorentz vector (it behaves like a term of the form @samp{e.mu gamma~mu}
2417 with a unique index whose dimension is given by the @code{dim} argument).
2418 Such slashed expressions are printed with a trailing backslash, e.g. @samp{e\}.
2420 In products of dirac gammas, superfluous unity elements are automatically
2421 removed, squares are replaced by their values, and @samp{gamma5}, @samp{gammaL}
2422 and @samp{gammaR} are moved to the front.
2424 The @code{simplify_indexed()} function performs contractions in gamma strings,
2430 symbol a("a"), b("b"), D("D");
2431 varidx mu(symbol("mu"), D);
2432 ex e = dirac_gamma(mu) * dirac_slash(a, D)
2433 * dirac_gamma(mu.toggle_variance());
2435 // -> gamma~mu*a\*gamma.mu
2436 e = e.simplify_indexed();
2439 cout << e.subs(D == 4) << endl;
2445 @cindex @code{dirac_trace()}
2446 To calculate the trace of an expression containing strings of Dirac gammas
2447 you use the function
2450 ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4);
2453 This function takes the trace of all gammas with the specified representation
2454 label; gammas with other labels are left standing. The last argument to
2455 @code{dirac_trace()} is the value to be returned for the trace of the unity
2456 element, which defaults to 4. The @code{dirac_trace()} function is a linear
2457 functional that is equal to the usual trace only in @math{D = 4} dimensions.
2458 In particular, the functional is not cyclic in @math{D != 4} dimensions when
2459 acting on expressions containing @samp{gamma5}, so it's not a proper trace.
2460 This @samp{gamma5} scheme is described in greater detail in
2461 @cite{The Role of gamma5 in Dimensional Regularization}.
2463 The value of the trace itself is also usually different in 4 and in
2464 @math{D != 4} dimensions:
2469 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
2470 ex e = dirac_gamma(mu) * dirac_gamma(nu) *
2471 dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
2472 cout << dirac_trace(e).simplify_indexed() << endl;
2479 varidx mu(symbol("mu"), D), nu(symbol("nu"), D), rho(symbol("rho"), D);
2480 ex e = dirac_gamma(mu) * dirac_gamma(nu) *
2481 dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
2482 cout << dirac_trace(e).simplify_indexed() << endl;
2483 // -> 8*eta~rho~nu-4*eta~rho~nu*D
2487 Here is an example for using @code{dirac_trace()} to compute a value that
2488 appears in the calculation of the one-loop vacuum polarization amplitude in
2493 symbol q("q"), l("l"), m("m"), ldotq("ldotq"), D("D");
2494 varidx mu(symbol("mu"), D), nu(symbol("nu"), D);
2497 sp.add(l, l, pow(l, 2));
2498 sp.add(l, q, ldotq);
2500 ex e = dirac_gamma(mu) *
2501 (dirac_slash(l, D) + dirac_slash(q, D) + m * dirac_ONE()) *
2502 dirac_gamma(mu.toggle_variance()) *
2503 (dirac_slash(l, D) + m * dirac_ONE());
2504 e = dirac_trace(e).simplify_indexed(sp);
2505 e = e.collect(lst(l, ldotq, m));
2507 // -> (8-4*D)*l^2+(8-4*D)*ldotq+4*D*m^2
2511 The @code{canonicalize_clifford()} function reorders all gamma products that
2512 appear in an expression to a canonical (but not necessarily simple) form.
2513 You can use this to compare two expressions or for further simplifications:
2517 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
2518 ex e = dirac_gamma(mu) * dirac_gamma(nu) + dirac_gamma(nu) * dirac_gamma(mu);
2520 // -> gamma~mu*gamma~nu+gamma~nu*gamma~mu
2522 e = canonicalize_clifford(e);
2529 @cindex @code{color} (class)
2530 @subsection Color algebra
2532 @cindex @code{color_T()}
2533 For computations in quantum chromodynamics, GiNaC implements the base elements
2534 and structure constants of the su(3) Lie algebra (color algebra). The base
2535 elements @math{T_a} are constructed by the function
2538 ex color_T(const ex & a, unsigned char rl = 0);
2541 which takes two arguments: the index and a @dfn{representation label} in the
2542 range 0 to 255 which is used to distinguish elements of different color
2543 algebras. Objects with different labels commute with each other. The
2544 dimension of the index must be exactly 8 and it should be of class @code{idx},
2547 @cindex @code{color_ONE()}
2548 The unity element of a color algebra is constructed by
2551 ex color_ONE(unsigned char rl = 0);
2554 @strong{Note:} You must always use @code{color_ONE()} when referring to
2555 multiples of the unity element, even though it's customary to omit it.
2556 E.g. instead of @code{color_T(a)*(color_T(b)*indexed(X,b)+1)} you have to
2557 write @code{color_T(a)*(color_T(b)*indexed(X,b)+color_ONE())}. Otherwise,
2558 GiNaC may produce incorrect results.
2560 @cindex @code{color_d()}
2561 @cindex @code{color_f()}
2565 ex color_d(const ex & a, const ex & b, const ex & c);
2566 ex color_f(const ex & a, const ex & b, const ex & c);
2569 create the symmetric and antisymmetric structure constants @math{d_abc} and
2570 @math{f_abc} which satisfy @math{@{T_a, T_b@} = 1/3 delta_ab + d_abc T_c}
2571 and @math{[T_a, T_b] = i f_abc T_c}.
2573 @cindex @code{color_h()}
2574 There's an additional function
2577 ex color_h(const ex & a, const ex & b, const ex & c);
2580 which returns the linear combination @samp{color_d(a, b, c)+I*color_f(a, b, c)}.
2582 The function @code{simplify_indexed()} performs some simplifications on
2583 expressions containing color objects:
2588 idx a(symbol("a"), 8), b(symbol("b"), 8), c(symbol("c"), 8),
2589 k(symbol("k"), 8), l(symbol("l"), 8);
2591 e = color_d(a, b, l) * color_f(a, b, k);
2592 cout << e.simplify_indexed() << endl;
2595 e = color_d(a, b, l) * color_d(a, b, k);
2596 cout << e.simplify_indexed() << endl;
2599 e = color_f(l, a, b) * color_f(a, b, k);
2600 cout << e.simplify_indexed() << endl;
2603 e = color_h(a, b, c) * color_h(a, b, c);
2604 cout << e.simplify_indexed() << endl;
2607 e = color_h(a, b, c) * color_T(b) * color_T(c);
2608 cout << e.simplify_indexed() << endl;
2611 e = color_h(a, b, c) * color_T(a) * color_T(b) * color_T(c);
2612 cout << e.simplify_indexed() << endl;
2615 e = color_T(k) * color_T(a) * color_T(b) * color_T(k);
2616 cout << e.simplify_indexed() << endl;
2617 // -> 1/4*delta.b.a*ONE-1/6*T.a*T.b
2621 @cindex @code{color_trace()}
2622 To calculate the trace of an expression containing color objects you use the
2626 ex color_trace(const ex & e, unsigned char rl = 0);
2629 This function takes the trace of all color @samp{T} objects with the
2630 specified representation label; @samp{T}s with other labels are left
2631 standing. For example:
2635 e = color_trace(4 * color_T(a) * color_T(b) * color_T(c));
2637 // -> -I*f.a.c.b+d.a.c.b
2642 @node Methods and Functions, Information About Expressions, Non-commutative objects, Top
2643 @c node-name, next, previous, up
2644 @chapter Methods and Functions
2647 In this chapter the most important algorithms provided by GiNaC will be
2648 described. Some of them are implemented as functions on expressions,
2649 others are implemented as methods provided by expression objects. If
2650 they are methods, there exists a wrapper function around it, so you can
2651 alternatively call it in a functional way as shown in the simple
2656 cout << "As method: " << sin(1).evalf() << endl;
2657 cout << "As function: " << evalf(sin(1)) << endl;
2661 @cindex @code{subs()}
2662 The general rule is that wherever methods accept one or more parameters
2663 (@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
2664 wrapper accepts is the same but preceded by the object to act on
2665 (@var{object}, @var{arg1}, @var{arg2}, @dots{}). This approach is the
2666 most natural one in an OO model but it may lead to confusion for MapleV
2667 users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
2668 would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
2669 @code{A} and @code{x}). On the other hand, since MapleV returns 3 on
2670 @code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
2671 coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
2672 here. Also, users of MuPAD will in most cases feel more comfortable
2673 with GiNaC's convention. All function wrappers are implemented
2674 as simple inline functions which just call the corresponding method and
2675 are only provided for users uncomfortable with OO who are dead set to
2676 avoid method invocations. Generally, nested function wrappers are much
2677 harder to read than a sequence of methods and should therefore be
2678 avoided if possible. On the other hand, not everything in GiNaC is a
2679 method on class @code{ex} and sometimes calling a function cannot be
2683 * Information About Expressions::
2684 * Substituting Expressions::
2685 * Pattern Matching and Advanced Substitutions::
2686 * Applying a Function on Subexpressions::
2687 * Polynomial Arithmetic:: Working with polynomials.
2688 * Rational Expressions:: Working with rational functions.
2689 * Symbolic Differentiation::
2690 * Series Expansion:: Taylor and Laurent expansion.
2692 * Built-in Functions:: List of predefined mathematical functions.
2693 * Input/Output:: Input and output of expressions.
2697 @node Information About Expressions, Substituting Expressions, Methods and Functions, Methods and Functions
2698 @c node-name, next, previous, up
2699 @section Getting information about expressions
2701 @subsection Checking expression types
2702 @cindex @code{is_a<@dots{}>()}
2703 @cindex @code{is_exactly_a<@dots{}>()}
2704 @cindex @code{ex_to<@dots{}>()}
2705 @cindex Converting @code{ex} to other classes
2706 @cindex @code{info()}
2707 @cindex @code{return_type()}
2708 @cindex @code{return_type_tinfo()}
2710 Sometimes it's useful to check whether a given expression is a plain number,
2711 a sum, a polynomial with integer coefficients, or of some other specific type.
2712 GiNaC provides a couple of functions for this:
2715 bool is_a<T>(const ex & e);
2716 bool is_exactly_a<T>(const ex & e);
2717 bool ex::info(unsigned flag);
2718 unsigned ex::return_type(void) const;
2719 unsigned ex::return_type_tinfo(void) const;
2722 When the test made by @code{is_a<T>()} returns true, it is safe to call
2723 one of the functions @code{ex_to<T>()}, where @code{T} is one of the
2724 class names (@xref{The Class Hierarchy}, for a list of all classes). For
2725 example, assuming @code{e} is an @code{ex}:
2730 if (is_a<numeric>(e))
2731 numeric n = ex_to<numeric>(e);
2736 @code{is_a<T>(e)} allows you to check whether the top-level object of
2737 an expression @samp{e} is an instance of the GiNaC class @samp{T}
2738 (@xref{The Class Hierarchy}, for a list of all classes). This is most useful,
2739 e.g., for checking whether an expression is a number, a sum, or a product:
2746 is_a<numeric>(e1); // true
2747 is_a<numeric>(e2); // false
2748 is_a<add>(e1); // false
2749 is_a<add>(e2); // true
2750 is_a<mul>(e1); // false
2751 is_a<mul>(e2); // false
2755 In contrast, @code{is_exactly_a<T>(e)} allows you to check whether the
2756 top-level object of an expression @samp{e} is an instance of the GiNaC
2757 class @samp{T}, not including parent classes.
2759 The @code{info()} method is used for checking certain attributes of
2760 expressions. The possible values for the @code{flag} argument are defined
2761 in @file{ginac/flags.h}, the most important being explained in the following
2765 @multitable @columnfractions .30 .70
2766 @item @strong{Flag} @tab @strong{Returns true if the object is@dots{}}
2767 @item @code{numeric}
2768 @tab @dots{}a number (same as @code{is_<numeric>(...)})
2770 @tab @dots{}a real integer, rational or float (i.e. is not complex)
2771 @item @code{rational}
2772 @tab @dots{}an exact rational number (integers are rational, too)
2773 @item @code{integer}
2774 @tab @dots{}a (non-complex) integer
2775 @item @code{crational}
2776 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
2777 @item @code{cinteger}
2778 @tab @dots{}a (complex) integer (such as @math{2-3*I})
2779 @item @code{positive}
2780 @tab @dots{}not complex and greater than 0
2781 @item @code{negative}
2782 @tab @dots{}not complex and less than 0
2783 @item @code{nonnegative}
2784 @tab @dots{}not complex and greater than or equal to 0
2786 @tab @dots{}an integer greater than 0
2788 @tab @dots{}an integer less than 0
2789 @item @code{nonnegint}
2790 @tab @dots{}an integer greater than or equal to 0
2792 @tab @dots{}an even integer
2794 @tab @dots{}an odd integer
2796 @tab @dots{}a prime integer (probabilistic primality test)
2797 @item @code{relation}
2798 @tab @dots{}a relation (same as @code{is_a<relational>(...)})
2799 @item @code{relation_equal}
2800 @tab @dots{}a @code{==} relation
2801 @item @code{relation_not_equal}
2802 @tab @dots{}a @code{!=} relation
2803 @item @code{relation_less}
2804 @tab @dots{}a @code{<} relation
2805 @item @code{relation_less_or_equal}
2806 @tab @dots{}a @code{<=} relation
2807 @item @code{relation_greater}
2808 @tab @dots{}a @code{>} relation
2809 @item @code{relation_greater_or_equal}
2810 @tab @dots{}a @code{>=} relation
2812 @tab @dots{}a symbol (same as @code{is_a<symbol>(...)})
2814 @tab @dots{}a list (same as @code{is_a<lst>(...)})
2815 @item @code{polynomial}
2816 @tab @dots{}a polynomial (i.e. only consists of sums and products of numbers and symbols with positive integer powers)
2817 @item @code{integer_polynomial}
2818 @tab @dots{}a polynomial with (non-complex) integer coefficients
2819 @item @code{cinteger_polynomial}
2820 @tab @dots{}a polynomial with (possibly complex) integer coefficients (such as @math{2-3*I})
2821 @item @code{rational_polynomial}
2822 @tab @dots{}a polynomial with (non-complex) rational coefficients
2823 @item @code{crational_polynomial}
2824 @tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I})
2825 @item @code{rational_function}
2826 @tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)})
2827 @item @code{algebraic}
2828 @tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1})
2832 To determine whether an expression is commutative or non-commutative and if
2833 so, with which other expressions it would commute, you use the methods
2834 @code{return_type()} and @code{return_type_tinfo()}. @xref{Non-commutative objects},
2835 for an explanation of these.
2838 @subsection Accessing subexpressions
2839 @cindex @code{nops()}
2842 @cindex @code{relational} (class)
2844 GiNaC provides the two methods
2847 unsigned ex::nops();
2848 ex ex::op(unsigned i);
2851 for accessing the subexpressions in the container-like GiNaC classes like
2852 @code{add}, @code{mul}, @code{lst}, and @code{function}. @code{nops()}
2853 determines the number of subexpressions (@samp{operands}) contained, while
2854 @code{op()} returns the @code{i}-th (0..@code{nops()-1}) subexpression.
2855 In the case of a @code{power} object, @code{op(0)} will return the basis
2856 and @code{op(1)} the exponent. For @code{indexed} objects, @code{op(0)}
2857 is the base expression and @code{op(i)}, @math{i>0} are the indices.
2859 The left-hand and right-hand side expressions of objects of class
2860 @code{relational} (and only of these) can also be accessed with the methods
2868 @subsection Comparing expressions
2869 @cindex @code{is_equal()}
2870 @cindex @code{is_zero()}
2872 Expressions can be compared with the usual C++ relational operators like
2873 @code{==}, @code{>}, and @code{<} but if the expressions contain symbols,
2874 the result is usually not determinable and the result will be @code{false},
2875 except in the case of the @code{!=} operator. You should also be aware that
2876 GiNaC will only do the most trivial test for equality (subtracting both
2877 expressions), so something like @code{(pow(x,2)+x)/x==x+1} will return
2880 Actually, if you construct an expression like @code{a == b}, this will be
2881 represented by an object of the @code{relational} class (@pxref{Relations})
2882 which is not evaluated until (explicitly or implicitly) cast to a @code{bool}.
2884 There are also two methods
2887 bool ex::is_equal(const ex & other);
2891 for checking whether one expression is equal to another, or equal to zero,
2894 @strong{Warning:} You will also find an @code{ex::compare()} method in the
2895 GiNaC header files. This method is however only to be used internally by
2896 GiNaC to establish a canonical sort order for terms, and using it to compare
2897 expressions will give very surprising results.
2900 @node Substituting Expressions, Pattern Matching and Advanced Substitutions, Information About Expressions, Methods and Functions
2901 @c node-name, next, previous, up
2902 @section Substituting expressions
2903 @cindex @code{subs()}
2905 Algebraic objects inside expressions can be replaced with arbitrary
2906 expressions via the @code{.subs()} method:
2909 ex ex::subs(const ex & e);
2910 ex ex::subs(const lst & syms, const lst & repls);
2913 In the first form, @code{subs()} accepts a relational of the form
2914 @samp{object == expression} or a @code{lst} of such relationals:
2918 symbol x("x"), y("y");
2920 ex e1 = 2*x^2-4*x+3;
2921 cout << "e1(7) = " << e1.subs(x == 7) << endl;
2925 cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
2930 If you specify multiple substitutions, they are performed in parallel, so e.g.
2931 @code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
2933 The second form of @code{subs()} takes two lists, one for the objects to be
2934 replaced and one for the expressions to be substituted (both lists must
2935 contain the same number of elements). Using this form, you would write
2936 @code{subs(lst(x, y), lst(y, x))} to exchange @samp{x} and @samp{y}.
2938 @code{subs()} performs syntactic substitution of any complete algebraic
2939 object; it does not try to match sub-expressions as is demonstrated by the
2944 symbol x("x"), y("y"), z("z");
2946 ex e1 = pow(x+y, 2);
2947 cout << e1.subs(x+y == 4) << endl;
2950 ex e2 = sin(x)*sin(y)*cos(x);
2951 cout << e2.subs(sin(x) == cos(x)) << endl;
2952 // -> cos(x)^2*sin(y)
2955 cout << e3.subs(x+y == 4) << endl;
2957 // (and not 4+z as one might expect)
2961 A more powerful form of substitution using wildcards is described in the
2965 @node Pattern Matching and Advanced Substitutions, Applying a Function on Subexpressions, Substituting Expressions, Methods and Functions
2966 @c node-name, next, previous, up
2967 @section Pattern matching and advanced substitutions
2968 @cindex @code{wildcard} (class)
2969 @cindex Pattern matching
2971 GiNaC allows the use of patterns for checking whether an expression is of a
2972 certain form or contains subexpressions of a certain form, and for
2973 substituting expressions in a more general way.
2975 A @dfn{pattern} is an algebraic expression that optionally contains wildcards.
2976 A @dfn{wildcard} is a special kind of object (of class @code{wildcard}) that
2977 represents an arbitrary expression. Every wildcard has a @dfn{label} which is
2978 an unsigned integer number to allow having multiple different wildcards in a
2979 pattern. Wildcards are printed as @samp{$label} (this is also the way they
2980 are specified in @command{ginsh}). In C++ code, wildcard objects are created
2984 ex wild(unsigned label = 0);
2987 which is simply a wrapper for the @code{wildcard()} constructor with a shorter
2990 Some examples for patterns:
2992 @multitable @columnfractions .5 .5
2993 @item @strong{Constructed as} @tab @strong{Output as}
2994 @item @code{wild()} @tab @samp{$0}
2995 @item @code{pow(x,wild())} @tab @samp{x^$0}
2996 @item @code{atan2(wild(1),wild(2))} @tab @samp{atan2($1,$2)}
2997 @item @code{indexed(A,idx(wild(),3))} @tab @samp{A.$0}
3003 @item Wildcards behave like symbols and are subject to the same algebraic
3004 rules. E.g., @samp{$0+2*$0} is automatically transformed to @samp{3*$0}.
3005 @item As shown in the last example, to use wildcards for indices you have to
3006 use them as the value of an @code{idx} object. This is because indices must
3007 always be of class @code{idx} (or a subclass).
3008 @item Wildcards only represent expressions or subexpressions. It is not
3009 possible to use them as placeholders for other properties like index
3010 dimension or variance, representation labels, symmetry of indexed objects
3012 @item Because wildcards are commutative, it is not possible to use wildcards
3013 as part of noncommutative products.
3014 @item A pattern does not have to contain wildcards. @samp{x} and @samp{x+y}
3015 are also valid patterns.
3018 @cindex @code{match()}
3019 The most basic application of patterns is to check whether an expression
3020 matches a given pattern. This is done by the function
3023 bool ex::match(const ex & pattern);
3024 bool ex::match(const ex & pattern, lst & repls);
3027 This function returns @code{true} when the expression matches the pattern
3028 and @code{false} if it doesn't. If used in the second form, the actual
3029 subexpressions matched by the wildcards get returned in the @code{repls}
3030 object as a list of relations of the form @samp{wildcard == expression}.
3031 If @code{match()} returns false, the state of @code{repls} is undefined.
3032 For reproducible results, the list should be empty when passed to
3033 @code{match()}, but it is also possible to find similarities in multiple
3034 expressions by passing in the result of a previous match.
3036 The matching algorithm works as follows:
3039 @item A single wildcard matches any expression. If one wildcard appears
3040 multiple times in a pattern, it must match the same expression in all
3041 places (e.g. @samp{$0} matches anything, and @samp{$0*($0+1)} matches
3042 @samp{x*(x+1)} but not @samp{x*(y+1)}).
3043 @item If the expression is not of the same class as the pattern, the match
3044 fails (i.e. a sum only matches a sum, a function only matches a function,
3046 @item If the pattern is a function, it only matches the same function
3047 (i.e. @samp{sin($0)} matches @samp{sin(x)} but doesn't match @samp{exp(x)}).
3048 @item Except for sums and products, the match fails if the number of
3049 subexpressions (@code{nops()}) is not equal to the number of subexpressions
3051 @item If there are no subexpressions, the expressions and the pattern must
3052 be equal (in the sense of @code{is_equal()}).
3053 @item Except for sums and products, each subexpression (@code{op()}) must
3054 match the corresponding subexpression of the pattern.
3057 Sums (@code{add}) and products (@code{mul}) are treated in a special way to
3058 account for their commutativity and associativity:
3061 @item If the pattern contains a term or factor that is a single wildcard,
3062 this one is used as the @dfn{global wildcard}. If there is more than one
3063 such wildcard, one of them is chosen as the global wildcard in a random
3065 @item Every term/factor of the pattern, except the global wildcard, is
3066 matched against every term of the expression in sequence. If no match is
3067 found, the whole match fails. Terms that did match are not considered in
3069 @item If there are no unmatched terms left, the match succeeds. Otherwise
3070 the match fails unless there is a global wildcard in the pattern, in
3071 which case this wildcard matches the remaining terms.
3074 In general, having more than one single wildcard as a term of a sum or a
3075 factor of a product (such as @samp{a+$0+$1}) will lead to unpredictable or
3078 Here are some examples in @command{ginsh} to demonstrate how it works (the
3079 @code{match()} function in @command{ginsh} returns @samp{FAIL} if the
3080 match fails, and the list of wildcard replacements otherwise):
3083 > match((x+y)^a,(x+y)^a);
3085 > match((x+y)^a,(x+y)^b);
3087 > match((x+y)^a,$1^$2);
3089 > match((x+y)^a,$1^$1);
3091 > match((x+y)^(x+y),$1^$1);
3093 > match((x+y)^(x+y),$1^$2);
3095 > match((a+b)*(a+c),($1+b)*($1+c));
3097 > match((a+b)*(a+c),(a+$1)*(a+$2));
3099 (Unpredictable. The result might also be [$1==c,$2==b].)
3100 > match((a+b)*(a+c),($1+$2)*($1+$3));
3101 (The result is undefined. Due to the sequential nature of the algorithm
3102 and the re-ordering of terms in GiNaC, the match for the first factor
3103 may be @{$1==a,$2==b@} in which case the match for the second factor
3104 succeeds, or it may be @{$1==b,$2==a@} which causes the second match to
3106 > match(a*(x+y)+a*z+b,a*$1+$2);
3107 (This is also ambiguous and may return either @{$1==z,$2==a*(x+y)+b@} or
3108 @{$1=x+y,$2=a*z+b@}.)
3109 > match(a+b+c+d+e+f,c);
3111 > match(a+b+c+d+e+f,c+$0);
3113 > match(a+b+c+d+e+f,c+e+$0);
3115 > match(a+b,a+b+$0);
3117 > match(a*b^2,a^$1*b^$2);
3119 (The matching is syntactic, not algebraic, and "a" doesn't match "a^$1"
3120 even though a==a^1.)
3121 > match(x*atan2(x,x^2),$0*atan2($0,$0^2));
3123 > match(atan2(y,x^2),atan2(y,$0));
3127 @cindex @code{has()}
3128 A more general way to look for patterns in expressions is provided by the
3132 bool ex::has(const ex & pattern);
3135 This function checks whether a pattern is matched by an expression itself or
3136 by any of its subexpressions.
3138 Again some examples in @command{ginsh} for illustration (in @command{ginsh},
3139 @code{has()} returns @samp{1} for @code{true} and @samp{0} for @code{false}):
3142 > has(x*sin(x+y+2*a),y);
3144 > has(x*sin(x+y+2*a),x+y);
3146 (This is because in GiNaC, "x+y" is not a subexpression of "x+y+2*a" (which
3147 has the subexpressions "x", "y" and "2*a".)
3148 > has(x*sin(x+y+2*a),x+y+$1);
3150 (But this is possible.)
3151 > has(x*sin(2*(x+y)+2*a),x+y);
3153 (This fails because "2*(x+y)" automatically gets converted to "2*x+2*y" of
3154 which "x+y" is not a subexpression.)
3157 (Although x^1==x and x^0==1, neither "x" nor "1" are actually of the form
3159 > has(4*x^2-x+3,$1*x);
3161 > has(4*x^2+x+3,$1*x);
3163 (Another possible pitfall. The first expression matches because the term
3164 "-x" has the form "(-1)*x" in GiNaC. To check whether a polynomial
3165 contains a linear term you should use the coeff() function instead.)
3168 @cindex @code{find()}
3172 bool ex::find(const ex & pattern, lst & found);
3175 works a bit like @code{has()} but it doesn't stop upon finding the first
3176 match. Instead, it appends all found matches to the specified list. If there
3177 are multiple occurrences of the same expression, it is entered only once to
3178 the list. @code{find()} returns false if no matches were found (in
3179 @command{ginsh}, it returns an empty list):
3182 > find(1+x+x^2+x^3,x);
3184 > find(1+x+x^2+x^3,y);
3186 > find(1+x+x^2+x^3,x^$1);
3188 (Note the absence of "x".)
3189 > expand((sin(x)+sin(y))*(a+b));
3190 sin(y)*a+sin(x)*b+sin(x)*a+sin(y)*b
3195 @cindex @code{subs()}
3196 Probably the most useful application of patterns is to use them for
3197 substituting expressions with the @code{subs()} method. Wildcards can be
3198 used in the search patterns as well as in the replacement expressions, where
3199 they get replaced by the expressions matched by them. @code{subs()} doesn't
3200 know anything about algebra; it performs purely syntactic substitutions.
3205 > subs(a^2+b^2+(x+y)^2,$1^2==$1^3);
3207 > subs(a^4+b^4+(x+y)^4,$1^2==$1^3);
3209 > subs((a+b+c)^2,a+b==x);
3211 > subs((a+b+c)^2,a+b+$1==x+$1);
3213 > subs(a+2*b,a+b==x);
3215 > subs(4*x^3-2*x^2+5*x-1,x==a);
3217 > subs(4*x^3-2*x^2+5*x-1,x^$0==a^$0);
3219 > subs(sin(1+sin(x)),sin($1)==cos($1));
3221 > expand(subs(a*sin(x+y)^2+a*cos(x+y)^2+b,cos($1)^2==1-sin($1)^2));
3225 The last example would be written in C++ in this way:
3229 symbol a("a"), b("b"), x("x"), y("y");
3230 e = a*pow(sin(x+y), 2) + a*pow(cos(x+y), 2) + b;
3231 e = e.subs(pow(cos(wild()), 2) == 1-pow(sin(wild()), 2));
3232 cout << e.expand() << endl;
3238 @node Applying a Function on Subexpressions, Polynomial Arithmetic, Pattern Matching and Advanced Substitutions, Methods and Functions
3239 @c node-name, next, previous, up
3240 @section Applying a Function on Subexpressions
3241 @cindex Tree traversal
3242 @cindex @code{map()}
3244 Sometimes you may want to perform an operation on specific parts of an
3245 expression while leaving the general structure of it intact. An example
3246 of this would be a matrix trace operation: the trace of a sum is the sum
3247 of the traces of the individual terms. That is, the trace should @dfn{map}
3248 on the sum, by applying itself to each of the sum's operands. It is possible
3249 to do this manually which usually results in code like this:
3254 if (is_a<matrix>(e))
3255 return ex_to<matrix>(e).trace();
3256 else if (is_a<add>(e)) @{
3258 for (unsigned i=0; i<e.nops(); i++)
3259 sum += calc_trace(e.op(i));
3261 @} else if (is_a<mul>)(e)) @{
3269 This is, however, slightly inefficient (if the sum is very large it can take
3270 a long time to add the terms one-by-one), and its applicability is limited to
3271 a rather small class of expressions. If @code{calc_trace()} is called with
3272 a relation or a list as its argument, you will probably want the trace to
3273 be taken on both sides of the relation or of all elements of the list.
3275 GiNaC offers the @code{map()} method to aid in the implementation of such
3279 ex ex::map(map_function & f) const;
3280 ex ex::map(ex (*f)(const ex & e)) const;
3283 In the first (preferred) form, @code{map()} takes a function object that
3284 is subclassed from the @code{map_function} class. In the second form, it
3285 takes a pointer to a function that accepts and returns an expression.
3286 @code{map()} constructs a new expression of the same type, applying the
3287 specified function on all subexpressions (in the sense of @code{op()}),
3290 The use of a function object makes it possible to supply more arguments to
3291 the function that is being mapped, or to keep local state information.
3292 The @code{map_function} class declares a virtual function call operator
3293 that you can overload. Here is a sample implementation of @code{calc_trace()}
3294 that uses @code{map()} in a recursive fashion:
3297 struct calc_trace : public map_function @{
3298 ex operator()(const ex &e)
3300 if (is_a<matrix>(e))
3301 return ex_to<matrix>(e).trace();
3302 else if (is_a<mul>(e)) @{
3305 return e.map(*this);
3310 This function object could then be used like this:
3314 ex M = ... // expression with matrices
3315 calc_trace do_trace;
3316 ex tr = do_trace(M);
3320 Here is another example for you to meditate over. It removes quadratic
3321 terms in a variable from an expanded polynomial:
3324 struct map_rem_quad : public map_function @{
3326 map_rem_quad(const ex & var_) : var(var_) @{@}
3328 ex operator()(const ex & e)
3330 if (is_a<add>(e) || is_a<mul>(e))
3331 return e.map(*this);
3332 else if (is_a<power>(e) &&
3333 e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
3343 symbol x("x"), y("y");
3346 for (int i=0; i<8; i++)
3347 e += pow(x, i) * pow(y, 8-i) * (i+1);
3349 // -> 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
3351 map_rem_quad rem_quad(x);
3352 cout << rem_quad(e) << endl;
3353 // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
3357 @command{ginsh} offers a slightly different implementation of @code{map()}
3358 that allows applying algebraic functions to operands. The second argument
3359 to @code{map()} is an expression containing the wildcard @samp{$0} which
3360 acts as the placeholder for the operands:
3365 > map(a+2*b,sin($0));
3367 > map(@{a,b,c@},$0^2+$0);
3368 @{a^2+a,b^2+b,c^2+c@}
3371 Note that it is only possible to use algebraic functions in the second
3372 argument. You can not use functions like @samp{diff()}, @samp{op()},
3373 @samp{subs()} etc. because these are evaluated immediately:
3376 > map(@{a,b,c@},diff($0,a));
3378 This is because "diff($0,a)" evaluates to "0", so the command is equivalent
3379 to "map(@{a,b,c@},0)".
3383 @node Polynomial Arithmetic, Rational Expressions, Applying a Function on Subexpressions, Methods and Functions
3384 @c node-name, next, previous, up
3385 @section Polynomial arithmetic
3387 @subsection Expanding and collecting
3388 @cindex @code{expand()}
3389 @cindex @code{collect()}
3391 A polynomial in one or more variables has many equivalent
3392 representations. Some useful ones serve a specific purpose. Consider
3393 for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
3394 21*y*z + 4*z^2} (written down here in output-style). It is equivalent
3395 to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
3396 representations are the recursive ones where one collects for exponents
3397 in one of the three variable. Since the factors are themselves
3398 polynomials in the remaining two variables the procedure can be
3399 repeated. In our example, two possibilities would be @math{(4*y + z)*x
3400 + 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
3403 To bring an expression into expanded form, its method
3409 may be called. In our example above, this corresponds to @math{4*x*y +
3410 x*z + 20*y^2 + 21*y*z + 4*z^2}. Again, since the canonical form in
3411 GiNaC is not easily guessable you should be prepared to see different
3412 orderings of terms in such sums!
3414 Another useful representation of multivariate polynomials is as a
3415 univariate polynomial in one of the variables with the coefficients
3416 being polynomials in the remaining variables. The method
3417 @code{collect()} accomplishes this task:
3420 ex ex::collect(const ex & s, bool distributed = false);
3423 The first argument to @code{collect()} can also be a list of objects in which
3424 case the result is either a recursively collected polynomial, or a polynomial
3425 in a distributed form with terms like @math{c*x1^e1*...*xn^en}, as specified
3426 by the @code{distributed} flag.
3428 Note that the original polynomial needs to be in expanded form (for the
3429 variables concerned) in order for @code{collect()} to be able to find the
3430 coefficients properly.
3432 The following @command{ginsh} transcript shows an application of @code{collect()}
3433 together with @code{find()}:
3436 > a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
3437 d*p*sin(x)+p*sin(x)+q*d*sin(x)+q*sin(y)+d*sin(x)+q*d*sin(y)+sin(y)+d*sin(y)+q*sin(x)+d*sin(y)*p+sin(x)+sin(y)*p
3438 > collect(a,@{p,q@});
3439 d*sin(x)+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*p+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*q+sin(y)+d*sin(y)+sin(x)
3440 > collect(a,find(a,sin($1)));
3441 (1+q+d+q*d+d*p+p)*sin(y)+(1+q+d+q*d+d*p+p)*sin(x)
3442 > collect(a,@{find(a,sin($1)),p,q@});
3443 (1+(1+d)*p+d+q*(1+d))*sin(x)+(1+(1+d)*p+d+q*(1+d))*sin(y)
3444 > collect(a,@{find(a,sin($1)),d@});
3445 (1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
3448 @subsection Degree and coefficients
3449 @cindex @code{degree()}
3450 @cindex @code{ldegree()}
3451 @cindex @code{coeff()}
3453 The degree and low degree of a polynomial can be obtained using the two
3457 int ex::degree(const ex & s);
3458 int ex::ldegree(const ex & s);
3461 These functions only work reliably if the input polynomial is collected in
3462 terms of the object @samp{s}. Otherwise, they are only guaranteed to return
3463 the upper/lower bounds of the exponents. If you need accurate results, you
3464 have to call @code{expand()} and/or @code{collect()} on the input polynomial.
3472 > degree(expand(a),x);
3476 @code{degree()} also works on rational functions, returning the asymptotic
3480 > degree((x+1)/(x^3+1),x);
3484 If the input is not a polynomial or rational function in the variable @samp{s},
3485 the behavior of @code{degree()} and @code{ldegree()} is undefined.
3487 To extract a coefficient with a certain power from an expanded
3491 ex ex::coeff(const ex & s, int n);
3494 You can also obtain the leading and trailing coefficients with the methods
3497 ex ex::lcoeff(const ex & s);
3498 ex ex::tcoeff(const ex & s);
3501 which are equivalent to @code{coeff(s, degree(s))} and @code{coeff(s, ldegree(s))},
3504 An application is illustrated in the next example, where a multivariate
3505 polynomial is analyzed:
3509 symbol x("x"), y("y");
3510 ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
3511 - pow(x+y,2) + 2*pow(y+2,2) - 8;
3512 ex Poly = PolyInp.expand();
3514 for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) @{
3515 cout << "The x^" << i << "-coefficient is "
3516 << Poly.coeff(x,i) << endl;
3518 cout << "As polynomial in y: "
3519 << Poly.collect(y) << endl;
3523 When run, it returns an output in the following fashion:
3526 The x^0-coefficient is y^2+11*y
3527 The x^1-coefficient is 5*y^2-2*y
3528 The x^2-coefficient is -1
3529 The x^3-coefficient is 4*y
3530 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
3533 As always, the exact output may vary between different versions of GiNaC
3534 or even from run to run since the internal canonical ordering is not
3535 within the user's sphere of influence.
3537 @code{degree()}, @code{ldegree()}, @code{coeff()}, @code{lcoeff()},
3538 @code{tcoeff()} and @code{collect()} can also be used to a certain degree
3539 with non-polynomial expressions as they not only work with symbols but with
3540 constants, functions and indexed objects as well:
3544 symbol a("a"), b("b"), c("c");
3545 idx i(symbol("i"), 3);
3547 ex e = pow(sin(x) - cos(x), 4);
3548 cout << e.degree(cos(x)) << endl;
3550 cout << e.expand().coeff(sin(x), 3) << endl;
3553 e = indexed(a+b, i) * indexed(b+c, i);
3554 e = e.expand(expand_options::expand_indexed);
3555 cout << e.collect(indexed(b, i)) << endl;
3556 // -> a.i*c.i+(a.i+c.i)*b.i+b.i^2
3561 @subsection Polynomial division
3562 @cindex polynomial division
3565 @cindex pseudo-remainder
3566 @cindex @code{quo()}
3567 @cindex @code{rem()}
3568 @cindex @code{prem()}
3569 @cindex @code{divide()}
3574 ex quo(const ex & a, const ex & b, const symbol & x);
3575 ex rem(const ex & a, const ex & b, const symbol & x);
3578 compute the quotient and remainder of univariate polynomials in the variable
3579 @samp{x}. The results satisfy @math{a = b*quo(a, b, x) + rem(a, b, x)}.
3581 The additional function
3584 ex prem(const ex & a, const ex & b, const symbol & x);
3587 computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
3588 @math{c*a = b*q + prem(a, b, x)}, where @math{c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1)}.
3590 Exact division of multivariate polynomials is performed by the function
3593 bool divide(const ex & a, const ex & b, ex & q);
3596 If @samp{b} divides @samp{a} over the rationals, this function returns @code{true}
3597 and returns the quotient in the variable @code{q}. Otherwise it returns @code{false}
3598 in which case the value of @code{q} is undefined.
3601 @subsection Unit, content and primitive part
3602 @cindex @code{unit()}
3603 @cindex @code{content()}
3604 @cindex @code{primpart()}
3609 ex ex::unit(const symbol & x);
3610 ex ex::content(const symbol & x);
3611 ex ex::primpart(const symbol & x);
3614 return the unit part, content part, and primitive polynomial of a multivariate
3615 polynomial with respect to the variable @samp{x} (the unit part being the sign
3616 of the leading coefficient, the content part being the GCD of the coefficients,
3617 and the primitive polynomial being the input polynomial divided by the unit and
3618 content parts). The product of unit, content, and primitive part is the
3619 original polynomial.
3622 @subsection GCD and LCM
3625 @cindex @code{gcd()}
3626 @cindex @code{lcm()}
3628 The functions for polynomial greatest common divisor and least common
3629 multiple have the synopsis
3632 ex gcd(const ex & a, const ex & b);
3633 ex lcm(const ex & a, const ex & b);
3636 The functions @code{gcd()} and @code{lcm()} accept two expressions
3637 @code{a} and @code{b} as arguments and return a new expression, their
3638 greatest common divisor or least common multiple, respectively. If the
3639 polynomials @code{a} and @code{b} are coprime @code{gcd(a,b)} returns 1
3640 and @code{lcm(a,b)} returns the product of @code{a} and @code{b}.
3643 #include <ginac/ginac.h>
3644 using namespace GiNaC;
3648 symbol x("x"), y("y"), z("z");
3649 ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
3650 ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
3652 ex P_gcd = gcd(P_a, P_b);
3654 ex P_lcm = lcm(P_a, P_b);
3655 // 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
3660 @subsection Square-free decomposition
3661 @cindex square-free decomposition
3662 @cindex factorization
3663 @cindex @code{sqrfree()}
3665 GiNaC still lacks proper factorization support. Some form of
3666 factorization is, however, easily implemented by noting that factors
3667 appearing in a polynomial with power two or more also appear in the
3668 derivative and hence can easily be found by computing the GCD of the
3669 original polynomial and its derivatives. Any decent system has an
3670 interface for this so called square-free factorization. So we provide
3673 ex sqrfree(const ex & a, const lst & l = lst());
3675 Here is an example that by the way illustrates how the exact form of the
3676 result may slightly depend on the order of differentiation, calling for
3677 some care with subsequent processing of the result:
3680 symbol x("x"), y("y");
3681 ex BiVarPol = expand(pow(2-2*y,3) * pow(1+x*y,2) * pow(x-2*y,2) * (x+y));
3683 cout << sqrfree(BiVarPol, lst(x,y)) << endl;
3684 // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)
3686 cout << sqrfree(BiVarPol, lst(y,x)) << endl;
3687 // -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x)
3689 cout << sqrfree(BiVarPol) << endl;
3690 // -> depending on luck, any of the above
3693 Note also, how factors with the same exponents are not fully factorized
3697 @node Rational Expressions, Symbolic Differentiation, Polynomial Arithmetic, Methods and Functions
3698 @c node-name, next, previous, up
3699 @section Rational expressions
3701 @subsection The @code{normal} method
3702 @cindex @code{normal()}
3703 @cindex simplification
3704 @cindex temporary replacement
3706 Some basic form of simplification of expressions is called for frequently.
3707 GiNaC provides the method @code{.normal()}, which converts a rational function
3708 into an equivalent rational function of the form @samp{numerator/denominator}
3709 where numerator and denominator are coprime. If the input expression is already
3710 a fraction, it just finds the GCD of numerator and denominator and cancels it,
3711 otherwise it performs fraction addition and multiplication.
3713 @code{.normal()} can also be used on expressions which are not rational functions
3714 as it will replace all non-rational objects (like functions or non-integer
3715 powers) by temporary symbols to bring the expression to the domain of rational
3716 functions before performing the normalization, and re-substituting these
3717 symbols afterwards. This algorithm is also available as a separate method
3718 @code{.to_rational()}, described below.
3720 This means that both expressions @code{t1} and @code{t2} are indeed
3721 simplified in this little code snippet:
3726 ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
3727 ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
3728 std::cout << "t1 is " << t1.normal() << std::endl;
3729 std::cout << "t2 is " << t2.normal() << std::endl;
3733 Of course this works for multivariate polynomials too, so the ratio of
3734 the sample-polynomials from the section about GCD and LCM above would be
3735 normalized to @code{P_a/P_b} = @code{(4*y+z)/(y+3*z)}.
3738 @subsection Numerator and denominator
3741 @cindex @code{numer()}
3742 @cindex @code{denom()}
3743 @cindex @code{numer_denom()}
3745 The numerator and denominator of an expression can be obtained with
3750 ex ex::numer_denom();
3753 These functions will first normalize the expression as described above and
3754 then return the numerator, denominator, or both as a list, respectively.
3755 If you need both numerator and denominator, calling @code{numer_denom()} is
3756 faster than using @code{numer()} and @code{denom()} separately.
3759 @subsection Converting to a rational expression
3760 @cindex @code{to_rational()}
3762 Some of the methods described so far only work on polynomials or rational
3763 functions. GiNaC provides a way to extend the domain of these functions to
3764 general expressions by using the temporary replacement algorithm described
3765 above. You do this by calling
3768 ex ex::to_rational(lst &l);
3771 on the expression to be converted. The supplied @code{lst} will be filled
3772 with the generated temporary symbols and their replacement expressions in
3773 a format that can be used directly for the @code{subs()} method. It can also
3774 already contain a list of replacements from an earlier application of
3775 @code{.to_rational()}, so it's possible to use it on multiple expressions
3776 and get consistent results.
3783 ex a = pow(sin(x), 2) - pow(cos(x), 2);
3784 ex b = sin(x) + cos(x);
3787 divide(a.to_rational(l), b.to_rational(l), q);
3788 cout << q.subs(l) << endl;
3792 will print @samp{sin(x)-cos(x)}.
3795 @node Symbolic Differentiation, Series Expansion, Rational Expressions, Methods and Functions
3796 @c node-name, next, previous, up
3797 @section Symbolic differentiation
3798 @cindex differentiation
3799 @cindex @code{diff()}
3801 @cindex product rule
3803 GiNaC's objects know how to differentiate themselves. Thus, a
3804 polynomial (class @code{add}) knows that its derivative is the sum of
3805 the derivatives of all the monomials:
3809 symbol x("x"), y("y"), z("z");
3810 ex P = pow(x, 5) + pow(x, 2) + y;
3812 cout << P.diff(x,2) << endl;
3814 cout << P.diff(y) << endl; // 1
3816 cout << P.diff(z) << endl; // 0
3821 If a second integer parameter @var{n} is given, the @code{diff} method
3822 returns the @var{n}th derivative.
3824 If @emph{every} object and every function is told what its derivative
3825 is, all derivatives of composed objects can be calculated using the
3826 chain rule and the product rule. Consider, for instance the expression
3827 @code{1/cosh(x)}. Since the derivative of @code{cosh(x)} is
3828 @code{sinh(x)} and the derivative of @code{pow(x,-1)} is
3829 @code{-pow(x,-2)}, GiNaC can readily compute the composition. It turns
3830 out that the composition is the generating function for Euler Numbers,
3831 i.e. the so called @var{n}th Euler number is the coefficient of
3832 @code{x^n/n!} in the expansion of @code{1/cosh(x)}. We may use this
3833 identity to code a function that generates Euler numbers in just three
3836 @cindex Euler numbers
3838 #include <ginac/ginac.h>
3839 using namespace GiNaC;
3841 ex EulerNumber(unsigned n)
3844 const ex generator = pow(cosh(x),-1);
3845 return generator.diff(x,n).subs(x==0);
3850 for (unsigned i=0; i<11; i+=2)
3851 std::cout << EulerNumber(i) << std::endl;
3856 When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
3857 @code{-61}, @code{1385}, @code{-50521}. We increment the loop variable
3858 @code{i} by two since all odd Euler numbers vanish anyways.
3861 @node Series Expansion, Symmetrization, Symbolic Differentiation, Methods and Functions
3862 @c node-name, next, previous, up
3863 @section Series expansion
3864 @cindex @code{series()}
3865 @cindex Taylor expansion
3866 @cindex Laurent expansion
3867 @cindex @code{pseries} (class)
3868 @cindex @code{Order()}
3870 Expressions know how to expand themselves as a Taylor series or (more
3871 generally) a Laurent series. As in most conventional Computer Algebra
3872 Systems, no distinction is made between those two. There is a class of
3873 its own for storing such series (@code{class pseries}) and a built-in
3874 function (called @code{Order}) for storing the order term of the series.
3875 As a consequence, if you want to work with series, i.e. multiply two
3876 series, you need to call the method @code{ex::series} again to convert
3877 it to a series object with the usual structure (expansion plus order
3878 term). A sample application from special relativity could read:
3881 #include <ginac/ginac.h>
3882 using namespace std;
3883 using namespace GiNaC;
3887 symbol v("v"), c("c");
3889 ex gamma = 1/sqrt(1 - pow(v/c,2));
3890 ex mass_nonrel = gamma.series(v==0, 10);
3892 cout << "the relativistic mass increase with v is " << endl
3893 << mass_nonrel << endl;
3895 cout << "the inverse square of this series is " << endl
3896 << pow(mass_nonrel,-2).series(v==0, 10) << endl;
3900 Only calling the series method makes the last output simplify to
3901 @math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
3902 series raised to the power @math{-2}.
3904 @cindex Machin's formula
3905 As another instructive application, let us calculate the numerical
3906 value of Archimedes' constant
3910 (for which there already exists the built-in constant @code{Pi})
3911 using Machin's amazing formula
3913 $\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
3916 @math{Pi==16*atan(1/5)-4*atan(1/239)}.
3918 We may expand the arcus tangent around @code{0} and insert the fractions
3919 @code{1/5} and @code{1/239}. But, as we have seen, a series in GiNaC
3920 carries an order term with it and the question arises what the system is
3921 supposed to do when the fractions are plugged into that order term. The
3922 solution is to use the function @code{series_to_poly()} to simply strip
3926 #include <ginac/ginac.h>
3927 using namespace GiNaC;
3929 ex machin_pi(int degr)
3932 ex pi_expansion = series_to_poly(atan(x).series(x,degr));
3933 ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
3934 -4*pi_expansion.subs(x==numeric(1,239));
3940 using std::cout; // just for fun, another way of...
3941 using std::endl; // ...dealing with this namespace std.
3943 for (int i=2; i<12; i+=2) @{
3944 pi_frac = machin_pi(i);
3945 cout << i << ":\t" << pi_frac << endl
3946 << "\t" << pi_frac.evalf() << endl;
3952 Note how we just called @code{.series(x,degr)} instead of
3953 @code{.series(x==0,degr)}. This is a simple shortcut for @code{ex}'s
3954 method @code{series()}: if the first argument is a symbol the expression
3955 is expanded in that symbol around point @code{0}. When you run this
3956 program, it will type out:
3960 3.1832635983263598326
3961 4: 5359397032/1706489875
3962 3.1405970293260603143
3963 6: 38279241713339684/12184551018734375
3964 3.141621029325034425
3965 8: 76528487109180192540976/24359780855939418203125
3966 3.141591772182177295
3967 10: 327853873402258685803048818236/104359128170408663038552734375
3968 3.1415926824043995174
3972 @node Symmetrization, Built-in Functions, Series Expansion, Methods and Functions
3973 @c node-name, next, previous, up
3974 @section Symmetrization
3975 @cindex @code{symmetrize()}
3976 @cindex @code{antisymmetrize()}
3977 @cindex @code{symmetrize_cyclic()}
3982 ex ex::symmetrize(const lst & l);
3983 ex ex::antisymmetrize(const lst & l);
3984 ex ex::symmetrize_cyclic(const lst & l);
3987 symmetrize an expression by returning the sum over all symmetric,
3988 antisymmetric or cyclic permutations of the specified list of objects,
3989 weighted by the number of permutations.
3991 The three additional methods
3994 ex ex::symmetrize();
3995 ex ex::antisymmetrize();
3996 ex ex::symmetrize_cyclic();
3999 symmetrize or antisymmetrize an expression over its free indices.
4001 Symmetrization is most useful with indexed expressions but can be used with
4002 almost any kind of object (anything that is @code{subs()}able):
4006 idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
4007 symbol A("A"), B("B"), a("a"), b("b"), c("c");
4009 cout << indexed(A, i, j).symmetrize() << endl;
4010 // -> 1/2*A.j.i+1/2*A.i.j
4011 cout << indexed(A, i, j, k).antisymmetrize(lst(i, j)) << endl;
4012 // -> -1/2*A.j.i.k+1/2*A.i.j.k
4013 cout << lst(a, b, c).symmetrize_cyclic(lst(a, b, c)) << endl;
4014 // -> 1/3*@{a,b,c@}+1/3*@{b,c,a@}+1/3*@{c,a,b@}
4019 @node Built-in Functions, Input/Output, Symmetrization, Methods and Functions
4020 @c node-name, next, previous, up
4021 @section Predefined mathematical functions
4023 GiNaC contains the following predefined mathematical functions:
4026 @multitable @columnfractions .30 .70
4027 @item @strong{Name} @tab @strong{Function}
4030 @cindex @code{abs()}
4031 @item @code{csgn(x)}
4033 @cindex @code{csgn()}
4034 @item @code{sqrt(x)}
4035 @tab square root (not a GiNaC function, rather an alias for @code{pow(x, numeric(1, 2))})
4036 @cindex @code{sqrt()}
4039 @cindex @code{sin()}
4042 @cindex @code{cos()}
4045 @cindex @code{tan()}
4046 @item @code{asin(x)}
4048 @cindex @code{asin()}
4049 @item @code{acos(x)}
4051 @cindex @code{acos()}
4052 @item @code{atan(x)}
4053 @tab inverse tangent
4054 @cindex @code{atan()}
4055 @item @code{atan2(y, x)}
4056 @tab inverse tangent with two arguments
4057 @item @code{sinh(x)}
4058 @tab hyperbolic sine
4059 @cindex @code{sinh()}
4060 @item @code{cosh(x)}
4061 @tab hyperbolic cosine
4062 @cindex @code{cosh()}
4063 @item @code{tanh(x)}
4064 @tab hyperbolic tangent
4065 @cindex @code{tanh()}
4066 @item @code{asinh(x)}
4067 @tab inverse hyperbolic sine
4068 @cindex @code{asinh()}
4069 @item @code{acosh(x)}
4070 @tab inverse hyperbolic cosine
4071 @cindex @code{acosh()}
4072 @item @code{atanh(x)}
4073 @tab inverse hyperbolic tangent
4074 @cindex @code{atanh()}
4076 @tab exponential function
4077 @cindex @code{exp()}
4079 @tab natural logarithm
4080 @cindex @code{log()}
4083 @cindex @code{Li2()}
4084 @item @code{zeta(x)}
4085 @tab Riemann's zeta function
4086 @cindex @code{zeta()}
4087 @item @code{zeta(n, x)}
4088 @tab derivatives of Riemann's zeta function
4089 @item @code{tgamma(x)}
4091 @cindex @code{tgamma()}
4092 @cindex Gamma function
4093 @item @code{lgamma(x)}
4094 @tab logarithm of Gamma function
4095 @cindex @code{lgamma()}
4096 @item @code{beta(x, y)}
4097 @tab Beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
4098 @cindex @code{beta()}
4100 @tab psi (digamma) function
4101 @cindex @code{psi()}
4102 @item @code{psi(n, x)}
4103 @tab derivatives of psi function (polygamma functions)
4104 @item @code{factorial(n)}
4105 @tab factorial function
4106 @cindex @code{factorial()}
4107 @item @code{binomial(n, m)}
4108 @tab binomial coefficients
4109 @cindex @code{binomial()}
4110 @item @code{Order(x)}
4111 @tab order term function in truncated power series
4112 @cindex @code{Order()}
4117 For functions that have a branch cut in the complex plane GiNaC follows
4118 the conventions for C++ as defined in the ANSI standard as far as
4119 possible. In particular: the natural logarithm (@code{log}) and the
4120 square root (@code{sqrt}) both have their branch cuts running along the
4121 negative real axis where the points on the axis itself belong to the
4122 upper part (i.e. continuous with quadrant II). The inverse
4123 trigonometric and hyperbolic functions are not defined for complex
4124 arguments by the C++ standard, however. In GiNaC we follow the
4125 conventions used by CLN, which in turn follow the carefully designed
4126 definitions in the Common Lisp standard. It should be noted that this
4127 convention is identical to the one used by the C99 standard and by most
4128 serious CAS. It is to be expected that future revisions of the C++
4129 standard incorporate these functions in the complex domain in a manner
4130 compatible with C99.
4133 @node Input/Output, Extending GiNaC, Built-in Functions, Methods and Functions
4134 @c node-name, next, previous, up
4135 @section Input and output of expressions
4138 @subsection Expression output
4140 @cindex output of expressions
4142 The easiest way to print an expression is to write it to a stream:
4147 ex e = 4.5+pow(x,2)*3/2;
4148 cout << e << endl; // prints '(4.5)+3/2*x^2'
4152 The output format is identical to the @command{ginsh} input syntax and
4153 to that used by most computer algebra systems, but not directly pastable
4154 into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
4155 is printed as @samp{x^2}).
4157 It is possible to print expressions in a number of different formats with
4161 void ex::print(const print_context & c, unsigned level = 0);
4164 @cindex @code{print_context} (class)
4165 The type of @code{print_context} object passed in determines the format
4166 of the output. The possible types are defined in @file{ginac/print.h}.
4167 All constructors of @code{print_context} and derived classes take an
4168 @code{ostream &} as their first argument.
4170 To print an expression in a way that can be directly used in a C or C++
4171 program, you pass a @code{print_csrc} object like this:
4175 cout << "float f = ";
4176 e.print(print_csrc_float(cout));
4179 cout << "double d = ";
4180 e.print(print_csrc_double(cout));
4183 cout << "cl_N n = ";
4184 e.print(print_csrc_cl_N(cout));
4189 The three possible types mostly affect the way in which floating point
4190 numbers are written.
4192 The above example will produce (note the @code{x^2} being converted to @code{x*x}):
4195 float f = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
4196 double d = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
4197 cl_N n = (cln::cl_F("3.0")/cln::cl_F("2.0"))*(x*x)+cln::cl_F("4.5");
4200 The @code{print_context} type @code{print_tree} provides a dump of the
4201 internal structure of an expression for debugging purposes:
4205 e.print(print_tree(cout));
4212 add, hash=0x0, flags=0x3, nops=2
4213 power, hash=0x9, flags=0x3, nops=2
4214 x (symbol), serial=3, hash=0x44a113a6, flags=0xf
4215 2 (numeric), hash=0x80000042, flags=0xf
4216 3/2 (numeric), hash=0x80000061, flags=0xf
4219 4.5L0 (numeric), hash=0x8000004b, flags=0xf
4223 This kind of output is also available in @command{ginsh} as the @code{print()}
4226 Another useful output format is for LaTeX parsing in mathematical mode.
4227 It is rather similar to the default @code{print_context} but provides
4228 some braces needed by LaTeX for delimiting boxes and also converts some
4229 common objects to conventional LaTeX names. It is possible to give symbols
4230 a special name for LaTeX output by supplying it as a second argument to
4231 the @code{symbol} constructor.
4233 For example, the code snippet
4238 ex foo = lgamma(x).series(x==0,3);
4239 foo.print(print_latex(std::cout));
4245 @{(-\ln(x))@}+@{(-\gamma_E)@} x+@{(1/12 \pi^2)@} x^@{2@}+\mathcal@{O@}(x^3)
4248 @cindex Tree traversal
4249 If you need any fancy special output format, e.g. for interfacing GiNaC
4250 with other algebra systems or for producing code for different
4251 programming languages, you can always traverse the expression tree yourself:
4254 static void my_print(const ex & e)
4256 if (is_a<function>(e))
4257 cout << ex_to<function>(e).get_name();
4259 cout << e.bp->class_name();
4261 unsigned n = e.nops();
4263 for (unsigned i=0; i<n; i++) @{
4275 my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
4283 add(power(numeric(3),symbol(x)),mul(sin(mul(power(constant(Pi),numeric(-1)),
4284 symbol(y))),numeric(-2)))
4287 If you need an output format that makes it possible to accurately
4288 reconstruct an expression by feeding the output to a suitable parser or
4289 object factory, you should consider storing the expression in an
4290 @code{archive} object and reading the object properties from there.
4291 See the section on archiving for more information.
4294 @subsection Expression input
4295 @cindex input of expressions
4297 GiNaC provides no way to directly read an expression from a stream because
4298 you will usually want the user to be able to enter something like @samp{2*x+sin(y)}
4299 and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
4300 @code{y} you defined in your program and there is no way to specify the
4301 desired symbols to the @code{>>} stream input operator.
4303 Instead, GiNaC lets you construct an expression from a string, specifying the
4304 list of symbols and indices to be used:
4308 symbol x("x"), y("y"), p("p");
4309 idx i(symbol("i"), 3);
4310 ex e("2*x+sin(y)+p.i", lst(x, y, p, i));
4314 The input syntax is the same as that used by @command{ginsh} and the stream
4315 output operator @code{<<}. The symbols and indices in the string are matched
4316 by name to the symbols and indices in the list and if GiNaC encounters a
4317 symbol or index not specified in the list it will throw an exception. Only
4318 indices whose values are single symbols can be used (i.e. numeric indices
4319 or compound indices as in "A.(2*n+1)" are not allowed).
4321 With this constructor, it's also easy to implement interactive GiNaC programs:
4326 #include <stdexcept>
4327 #include <ginac/ginac.h>
4328 using namespace std;
4329 using namespace GiNaC;
4336 cout << "Enter an expression containing 'x': ";
4341 cout << "The derivative of " << e << " with respect to x is ";
4342 cout << e.diff(x) << ".\n";
4343 @} catch (exception &p) @{
4344 cerr << p.what() << endl;
4350 @subsection Archiving
4351 @cindex @code{archive} (class)
4354 GiNaC allows creating @dfn{archives} of expressions which can be stored
4355 to or retrieved from files. To create an archive, you declare an object
4356 of class @code{archive} and archive expressions in it, giving each
4357 expression a unique name:
4361 using namespace std;
4362 #include <ginac/ginac.h>
4363 using namespace GiNaC;
4367 symbol x("x"), y("y"), z("z");
4369 ex foo = sin(x + 2*y) + 3*z + 41;
4373 a.archive_ex(foo, "foo");
4374 a.archive_ex(bar, "the second one");
4378 The archive can then be written to a file:
4382 ofstream out("foobar.gar");
4388 The file @file{foobar.gar} contains all information that is needed to
4389 reconstruct the expressions @code{foo} and @code{bar}.
4391 @cindex @command{viewgar}
4392 The tool @command{viewgar} that comes with GiNaC can be used to view
4393 the contents of GiNaC archive files:
4396 $ viewgar foobar.gar
4397 foo = 41+sin(x+2*y)+3*z
4398 the second one = 42+sin(x+2*y)+3*z
4401 The point of writing archive files is of course that they can later be
4407 ifstream in("foobar.gar");
4412 And the stored expressions can be retrieved by their name:
4418 ex ex1 = a2.unarchive_ex(syms, "foo");
4419 ex ex2 = a2.unarchive_ex(syms, "the second one");
4421 cout << ex1 << endl; // prints "41+sin(x+2*y)+3*z"
4422 cout << ex2 << endl; // prints "42+sin(x+2*y)+3*z"
4423 cout << ex1.subs(x == 2) << endl; // prints "41+sin(2+2*y)+3*z"
4427 Note that you have to supply a list of the symbols which are to be inserted
4428 in the expressions. Symbols in archives are stored by their name only and
4429 if you don't specify which symbols you have, unarchiving the expression will
4430 create new symbols with that name. E.g. if you hadn't included @code{x} in
4431 the @code{syms} list above, the @code{ex1.subs(x == 2)} statement would
4432 have had no effect because the @code{x} in @code{ex1} would have been a
4433 different symbol than the @code{x} which was defined at the beginning of
4434 the program, although both would appear as @samp{x} when printed.
4436 You can also use the information stored in an @code{archive} object to
4437 output expressions in a format suitable for exact reconstruction. The
4438 @code{archive} and @code{archive_node} classes have a couple of member
4439 functions that let you access the stored properties:
4442 static void my_print2(const archive_node & n)
4445 n.find_string("class", class_name);
4446 cout << class_name << "(";
4448 archive_node::propinfovector p;
4449 n.get_properties(p);
4451 unsigned num = p.size();
4452 for (unsigned i=0; i<num; i++) @{
4453 const string &name = p[i].name;
4454 if (name == "class")
4456 cout << name << "=";
4458 unsigned count = p[i].count;
4462 for (unsigned j=0; j<count; j++) @{
4463 switch (p[i].type) @{
4464 case archive_node::PTYPE_BOOL: @{
4466 n.find_bool(name, x, j);
4467 cout << (x ? "true" : "false");
4470 case archive_node::PTYPE_UNSIGNED: @{
4472 n.find_unsigned(name, x, j);
4476 case archive_node::PTYPE_STRING: @{
4478 n.find_string(name, x, j);
4479 cout << '\"' << x << '\"';
4482 case archive_node::PTYPE_NODE: @{
4483 const archive_node &x = n.find_ex_node(name, j);
4505 ex e = pow(2, x) - y;
4507 my_print2(ar.get_top_node(0)); cout << endl;
4515 add(rest=@{power(basis=numeric(number="2"),exponent=symbol(name="x")),
4516 symbol(name="y")@},coeff=@{numeric(number="1"),numeric(number="-1")@},
4517 overall_coeff=numeric(number="0"))
4520 Be warned, however, that the set of properties and their meaning for each
4521 class may change between GiNaC versions.
4524 @node Extending GiNaC, What does not belong into GiNaC, Input/Output, Top
4525 @c node-name, next, previous, up
4526 @chapter Extending GiNaC
4528 By reading so far you should have gotten a fairly good understanding of
4529 GiNaC's design-patterns. From here on you should start reading the
4530 sources. All we can do now is issue some recommendations how to tackle
4531 GiNaC's many loose ends in order to fulfill everybody's dreams. If you
4532 develop some useful extension please don't hesitate to contact the GiNaC
4533 authors---they will happily incorporate them into future versions.
4536 * What does not belong into GiNaC:: What to avoid.
4537 * Symbolic functions:: Implementing symbolic functions.
4538 * Adding classes:: Defining new algebraic classes.
4542 @node What does not belong into GiNaC, Symbolic functions, Extending GiNaC, Extending GiNaC
4543 @c node-name, next, previous, up
4544 @section What doesn't belong into GiNaC
4546 @cindex @command{ginsh}
4547 First of all, GiNaC's name must be read literally. It is designed to be
4548 a library for use within C++. The tiny @command{ginsh} accompanying
4549 GiNaC makes this even more clear: it doesn't even attempt to provide a
4550 language. There are no loops or conditional expressions in
4551 @command{ginsh}, it is merely a window into the library for the
4552 programmer to test stuff (or to show off). Still, the design of a
4553 complete CAS with a language of its own, graphical capabilities and all
4554 this on top of GiNaC is possible and is without doubt a nice project for
4557 There are many built-in functions in GiNaC that do not know how to
4558 evaluate themselves numerically to a precision declared at runtime
4559 (using @code{Digits}). Some may be evaluated at certain points, but not
4560 generally. This ought to be fixed. However, doing numerical
4561 computations with GiNaC's quite abstract classes is doomed to be
4562 inefficient. For this purpose, the underlying foundation classes
4563 provided by CLN are much better suited.
4566 @node Symbolic functions, Adding classes, What does not belong into GiNaC, Extending GiNaC
4567 @c node-name, next, previous, up
4568 @section Symbolic functions
4570 The easiest and most instructive way to start with is probably to
4571 implement your own function. GiNaC's functions are objects of class
4572 @code{function}. The preprocessor is then used to convert the function
4573 names to objects with a corresponding serial number that is used
4574 internally to identify them. You usually need not worry about this
4575 number. New functions may be inserted into the system via a kind of
4576 `registry'. It is your responsibility to care for some functions that
4577 are called when the user invokes certain methods. These are usual
4578 C++-functions accepting a number of @code{ex} as arguments and returning
4579 one @code{ex}. As an example, if we have a look at a simplified
4580 implementation of the cosine trigonometric function, we first need a
4581 function that is called when one wishes to @code{eval} it. It could
4582 look something like this:
4585 static ex cos_eval_method(const ex & x)
4587 // if (!x%(2*Pi)) return 1
4588 // if (!x%Pi) return -1
4589 // if (!x%Pi/2) return 0
4590 // care for other cases...
4591 return cos(x).hold();
4595 @cindex @code{hold()}
4597 The last line returns @code{cos(x)} if we don't know what else to do and
4598 stops a potential recursive evaluation by saying @code{.hold()}, which
4599 sets a flag to the expression signaling that it has been evaluated. We
4600 should also implement a method for numerical evaluation and since we are
4601 lazy we sweep the problem under the rug by calling someone else's
4602 function that does so, in this case the one in class @code{numeric}:
4605 static ex cos_evalf(const ex & x)
4607 if (is_a<numeric>(x))
4608 return cos(ex_to<numeric>(x));
4610 return cos(x).hold();
4614 Differentiation will surely turn up and so we need to tell @code{cos}
4615 what the first derivative is (higher derivatives (@code{.diff(x,3)} for
4616 instance are then handled automatically by @code{basic::diff} and
4620 static ex cos_deriv(const ex & x, unsigned diff_param)
4626 @cindex product rule
4627 The second parameter is obligatory but uninteresting at this point. It
4628 specifies which parameter to differentiate in a partial derivative in
4629 case the function has more than one parameter and its main application
4630 is for correct handling of the chain rule. For Taylor expansion, it is
4631 enough to know how to differentiate. But if the function you want to
4632 implement does have a pole somewhere in the complex plane, you need to
4633 write another method for Laurent expansion around that point.
4635 Now that all the ingredients for @code{cos} have been set up, we need
4636 to tell the system about it. This is done by a macro and we are not
4637 going to describe how it expands, please consult your preprocessor if you
4641 REGISTER_FUNCTION(cos, eval_func(cos_eval).
4642 evalf_func(cos_evalf).
4643 derivative_func(cos_deriv));
4646 The first argument is the function's name used for calling it and for
4647 output. The second binds the corresponding methods as options to this
4648 object. Options are separated by a dot and can be given in an arbitrary
4649 order. GiNaC functions understand several more options which are always
4650 specified as @code{.option(params)}, for example a method for series
4651 expansion @code{.series_func(cos_series)}. Again, if no series
4652 expansion method is given, GiNaC defaults to simple Taylor expansion,
4653 which is correct if there are no poles involved as is the case for the
4654 @code{cos} function. The way GiNaC handles poles in case there are any
4655 is best understood by studying one of the examples, like the Gamma
4656 (@code{tgamma}) function for instance. (In essence the function first
4657 checks if there is a pole at the evaluation point and falls back to
4658 Taylor expansion if there isn't. Then, the pole is regularized by some
4659 suitable transformation.) Also, the new function needs to be declared
4660 somewhere. This may also be done by a convenient preprocessor macro:
4663 DECLARE_FUNCTION_1P(cos)
4666 The suffix @code{_1P} stands for @emph{one parameter}. Of course, this
4667 implementation of @code{cos} is very incomplete and lacks several safety
4668 mechanisms. Please, have a look at the real implementation in GiNaC.
4669 (By the way: in case you are worrying about all the macros above we can
4670 assure you that functions are GiNaC's most macro-intense classes. We
4671 have done our best to avoid macros where we can.)
4674 @node Adding classes, A Comparison With Other CAS, Symbolic functions, Extending GiNaC
4675 @c node-name, next, previous, up
4676 @section Adding classes
4678 If you are doing some very specialized things with GiNaC you may find that
4679 you have to implement your own algebraic classes to fit your needs. This
4680 section will explain how to do this by giving the example of a simple
4681 'string' class. After reading this section you will know how to properly
4682 declare a GiNaC class and what the minimum required member functions are
4683 that you have to implement. We only cover the implementation of a 'leaf'
4684 class here (i.e. one that doesn't contain subexpressions). Creating a
4685 container class like, for example, a class representing tensor products is
4686 more involved but this section should give you enough information so you can
4687 consult the source to GiNaC's predefined classes if you want to implement
4688 something more complicated.
4690 @subsection GiNaC's run-time type information system
4692 @cindex hierarchy of classes
4694 All algebraic classes (that is, all classes that can appear in expressions)
4695 in GiNaC are direct or indirect subclasses of the class @code{basic}. So a
4696 @code{basic *} (which is essentially what an @code{ex} is) represents a
4697 generic pointer to an algebraic class. Occasionally it is necessary to find
4698 out what the class of an object pointed to by a @code{basic *} really is.
4699 Also, for the unarchiving of expressions it must be possible to find the
4700 @code{unarchive()} function of a class given the class name (as a string). A
4701 system that provides this kind of information is called a run-time type
4702 information (RTTI) system. The C++ language provides such a thing (see the
4703 standard header file @file{<typeinfo>}) but for efficiency reasons GiNaC
4704 implements its own, simpler RTTI.
4706 The RTTI in GiNaC is based on two mechanisms:
4711 The @code{basic} class declares a member variable @code{tinfo_key} which
4712 holds an unsigned integer that identifies the object's class. These numbers
4713 are defined in the @file{tinfos.h} header file for the built-in GiNaC
4714 classes. They all start with @code{TINFO_}.
4717 By means of some clever tricks with static members, GiNaC maintains a list
4718 of information for all classes derived from @code{basic}. The information
4719 available includes the class names, the @code{tinfo_key}s, and pointers
4720 to the unarchiving functions. This class registry is defined in the
4721 @file{registrar.h} header file.
4725 The disadvantage of this proprietary RTTI implementation is that there's
4726 a little more to do when implementing new classes (C++'s RTTI works more
4727 or less automatic) but don't worry, most of the work is simplified by
4730 @subsection A minimalistic example
4732 Now we will start implementing a new class @code{mystring} that allows
4733 placing character strings in algebraic expressions (this is not very useful,
4734 but it's just an example). This class will be a direct subclass of
4735 @code{basic}. You can use this sample implementation as a starting point
4736 for your own classes.
4738 The code snippets given here assume that you have included some header files
4744 #include <stdexcept>
4745 using namespace std;
4747 #include <ginac/ginac.h>
4748 using namespace GiNaC;
4751 The first thing we have to do is to define a @code{tinfo_key} for our new
4752 class. This can be any arbitrary unsigned number that is not already taken
4753 by one of the existing classes but it's better to come up with something
4754 that is unlikely to clash with keys that might be added in the future. The
4755 numbers in @file{tinfos.h} are modeled somewhat after the class hierarchy
4756 which is not a requirement but we are going to stick with this scheme:
4759 const unsigned TINFO_mystring = 0x42420001U;
4762 Now we can write down the class declaration. The class stores a C++
4763 @code{string} and the user shall be able to construct a @code{mystring}
4764 object from a C or C++ string:
4767 class mystring : public basic
4769 GINAC_DECLARE_REGISTERED_CLASS(mystring, basic)
4772 mystring(const string &s);
4773 mystring(const char *s);
4779 GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
4782 The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
4783 macros are defined in @file{registrar.h}. They take the name of the class
4784 and its direct superclass as arguments and insert all required declarations
4785 for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
4786 the first line after the opening brace of the class definition. The
4787 @code{GINAC_IMPLEMENT_REGISTERED_CLASS} may appear anywhere else in the
4788 source (at global scope, of course, not inside a function).
4790 @code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
4791 declarations of the default and copy constructor, the destructor, the
4792 assignment operator and a couple of other functions that are required. It
4793 also defines a type @code{inherited} which refers to the superclass so you
4794 don't have to modify your code every time you shuffle around the class
4795 hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} implements the copy
4796 constructor, the destructor and the assignment operator.
4798 Now there are nine member functions we have to implement to get a working
4804 @code{mystring()}, the default constructor.
4807 @code{void destroy(bool call_parent)}, which is used in the destructor and the
4808 assignment operator to free dynamically allocated members. The @code{call_parent}
4809 specifies whether the @code{destroy()} function of the superclass is to be
4813 @code{void copy(const mystring &other)}, which is used in the copy constructor
4814 and assignment operator to copy the member variables over from another
4815 object of the same class.
4818 @code{void archive(archive_node &n)}, the archiving function. This stores all
4819 information needed to reconstruct an object of this class inside an
4820 @code{archive_node}.
4823 @code{mystring(const archive_node &n, const lst &sym_lst)}, the unarchiving
4824 constructor. This constructs an instance of the class from the information
4825 found in an @code{archive_node}.
4828 @code{ex unarchive(const archive_node &n, const lst &sym_lst)}, the static
4829 unarchiving function. It constructs a new instance by calling the unarchiving
4833 @code{int compare_same_type(const basic &other)}, which is used internally
4834 by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
4835 -1, depending on the relative order of this object and the @code{other}
4836 object. If it returns 0, the objects are considered equal.
4837 @strong{Note:} This has nothing to do with the (numeric) ordering
4838 relationship expressed by @code{<}, @code{>=} etc (which cannot be defined
4839 for non-numeric classes). For example, @code{numeric(1).compare_same_type(numeric(2))}
4840 may return +1 even though 1 is clearly smaller than 2. Every GiNaC class
4841 must provide a @code{compare_same_type()} function, even those representing
4842 objects for which no reasonable algebraic ordering relationship can be
4846 And, of course, @code{mystring(const string &s)} and @code{mystring(const char *s)}
4847 which are the two constructors we declared.
4851 Let's proceed step-by-step. The default constructor looks like this:
4854 mystring::mystring() : inherited(TINFO_mystring)
4856 // dynamically allocate resources here if required
4860 The golden rule is that in all constructors you have to set the
4861 @code{tinfo_key} member to the @code{TINFO_*} value of your class. Otherwise
4862 it will be set by the constructor of the superclass and all hell will break
4863 loose in the RTTI. For your convenience, the @code{basic} class provides
4864 a constructor that takes a @code{tinfo_key} value, which we are using here
4865 (remember that in our case @code{inherited = basic}). If the superclass
4866 didn't have such a constructor, we would have to set the @code{tinfo_key}
4867 to the right value manually.
4869 In the default constructor you should set all other member variables to
4870 reasonable default values (we don't need that here since our @code{str}
4871 member gets set to an empty string automatically). The constructor(s) are of
4872 course also the right place to allocate any dynamic resources you require.
4874 Next, the @code{destroy()} function:
4877 void mystring::destroy(bool call_parent)
4879 // free dynamically allocated resources here if required
4881 inherited::destroy(call_parent);
4885 This function is where we free all dynamically allocated resources. We
4886 don't have any so we're not doing anything here, but if we had, for
4887 example, used a C-style @code{char *} to store our string, this would be
4888 the place to @code{delete[]} the string storage. If @code{call_parent}
4889 is true, we have to call the @code{destroy()} function of the superclass
4890 after we're done (to mimic C++'s automatic invocation of superclass
4891 destructors where @code{destroy()} is called from outside a destructor).
4893 The @code{copy()} function just copies over the member variables from
4897 void mystring::copy(const mystring &other)
4899 inherited::copy(other);
4904 We can simply overwrite the member variables here. There's no need to worry
4905 about dynamically allocated storage. The assignment operator (which is
4906 automatically defined by @code{GINAC_IMPLEMENT_REGISTERED_CLASS}, as you
4907 recall) calls @code{destroy()} before it calls @code{copy()}. You have to
4908 explicitly call the @code{copy()} function of the superclass here so
4909 all the member variables will get copied.
4911 Next are the three functions for archiving. You have to implement them even
4912 if you don't plan to use archives, but the minimum required implementation
4913 is really simple. First, the archiving function:
4916 void mystring::archive(archive_node &n) const
4918 inherited::archive(n);
4919 n.add_string("string", str);
4923 The only thing that is really required is calling the @code{archive()}
4924 function of the superclass. Optionally, you can store all information you
4925 deem necessary for representing the object into the passed
4926 @code{archive_node}. We are just storing our string here. For more
4927 information on how the archiving works, consult the @file{archive.h} header
4930 The unarchiving constructor is basically the inverse of the archiving
4934 mystring::mystring(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
4936 n.find_string("string", str);
4940 If you don't need archiving, just leave this function empty (but you must
4941 invoke the unarchiving constructor of the superclass). Note that we don't
4942 have to set the @code{tinfo_key} here because it is done automatically
4943 by the unarchiving constructor of the @code{basic} class.
4945 Finally, the unarchiving function:
4948 ex mystring::unarchive(const archive_node &n, const lst &sym_lst)
4950 return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
4954 You don't have to understand how exactly this works. Just copy these
4955 four lines into your code literally (replacing the class name, of
4956 course). It calls the unarchiving constructor of the class and unless
4957 you are doing something very special (like matching @code{archive_node}s
4958 to global objects) you don't need a different implementation. For those
4959 who are interested: setting the @code{dynallocated} flag puts the object
4960 under the control of GiNaC's garbage collection. It will get deleted
4961 automatically once it is no longer referenced.
4963 Our @code{compare_same_type()} function uses a provided function to compare
4967 int mystring::compare_same_type(const basic &other) const
4969 const mystring &o = static_cast<const mystring &>(other);
4970 int cmpval = str.compare(o.str);
4973 else if (cmpval < 0)
4980 Although this function takes a @code{basic &}, it will always be a reference
4981 to an object of exactly the same class (objects of different classes are not
4982 comparable), so the cast is safe. If this function returns 0, the two objects
4983 are considered equal (in the sense that @math{A-B=0}), so you should compare
4984 all relevant member variables.
4986 Now the only thing missing is our two new constructors:
4989 mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s)
4991 // dynamically allocate resources here if required
4994 mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s)
4996 // dynamically allocate resources here if required
5000 No surprises here. We set the @code{str} member from the argument and
5001 remember to pass the right @code{tinfo_key} to the @code{basic} constructor.
5003 That's it! We now have a minimal working GiNaC class that can store
5004 strings in algebraic expressions. Let's confirm that the RTTI works:
5007 ex e = mystring("Hello, world!");
5008 cout << is_a<mystring>(e) << endl;
5011 cout << e.bp->class_name() << endl;
5015 Obviously it does. Let's see what the expression @code{e} looks like:
5019 // -> [mystring object]
5022 Hm, not exactly what we expect, but of course the @code{mystring} class
5023 doesn't yet know how to print itself. This is done in the @code{print()}
5024 member function. Let's say that we wanted to print the string surrounded
5028 class mystring : public basic
5032 void print(const print_context &c, unsigned level = 0) const;
5036 void mystring::print(const print_context &c, unsigned level) const
5038 // print_context::s is a reference to an ostream
5039 c.s << '\"' << str << '\"';
5043 The @code{level} argument is only required for container classes to
5044 correctly parenthesize the output. Let's try again to print the expression:
5048 // -> "Hello, world!"
5051 Much better. The @code{mystring} class can be used in arbitrary expressions:
5054 e += mystring("GiNaC rulez");
5056 // -> "GiNaC rulez"+"Hello, world!"
5059 (GiNaC's automatic term reordering is in effect here), or even
5062 e = pow(mystring("One string"), 2*sin(Pi-mystring("Another string")));
5064 // -> "One string"^(2*sin(-"Another string"+Pi))
5067 Whether this makes sense is debatable but remember that this is only an
5068 example. At least it allows you to implement your own symbolic algorithms
5071 Note that GiNaC's algebraic rules remain unchanged:
5074 e = mystring("Wow") * mystring("Wow");
5078 e = pow(mystring("First")-mystring("Second"), 2);
5079 cout << e.expand() << endl;
5080 // -> -2*"First"*"Second"+"First"^2+"Second"^2
5083 There's no way to, for example, make GiNaC's @code{add} class perform string
5084 concatenation. You would have to implement this yourself.
5086 @subsection Automatic evaluation
5088 @cindex @code{hold()}
5089 @cindex @code{eval()}
5091 When dealing with objects that are just a little more complicated than the
5092 simple string objects we have implemented, chances are that you will want to
5093 have some automatic simplifications or canonicalizations performed on them.
5094 This is done in the evaluation member function @code{eval()}. Let's say that
5095 we wanted all strings automatically converted to lowercase with
5096 non-alphabetic characters stripped, and empty strings removed:
5099 class mystring : public basic
5103 ex eval(int level = 0) const;
5107 ex mystring::eval(int level) const
5110 for (int i=0; i<str.length(); i++) @{
5112 if (c >= 'A' && c <= 'Z')
5113 new_str += tolower(c);
5114 else if (c >= 'a' && c <= 'z')
5118 if (new_str.length() == 0)
5121 return mystring(new_str).hold();
5125 The @code{level} argument is used to limit the recursion depth of the
5126 evaluation. We don't have any subexpressions in the @code{mystring}
5127 class so we are not concerned with this. If we had, we would call the
5128 @code{eval()} functions of the subexpressions with @code{level - 1} as
5129 the argument if @code{level != 1}. The @code{hold()} member function
5130 sets a flag in the object that prevents further evaluation. Otherwise
5131 we might end up in an endless loop. When you want to return the object
5132 unmodified, use @code{return this->hold();}.
5134 Let's confirm that it works:
5137 ex e = mystring("Hello, world!") + mystring("!?#");
5141 e = mystring("Wow!") + mystring("WOW") + mystring(" W ** o ** W");
5146 @subsection Other member functions
5148 We have implemented only a small set of member functions to make the class
5149 work in the GiNaC framework. For a real algebraic class, there are probably
5150 some more functions that you will want to re-implement, such as
5151 @code{evalf()}, @code{series()} or @code{op()}. Have a look at @file{basic.h}
5152 or the header file of the class you want to make a subclass of to see
5153 what's there. One member function that you will most likely want to
5154 implement for terminal classes like the described string class is
5155 @code{calcchash()} that returns an @code{unsigned} hash value for the object
5156 which will allow GiNaC to compare and canonicalize expressions much more
5159 You can, of course, also add your own new member functions. Remember,
5160 that the RTTI may be used to get information about what kinds of objects
5161 you are dealing with (the position in the class hierarchy) and that you
5162 can always extract the bare object from an @code{ex} by stripping the
5163 @code{ex} off using the @code{ex_to<mystring>(e)} function when that
5164 should become a need.
5166 That's it. May the source be with you!
5169 @node A Comparison With Other CAS, Advantages, Adding classes, Top
5170 @c node-name, next, previous, up
5171 @chapter A Comparison With Other CAS
5174 This chapter will give you some information on how GiNaC compares to
5175 other, traditional Computer Algebra Systems, like @emph{Maple},
5176 @emph{Mathematica} or @emph{Reduce}, where it has advantages and
5177 disadvantages over these systems.
5180 * Advantages:: Strengths of the GiNaC approach.
5181 * Disadvantages:: Weaknesses of the GiNaC approach.
5182 * Why C++?:: Attractiveness of C++.
5185 @node Advantages, Disadvantages, A Comparison With Other CAS, A Comparison With Other CAS
5186 @c node-name, next, previous, up
5189 GiNaC has several advantages over traditional Computer
5190 Algebra Systems, like
5195 familiar language: all common CAS implement their own proprietary
5196 grammar which you have to learn first (and maybe learn again when your
5197 vendor decides to `enhance' it). With GiNaC you can write your program
5198 in common C++, which is standardized.
5202 structured data types: you can build up structured data types using
5203 @code{struct}s or @code{class}es together with STL features instead of
5204 using unnamed lists of lists of lists.
5207 strongly typed: in CAS, you usually have only one kind of variables
5208 which can hold contents of an arbitrary type. This 4GL like feature is
5209 nice for novice programmers, but dangerous.
5212 development tools: powerful development tools exist for C++, like fancy
5213 editors (e.g. with automatic indentation and syntax highlighting),
5214 debuggers, visualization tools, documentation generators@dots{}
5217 modularization: C++ programs can easily be split into modules by
5218 separating interface and implementation.
5221 price: GiNaC is distributed under the GNU Public License which means
5222 that it is free and available with source code. And there are excellent
5223 C++-compilers for free, too.
5226 extendable: you can add your own classes to GiNaC, thus extending it on
5227 a very low level. Compare this to a traditional CAS that you can
5228 usually only extend on a high level by writing in the language defined
5229 by the parser. In particular, it turns out to be almost impossible to
5230 fix bugs in a traditional system.
5233 multiple interfaces: Though real GiNaC programs have to be written in
5234 some editor, then be compiled, linked and executed, there are more ways
5235 to work with the GiNaC engine. Many people want to play with
5236 expressions interactively, as in traditional CASs. Currently, two such
5237 windows into GiNaC have been implemented and many more are possible: the
5238 tiny @command{ginsh} that is part of the distribution exposes GiNaC's
5239 types to a command line and second, as a more consistent approach, an
5240 interactive interface to the Cint C++ interpreter has been put together
5241 (called GiNaC-cint) that allows an interactive scripting interface
5242 consistent with the C++ language. It is available from the usual GiNaC
5246 seamless integration: it is somewhere between difficult and impossible
5247 to call CAS functions from within a program written in C++ or any other
5248 programming language and vice versa. With GiNaC, your symbolic routines
5249 are part of your program. You can easily call third party libraries,
5250 e.g. for numerical evaluation or graphical interaction. All other
5251 approaches are much more cumbersome: they range from simply ignoring the
5252 problem (i.e. @emph{Maple}) to providing a method for `embedding' the
5253 system (i.e. @emph{Yacas}).
5256 efficiency: often large parts of a program do not need symbolic
5257 calculations at all. Why use large integers for loop variables or
5258 arbitrary precision arithmetics where @code{int} and @code{double} are
5259 sufficient? For pure symbolic applications, GiNaC is comparable in
5260 speed with other CAS.
5265 @node Disadvantages, Why C++?, Advantages, A Comparison With Other CAS
5266 @c node-name, next, previous, up
5267 @section Disadvantages
5269 Of course it also has some disadvantages:
5274 advanced features: GiNaC cannot compete with a program like
5275 @emph{Reduce} which exists for more than 30 years now or @emph{Maple}
5276 which grows since 1981 by the work of dozens of programmers, with
5277 respect to mathematical features. Integration, factorization,
5278 non-trivial simplifications, limits etc. are missing in GiNaC (and are
5279 not planned for the near future).
5282 portability: While the GiNaC library itself is designed to avoid any
5283 platform dependent features (it should compile on any ANSI compliant C++
5284 compiler), the currently used version of the CLN library (fast large
5285 integer and arbitrary precision arithmetics) can only by compiled
5286 without hassle on systems with the C++ compiler from the GNU Compiler
5287 Collection (GCC).@footnote{This is because CLN uses PROVIDE/REQUIRE like
5288 macros to let the compiler gather all static initializations, which
5289 works for GNU C++ only. Feel free to contact the authors in case you
5290 really believe that you need to use a different compiler. We have
5291 occasionally used other compilers and may be able to give you advice.}
5292 GiNaC uses recent language features like explicit constructors, mutable
5293 members, RTTI, @code{dynamic_cast}s and STL, so ANSI compliance is meant
5294 literally. Recent GCC versions starting at 2.95.3, although itself not
5295 yet ANSI compliant, support all needed features.
5300 @node Why C++?, Internal Structures, Disadvantages, A Comparison With Other CAS
5301 @c node-name, next, previous, up
5304 Why did we choose to implement GiNaC in C++ instead of Java or any other
5305 language? C++ is not perfect: type checking is not strict (casting is
5306 possible), separation between interface and implementation is not
5307 complete, object oriented design is not enforced. The main reason is
5308 the often scolded feature of operator overloading in C++. While it may
5309 be true that operating on classes with a @code{+} operator is rarely
5310 meaningful, it is perfectly suited for algebraic expressions. Writing
5311 @math{3x+5y} as @code{3*x+5*y} instead of
5312 @code{x.times(3).plus(y.times(5))} looks much more natural.
5313 Furthermore, the main developers are more familiar with C++ than with
5314 any other programming language.
5317 @node Internal Structures, Expressions are reference counted, Why C++? , Top
5318 @c node-name, next, previous, up
5319 @appendix Internal Structures
5322 * Expressions are reference counted::
5323 * Internal representation of products and sums::
5326 @node Expressions are reference counted, Internal representation of products and sums, Internal Structures, Internal Structures
5327 @c node-name, next, previous, up
5328 @appendixsection Expressions are reference counted
5330 @cindex reference counting
5331 @cindex copy-on-write
5332 @cindex garbage collection
5333 An expression is extremely light-weight since internally it works like a
5334 handle to the actual representation and really holds nothing more than a
5335 pointer to some other object. What this means in practice is that
5336 whenever you create two @code{ex} and set the second equal to the first
5337 no copying process is involved. Instead, the copying takes place as soon
5338 as you try to change the second. Consider the simple sequence of code:
5342 #include <ginac/ginac.h>
5343 using namespace std;
5344 using namespace GiNaC;
5348 symbol x("x"), y("y"), z("z");
5351 e1 = sin(x + 2*y) + 3*z + 41;
5352 e2 = e1; // e2 points to same object as e1
5353 cout << e2 << endl; // prints sin(x+2*y)+3*z+41
5354 e2 += 1; // e2 is copied into a new object
5355 cout << e2 << endl; // prints sin(x+2*y)+3*z+42
5359 The line @code{e2 = e1;} creates a second expression pointing to the
5360 object held already by @code{e1}. The time involved for this operation
5361 is therefore constant, no matter how large @code{e1} was. Actual
5362 copying, however, must take place in the line @code{e2 += 1;} because
5363 @code{e1} and @code{e2} are not handles for the same object any more.
5364 This concept is called @dfn{copy-on-write semantics}. It increases
5365 performance considerably whenever one object occurs multiple times and
5366 represents a simple garbage collection scheme because when an @code{ex}
5367 runs out of scope its destructor checks whether other expressions handle
5368 the object it points to too and deletes the object from memory if that
5369 turns out not to be the case. A slightly less trivial example of
5370 differentiation using the chain-rule should make clear how powerful this
5375 symbol x("x"), y("y");
5379 ex e3 = diff(sin(e2), x); // first derivative of sin(e2) by x
5380 cout << e1 << endl // prints x+3*y
5381 << e2 << endl // prints (x+3*y)^3
5382 << e3 << endl; // prints 3*(x+3*y)^2*cos((x+3*y)^3)
5386 Here, @code{e1} will actually be referenced three times while @code{e2}
5387 will be referenced two times. When the power of an expression is built,
5388 that expression needs not be copied. Likewise, since the derivative of
5389 a power of an expression can be easily expressed in terms of that
5390 expression, no copying of @code{e1} is involved when @code{e3} is
5391 constructed. So, when @code{e3} is constructed it will print as
5392 @code{3*(x+3*y)^2*cos((x+3*y)^3)} but the argument of @code{cos()} only
5393 holds a reference to @code{e2} and the factor in front is just
5396 As a user of GiNaC, you cannot see this mechanism of copy-on-write
5397 semantics. When you insert an expression into a second expression, the
5398 result behaves exactly as if the contents of the first expression were
5399 inserted. But it may be useful to remember that this is not what
5400 happens. Knowing this will enable you to write much more efficient
5401 code. If you still have an uncertain feeling with copy-on-write
5402 semantics, we recommend you have a look at the
5403 @uref{http://www.cerfnet.com/~mpcline/c++-faq-lite/, C++-FAQ lite} by
5404 Marshall Cline. Chapter 16 covers this issue and presents an
5405 implementation which is pretty close to the one in GiNaC.
5408 @node Internal representation of products and sums, Package Tools, Expressions are reference counted, Internal Structures
5409 @c node-name, next, previous, up
5410 @appendixsection Internal representation of products and sums
5412 @cindex representation
5415 @cindex @code{power}
5416 Although it should be completely transparent for the user of
5417 GiNaC a short discussion of this topic helps to understand the sources
5418 and also explain performance to a large degree. Consider the
5419 unexpanded symbolic expression
5421 $2d^3 \left( 4a + 5b - 3 \right)$
5424 @math{2*d^3*(4*a+5*b-3)}
5426 which could naively be represented by a tree of linear containers for
5427 addition and multiplication, one container for exponentiation with base
5428 and exponent and some atomic leaves of symbols and numbers in this
5433 @cindex pair-wise representation
5434 However, doing so results in a rather deeply nested tree which will
5435 quickly become inefficient to manipulate. We can improve on this by
5436 representing the sum as a sequence of terms, each one being a pair of a
5437 purely numeric multiplicative coefficient and its rest. In the same
5438 spirit we can store the multiplication as a sequence of terms, each
5439 having a numeric exponent and a possibly complicated base, the tree
5440 becomes much more flat:
5444 The number @code{3} above the symbol @code{d} shows that @code{mul}
5445 objects are treated similarly where the coefficients are interpreted as
5446 @emph{exponents} now. Addition of sums of terms or multiplication of
5447 products with numerical exponents can be coded to be very efficient with
5448 such a pair-wise representation. Internally, this handling is performed
5449 by most CAS in this way. It typically speeds up manipulations by an
5450 order of magnitude. The overall multiplicative factor @code{2} and the
5451 additive term @code{-3} look somewhat out of place in this
5452 representation, however, since they are still carrying a trivial
5453 exponent and multiplicative factor @code{1} respectively. Within GiNaC,
5454 this is avoided by adding a field that carries an overall numeric
5455 coefficient. This results in the realistic picture of internal
5458 $2d^3 \left( 4a + 5b - 3 \right)$:
5461 @math{2*d^3*(4*a+5*b-3)}:
5467 This also allows for a better handling of numeric radicals, since
5468 @code{sqrt(2)} can now be carried along calculations. Now it should be
5469 clear, why both classes @code{add} and @code{mul} are derived from the
5470 same abstract class: the data representation is the same, only the
5471 semantics differs. In the class hierarchy, methods for polynomial
5472 expansion and the like are reimplemented for @code{add} and @code{mul},
5473 but the data structure is inherited from @code{expairseq}.
5476 @node Package Tools, ginac-config, Internal representation of products and sums, Top
5477 @c node-name, next, previous, up
5478 @appendix Package Tools
5480 If you are creating a software package that uses the GiNaC library,
5481 setting the correct command line options for the compiler and linker
5482 can be difficult. GiNaC includes two tools to make this process easier.
5485 * ginac-config:: A shell script to detect compiler and linker flags.
5486 * AM_PATH_GINAC:: Macro for GNU automake.
5490 @node ginac-config, AM_PATH_GINAC, Package Tools, Package Tools
5491 @c node-name, next, previous, up
5492 @section @command{ginac-config}
5493 @cindex ginac-config
5495 @command{ginac-config} is a shell script that you can use to determine
5496 the compiler and linker command line options required to compile and
5497 link a program with the GiNaC library.
5499 @command{ginac-config} takes the following flags:
5503 Prints out the version of GiNaC installed.
5505 Prints '-I' flags pointing to the installed header files.
5507 Prints out the linker flags necessary to link a program against GiNaC.
5508 @item --prefix[=@var{PREFIX}]
5509 If @var{PREFIX} is specified, overrides the configured value of @env{$prefix}.
5510 (And of exec-prefix, unless @code{--exec-prefix} is also specified)
5511 Otherwise, prints out the configured value of @env{$prefix}.
5512 @item --exec-prefix[=@var{PREFIX}]
5513 If @var{PREFIX} is specified, overrides the configured value of @env{$exec_prefix}.
5514 Otherwise, prints out the configured value of @env{$exec_prefix}.
5517 Typically, @command{ginac-config} will be used within a configure
5518 script, as described below. It, however, can also be used directly from
5519 the command line using backquotes to compile a simple program. For
5523 c++ -o simple `ginac-config --cppflags` simple.cpp `ginac-config --libs`
5526 This command line might expand to (for example):
5529 cc -o simple -I/usr/local/include simple.cpp -L/usr/local/lib \
5530 -lginac -lcln -lstdc++
5533 Not only is the form using @command{ginac-config} easier to type, it will
5534 work on any system, no matter how GiNaC was configured.
5537 @node AM_PATH_GINAC, Configure script options, ginac-config, Package Tools
5538 @c node-name, next, previous, up
5539 @section @samp{AM_PATH_GINAC}
5540 @cindex AM_PATH_GINAC
5542 For packages configured using GNU automake, GiNaC also provides
5543 a macro to automate the process of checking for GiNaC.
5546 AM_PATH_GINAC([@var{MINIMUM-VERSION}, [@var{ACTION-IF-FOUND} [, @var{ACTION-IF-NOT-FOUND}]]])
5554 Determines the location of GiNaC using @command{ginac-config}, which is
5555 either found in the user's path, or from the environment variable
5556 @env{GINACLIB_CONFIG}.
5559 Tests the installed libraries to make sure that their version
5560 is later than @var{MINIMUM-VERSION}. (A default version will be used
5564 If the required version was found, sets the @env{GINACLIB_CPPFLAGS} variable
5565 to the output of @command{ginac-config --cppflags} and the @env{GINACLIB_LIBS}
5566 variable to the output of @command{ginac-config --libs}, and calls
5567 @samp{AC_SUBST()} for these variables so they can be used in generated
5568 makefiles, and then executes @var{ACTION-IF-FOUND}.
5571 If the required version was not found, sets @env{GINACLIB_CPPFLAGS} and
5572 @env{GINACLIB_LIBS} to empty strings, and executes @var{ACTION-IF-NOT-FOUND}.
5576 This macro is in file @file{ginac.m4} which is installed in
5577 @file{$datadir/aclocal}. Note that if automake was installed with a
5578 different @samp{--prefix} than GiNaC, you will either have to manually
5579 move @file{ginac.m4} to automake's @file{$datadir/aclocal}, or give
5580 aclocal the @samp{-I} option when running it.
5583 * Configure script options:: Configuring a package that uses AM_PATH_GINAC.
5584 * Example package:: Example of a package using AM_PATH_GINAC.
5588 @node Configure script options, Example package, AM_PATH_GINAC, AM_PATH_GINAC
5589 @c node-name, next, previous, up
5590 @subsection Configuring a package that uses @samp{AM_PATH_GINAC}
5592 Simply make sure that @command{ginac-config} is in your path, and run
5593 the configure script.
5600 The directory where the GiNaC libraries are installed needs
5601 to be found by your system's dynamic linker.
5603 This is generally done by
5606 editing @file{/etc/ld.so.conf} and running @command{ldconfig}
5612 setting the environment variable @env{LD_LIBRARY_PATH},
5615 or, as a last resort,
5618 giving a @samp{-R} or @samp{-rpath} flag (depending on your linker) when
5619 running configure, for instance:
5622 LDFLAGS=-R/home/cbauer/lib ./configure
5627 You can also specify a @command{ginac-config} not in your path by
5628 setting the @env{GINACLIB_CONFIG} environment variable to the
5629 name of the executable
5632 If you move the GiNaC package from its installed location,
5633 you will either need to modify @command{ginac-config} script
5634 manually to point to the new location or rebuild GiNaC.
5645 --with-ginac-prefix=@var{PREFIX}
5646 --with-ginac-exec-prefix=@var{PREFIX}
5649 are provided to override the prefix and exec-prefix that were stored
5650 in the @command{ginac-config} shell script by GiNaC's configure. You are
5651 generally better off configuring GiNaC with the right path to begin with.
5655 @node Example package, Bibliography, Configure script options, AM_PATH_GINAC
5656 @c node-name, next, previous, up
5657 @subsection Example of a package using @samp{AM_PATH_GINAC}
5659 The following shows how to build a simple package using automake
5660 and the @samp{AM_PATH_GINAC} macro. The program used here is @file{simple.cpp}:
5663 #include <ginac/ginac.h>
5667 GiNaC::symbol x("x");
5668 GiNaC::ex a = GiNaC::sin(x);
5669 std::cout << "Derivative of " << a
5670 << " is " << a.diff(x) << std::endl;
5675 You should first read the introductory portions of the automake
5676 Manual, if you are not already familiar with it.
5678 Two files are needed, @file{configure.in}, which is used to build the
5682 dnl Process this file with autoconf to produce a configure script.
5684 AM_INIT_AUTOMAKE(simple.cpp, 1.0.0)
5690 AM_PATH_GINAC(0.9.0, [
5691 LIBS="$LIBS $GINACLIB_LIBS"
5692 CPPFLAGS="$CPPFLAGS $GINACLIB_CPPFLAGS"
5693 ], AC_MSG_ERROR([need to have GiNaC installed]))
5698 The only command in this which is not standard for automake
5699 is the @samp{AM_PATH_GINAC} macro.
5701 That command does the following: If a GiNaC version greater or equal
5702 than 0.7.0 is found, then it adds @env{$GINACLIB_LIBS} to @env{$LIBS}
5703 and @env{$GINACLIB_CPPFLAGS} to @env{$CPPFLAGS}. Otherwise, it dies with
5704 the error message `need to have GiNaC installed'
5706 And the @file{Makefile.am}, which will be used to build the Makefile.
5709 ## Process this file with automake to produce Makefile.in
5710 bin_PROGRAMS = simple
5711 simple_SOURCES = simple.cpp
5714 This @file{Makefile.am}, says that we are building a single executable,
5715 from a single sourcefile @file{simple.cpp}. Since every program
5716 we are building uses GiNaC we simply added the GiNaC options
5717 to @env{$LIBS} and @env{$CPPFLAGS}, but in other circumstances, we might
5718 want to specify them on a per-program basis: for instance by
5722 simple_LDADD = $(GINACLIB_LIBS)
5723 INCLUDES = $(GINACLIB_CPPFLAGS)
5726 to the @file{Makefile.am}.
5728 To try this example out, create a new directory and add the three
5731 Now execute the following commands:
5734 $ automake --add-missing
5739 You now have a package that can be built in the normal fashion
5748 @node Bibliography, Concept Index, Example package, Top
5749 @c node-name, next, previous, up
5750 @appendix Bibliography
5755 @cite{ISO/IEC 14882:1998: Programming Languages: C++}
5758 @cite{CLN: A Class Library for Numbers}, @email{haible@@ilog.fr, Bruno Haible}
5761 @cite{The C++ Programming Language}, Bjarne Stroustrup, 3rd Edition, ISBN 0-201-88954-4, Addison Wesley
5764 @cite{C++ FAQs}, Marshall Cline, ISBN 0-201-58958-3, 1995, Addison Wesley
5767 @cite{Algorithms for Computer Algebra}, Keith O. Geddes, Stephen R. Czapor,
5768 and George Labahn, ISBN 0-7923-9259-0, 1992, Kluwer Academic Publishers, Norwell, Massachusetts
5771 @cite{Computer Algebra: Systems and Algorithms for Algebraic Computation},
5772 James H. Davenport, Yvon Siret, and Evelyne Tournier, ISBN 0-12-204230-1, 1988,
5773 Academic Press, London
5776 @cite{Computer Algebra Systems - A Practical Guide},
5777 Michael J. Wester (editor), ISBN 0-471-98353-5, 1999, Wiley, Chichester
5780 @cite{The Art of Computer Programming, Vol 2: Seminumerical Algorithms},
5781 Donald E. Knuth, ISBN 0-201-89684-2, 1998, Addison Wesley
5784 @cite{The Role of gamma5 in Dimensional Regularization}, Dirk Kreimer, hep-ph/9401354
5789 @node Concept Index, , Bibliography, Top
5790 @c node-name, next, previous, up
5791 @unnumbered Concept Index