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()}
1217 @cindex @code{remove_all()}
1219 The GiNaC class @code{lst} serves for holding a @dfn{list} of arbitrary
1220 expressions. These are sometimes used to supply a variable number of
1221 arguments of the same type to GiNaC methods such as @code{subs()} and
1222 @code{to_rational()}, so you should have a basic understanding about them.
1224 Lists of up to 16 expressions can be directly constructed from single
1229 symbol x("x"), y("y");
1230 lst l(x, 2, y, x+y);
1231 // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y'
1235 Use the @code{nops()} method to determine the size (number of expressions) of
1236 a list and the @code{op()} method to access individual elements:
1240 cout << l.nops() << endl; // prints '4'
1241 cout << l.op(2) << " " << l.op(0) << endl; // prints 'y x'
1245 You can append or prepend an expression to a list with the @code{append()}
1246 and @code{prepend()} methods:
1250 l.append(4*x); // l is now @{x, 2, y, x+y, 4*x@}
1251 l.prepend(0); // l is now @{0, x, 2, y, x+y, 4*x@}
1255 You can remove the first or last element of a list with
1256 @code{remove_first()} and @code{remove_last()}:
1260 l.remove_first(); // l is now @{x, 2, y, x+y, 4*x@}
1261 l.remove_last(); // l is now @{x, 2, y, x+y@}
1264 Finally, you can remove all the elements of a list with
1265 @code{remove_all()}:
1269 l.remove_all(); // l is now empty
1274 @node Mathematical functions, Relations, Lists, Basic Concepts
1275 @c node-name, next, previous, up
1276 @section Mathematical functions
1277 @cindex @code{function} (class)
1278 @cindex trigonometric function
1279 @cindex hyperbolic function
1281 There are quite a number of useful functions hard-wired into GiNaC. For
1282 instance, all trigonometric and hyperbolic functions are implemented
1283 (@xref{Built-in Functions}, for a complete list).
1285 These functions (better called @emph{pseudofunctions}) are all objects
1286 of class @code{function}. They accept one or more expressions as
1287 arguments and return one expression. If the arguments are not
1288 numerical, the evaluation of the function may be halted, as it does in
1289 the next example, showing how a function returns itself twice and
1290 finally an expression that may be really useful:
1292 @cindex Gamma function
1293 @cindex @code{subs()}
1296 symbol x("x"), y("y");
1298 cout << tgamma(foo) << endl;
1299 // -> tgamma(x+(1/2)*y)
1300 ex bar = foo.subs(y==1);
1301 cout << tgamma(bar) << endl;
1303 ex foobar = bar.subs(x==7);
1304 cout << tgamma(foobar) << endl;
1305 // -> (135135/128)*Pi^(1/2)
1309 Besides evaluation most of these functions allow differentiation, series
1310 expansion and so on. Read the next chapter in order to learn more about
1313 It must be noted that these pseudofunctions are created by inline
1314 functions, where the argument list is templated. This means that
1315 whenever you call @code{GiNaC::sin(1)} it is equivalent to
1316 @code{sin(ex(1))} and will therefore not result in a floating point
1317 number. Unless of course the function prototype is explicitly
1318 overridden -- which is the case for arguments of type @code{numeric}
1319 (not wrapped inside an @code{ex}). Hence, in order to obtain a floating
1320 point number of class @code{numeric} you should call
1321 @code{sin(numeric(1))}. This is almost the same as calling
1322 @code{sin(1).evalf()} except that the latter will return a numeric
1323 wrapped inside an @code{ex}.
1326 @node Relations, Matrices, Mathematical functions, Basic Concepts
1327 @c node-name, next, previous, up
1329 @cindex @code{relational} (class)
1331 Sometimes, a relation holding between two expressions must be stored
1332 somehow. The class @code{relational} is a convenient container for such
1333 purposes. A relation is by definition a container for two @code{ex} and
1334 a relation between them that signals equality, inequality and so on.
1335 They are created by simply using the C++ operators @code{==}, @code{!=},
1336 @code{<}, @code{<=}, @code{>} and @code{>=} between two expressions.
1338 @xref{Mathematical functions}, for examples where various applications
1339 of the @code{.subs()} method show how objects of class relational are
1340 used as arguments. There they provide an intuitive syntax for
1341 substitutions. They are also used as arguments to the @code{ex::series}
1342 method, where the left hand side of the relation specifies the variable
1343 to expand in and the right hand side the expansion point. They can also
1344 be used for creating systems of equations that are to be solved for
1345 unknown variables. But the most common usage of objects of this class
1346 is rather inconspicuous in statements of the form @code{if
1347 (expand(pow(a+b,2))==a*a+2*a*b+b*b) @{...@}}. Here, an implicit
1348 conversion from @code{relational} to @code{bool} takes place. Note,
1349 however, that @code{==} here does not perform any simplifications, hence
1350 @code{expand()} must be called explicitly.
1353 @node Matrices, Indexed objects, Relations, Basic Concepts
1354 @c node-name, next, previous, up
1356 @cindex @code{matrix} (class)
1358 A @dfn{matrix} is a two-dimensional array of expressions. The elements of a
1359 matrix with @math{m} rows and @math{n} columns are accessed with two
1360 @code{unsigned} indices, the first one in the range 0@dots{}@math{m-1}, the
1361 second one in the range 0@dots{}@math{n-1}.
1363 There are a couple of ways to construct matrices, with or without preset
1367 matrix::matrix(unsigned r, unsigned c);
1368 matrix::matrix(unsigned r, unsigned c, const lst & l);
1369 ex lst_to_matrix(const lst & l);
1370 ex diag_matrix(const lst & l);
1373 The first two functions are @code{matrix} constructors which create a matrix
1374 with @samp{r} rows and @samp{c} columns. The matrix elements can be
1375 initialized from a (flat) list of expressions @samp{l}. Otherwise they are
1376 all set to zero. The @code{lst_to_matrix()} function constructs a matrix
1377 from a list of lists, each list representing a matrix row. Finally,
1378 @code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
1379 elements. Note that the last two functions return expressions, not matrix
1382 Matrix elements can be accessed and set using the parenthesis (function call)
1386 const ex & matrix::operator()(unsigned r, unsigned c) const;
1387 ex & matrix::operator()(unsigned r, unsigned c);
1390 It is also possible to access the matrix elements in a linear fashion with
1391 the @code{op()} method. But C++-style subscripting with square brackets
1392 @samp{[]} is not available.
1394 Here are a couple of examples that all construct the same 2x2 diagonal
1399 symbol a("a"), b("b");
1407 e = matrix(2, 2, lst(a, 0, 0, b));
1409 e = lst_to_matrix(lst(lst(a, 0), lst(0, b)));
1411 e = diag_matrix(lst(a, b));
1418 @cindex @code{transpose()}
1419 @cindex @code{inverse()}
1420 There are three ways to do arithmetic with matrices. The first (and most
1421 efficient one) is to use the methods provided by the @code{matrix} class:
1424 matrix matrix::add(const matrix & other) const;
1425 matrix matrix::sub(const matrix & other) const;
1426 matrix matrix::mul(const matrix & other) const;
1427 matrix matrix::mul_scalar(const ex & other) const;
1428 matrix matrix::pow(const ex & expn) const;
1429 matrix matrix::transpose(void) const;
1430 matrix matrix::inverse(void) const;
1433 All of these methods return the result as a new matrix object. Here is an
1434 example that calculates @math{A*B-2*C} for three matrices @math{A}, @math{B}
1439 matrix A(2, 2, lst(1, 2, 3, 4));
1440 matrix B(2, 2, lst(-1, 0, 2, 1));
1441 matrix C(2, 2, lst(8, 4, 2, 1));
1443 matrix result = A.mul(B).sub(C.mul_scalar(2));
1444 cout << result << endl;
1445 // -> [[-13,-6],[1,2]]
1450 @cindex @code{evalm()}
1451 The second (and probably the most natural) way is to construct an expression
1452 containing matrices with the usual arithmetic operators and @code{pow()}.
1453 For efficiency reasons, expressions with sums, products and powers of
1454 matrices are not automatically evaluated in GiNaC. You have to call the
1458 ex ex::evalm() const;
1461 to obtain the result:
1468 // -> [[1,2],[3,4]]*[[-1,0],[2,1]]-2*[[8,4],[2,1]]
1469 cout << e.evalm() << endl;
1470 // -> [[-13,-6],[1,2]]
1475 The non-commutativity of the product @code{A*B} in this example is
1476 automatically recognized by GiNaC. There is no need to use a special
1477 operator here. @xref{Non-commutative objects}, for more information about
1478 dealing with non-commutative expressions.
1480 Finally, you can work with indexed matrices and call @code{simplify_indexed()}
1481 to perform the arithmetic:
1486 idx i(symbol("i"), 2), j(symbol("j"), 2), k(symbol("k"), 2);
1487 e = indexed(A, i, k) * indexed(B, k, j) - 2 * indexed(C, i, j);
1489 // -> -2*[[8,4],[2,1]].i.j+[[-1,0],[2,1]].k.j*[[1,2],[3,4]].i.k
1490 cout << e.simplify_indexed() << endl;
1491 // -> [[-13,-6],[1,2]].i.j
1495 Using indices is most useful when working with rectangular matrices and
1496 one-dimensional vectors because you don't have to worry about having to
1497 transpose matrices before multiplying them. @xref{Indexed objects}, for
1498 more information about using matrices with indices, and about indices in
1501 The @code{matrix} class provides a couple of additional methods for
1502 computing determinants, traces, and characteristic polynomials:
1505 ex matrix::determinant(unsigned algo = determinant_algo::automatic) const;
1506 ex matrix::trace(void) const;
1507 ex matrix::charpoly(const symbol & lambda) const;
1510 The @samp{algo} argument of @code{determinant()} allows to select between
1511 different algorithms for calculating the determinant. The possible values
1512 are defined in the @file{flags.h} header file. By default, GiNaC uses a
1513 heuristic to automatically select an algorithm that is likely to give the
1514 result most quickly.
1517 @node Indexed objects, Non-commutative objects, Matrices, Basic Concepts
1518 @c node-name, next, previous, up
1519 @section Indexed objects
1521 GiNaC allows you to handle expressions containing general indexed objects in
1522 arbitrary spaces. It is also able to canonicalize and simplify such
1523 expressions and perform symbolic dummy index summations. There are a number
1524 of predefined indexed objects provided, like delta and metric tensors.
1526 There are few restrictions placed on indexed objects and their indices and
1527 it is easy to construct nonsense expressions, but our intention is to
1528 provide a general framework that allows you to implement algorithms with
1529 indexed quantities, getting in the way as little as possible.
1531 @cindex @code{idx} (class)
1532 @cindex @code{indexed} (class)
1533 @subsection Indexed quantities and their indices
1535 Indexed expressions in GiNaC are constructed of two special types of objects,
1536 @dfn{index objects} and @dfn{indexed objects}.
1540 @cindex contravariant
1543 @item Index objects are of class @code{idx} or a subclass. Every index has
1544 a @dfn{value} and a @dfn{dimension} (which is the dimension of the space
1545 the index lives in) which can both be arbitrary expressions but are usually
1546 a number or a simple symbol. In addition, indices of class @code{varidx} have
1547 a @dfn{variance} (they can be co- or contravariant), and indices of class
1548 @code{spinidx} have a variance and can be @dfn{dotted} or @dfn{undotted}.
1550 @item Indexed objects are of class @code{indexed} or a subclass. They
1551 contain a @dfn{base expression} (which is the expression being indexed), and
1552 one or more indices.
1556 @strong{Note:} when printing expressions, covariant indices and indices
1557 without variance are denoted @samp{.i} while contravariant indices are
1558 denoted @samp{~i}. Dotted indices have a @samp{*} in front of the index
1559 value. In the following, we are going to use that notation in the text so
1560 instead of @math{A^i_jk} we will write @samp{A~i.j.k}. Index dimensions are
1561 not visible in the output.
1563 A simple example shall illustrate the concepts:
1567 #include <ginac/ginac.h>
1568 using namespace std;
1569 using namespace GiNaC;
1573 symbol i_sym("i"), j_sym("j");
1574 idx i(i_sym, 3), j(j_sym, 3);
1577 cout << indexed(A, i, j) << endl;
1582 The @code{idx} constructor takes two arguments, the index value and the
1583 index dimension. First we define two index objects, @code{i} and @code{j},
1584 both with the numeric dimension 3. The value of the index @code{i} is the
1585 symbol @code{i_sym} (which prints as @samp{i}) and the value of the index
1586 @code{j} is the symbol @code{j_sym} (which prints as @samp{j}). Next we
1587 construct an expression containing one indexed object, @samp{A.i.j}. It has
1588 the symbol @code{A} as its base expression and the two indices @code{i} and
1591 Note the difference between the indices @code{i} and @code{j} which are of
1592 class @code{idx}, and the index values which are the symbols @code{i_sym}
1593 and @code{j_sym}. The indices of indexed objects cannot directly be symbols
1594 or numbers but must be index objects. For example, the following is not
1595 correct and will raise an exception:
1598 symbol i("i"), j("j");
1599 e = indexed(A, i, j); // ERROR: indices must be of type idx
1602 You can have multiple indexed objects in an expression, index values can
1603 be numeric, and index dimensions symbolic:
1607 symbol B("B"), dim("dim");
1608 cout << 4 * indexed(A, i)
1609 + indexed(B, idx(j_sym, 4), idx(2, 3), idx(i_sym, dim)) << endl;
1614 @code{B} has a 4-dimensional symbolic index @samp{k}, a 3-dimensional numeric
1615 index of value 2, and a symbolic index @samp{i} with the symbolic dimension
1616 @samp{dim}. Note that GiNaC doesn't automatically notify you that the free
1617 indices of @samp{A} and @samp{B} in the sum don't match (you have to call
1618 @code{simplify_indexed()} for that, see below).
1620 In fact, base expressions, index values and index dimensions can be
1621 arbitrary expressions:
1625 cout << indexed(A+B, idx(2*i_sym+1, dim/2)) << endl;
1630 It's also possible to construct nonsense like @samp{Pi.sin(x)}. You will not
1631 get an error message from this but you will probably not be able to do
1632 anything useful with it.
1634 @cindex @code{get_value()}
1635 @cindex @code{get_dimension()}
1639 ex idx::get_value(void);
1640 ex idx::get_dimension(void);
1643 return the value and dimension of an @code{idx} object. If you have an index
1644 in an expression, such as returned by calling @code{.op()} on an indexed
1645 object, you can get a reference to the @code{idx} object with the function
1646 @code{ex_to<idx>()} on the expression.
1648 There are also the methods
1651 bool idx::is_numeric(void);
1652 bool idx::is_symbolic(void);
1653 bool idx::is_dim_numeric(void);
1654 bool idx::is_dim_symbolic(void);
1657 for checking whether the value and dimension are numeric or symbolic
1658 (non-numeric). Using the @code{info()} method of an index (see @ref{Information
1659 About Expressions}) returns information about the index value.
1661 @cindex @code{varidx} (class)
1662 If you need co- and contravariant indices, use the @code{varidx} class:
1666 symbol mu_sym("mu"), nu_sym("nu");
1667 varidx mu(mu_sym, 4), nu(nu_sym, 4); // default is contravariant ~mu, ~nu
1668 varidx mu_co(mu_sym, 4, true); // covariant index .mu
1670 cout << indexed(A, mu, nu) << endl;
1672 cout << indexed(A, mu_co, nu) << endl;
1674 cout << indexed(A, mu.toggle_variance(), nu) << endl;
1679 A @code{varidx} is an @code{idx} with an additional flag that marks it as
1680 co- or contravariant. The default is a contravariant (upper) index, but
1681 this can be overridden by supplying a third argument to the @code{varidx}
1682 constructor. The two methods
1685 bool varidx::is_covariant(void);
1686 bool varidx::is_contravariant(void);
1689 allow you to check the variance of a @code{varidx} object (use @code{ex_to<varidx>()}
1690 to get the object reference from an expression). There's also the very useful
1694 ex varidx::toggle_variance(void);
1697 which makes a new index with the same value and dimension but the opposite
1698 variance. By using it you only have to define the index once.
1700 @cindex @code{spinidx} (class)
1701 The @code{spinidx} class provides dotted and undotted variant indices, as
1702 used in the Weyl-van-der-Waerden spinor formalism:
1706 symbol K("K"), C_sym("C"), D_sym("D");
1707 spinidx C(C_sym, 2), D(D_sym); // default is 2-dimensional,
1708 // contravariant, undotted
1709 spinidx C_co(C_sym, 2, true); // covariant index
1710 spinidx D_dot(D_sym, 2, false, true); // contravariant, dotted
1711 spinidx D_co_dot(D_sym, 2, true, true); // covariant, dotted
1713 cout << indexed(K, C, D) << endl;
1715 cout << indexed(K, C_co, D_dot) << endl;
1717 cout << indexed(K, D_co_dot, D) << endl;
1722 A @code{spinidx} is a @code{varidx} with an additional flag that marks it as
1723 dotted or undotted. The default is undotted but this can be overridden by
1724 supplying a fourth argument to the @code{spinidx} constructor. The two
1728 bool spinidx::is_dotted(void);
1729 bool spinidx::is_undotted(void);
1732 allow you to check whether or not a @code{spinidx} object is dotted (use
1733 @code{ex_to<spinidx>()} to get the object reference from an expression).
1734 Finally, the two methods
1737 ex spinidx::toggle_dot(void);
1738 ex spinidx::toggle_variance_dot(void);
1741 create a new index with the same value and dimension but opposite dottedness
1742 and the same or opposite variance.
1744 @subsection Substituting indices
1746 @cindex @code{subs()}
1747 Sometimes you will want to substitute one symbolic index with another
1748 symbolic or numeric index, for example when calculating one specific element
1749 of a tensor expression. This is done with the @code{.subs()} method, as it
1750 is done for symbols (see @ref{Substituting Expressions}).
1752 You have two possibilities here. You can either substitute the whole index
1753 by another index or expression:
1757 ex e = indexed(A, mu_co);
1758 cout << e << " becomes " << e.subs(mu_co == nu) << endl;
1759 // -> A.mu becomes A~nu
1760 cout << e << " becomes " << e.subs(mu_co == varidx(0, 4)) << endl;
1761 // -> A.mu becomes A~0
1762 cout << e << " becomes " << e.subs(mu_co == 0) << endl;
1763 // -> A.mu becomes A.0
1767 The third example shows that trying to replace an index with something that
1768 is not an index will substitute the index value instead.
1770 Alternatively, you can substitute the @emph{symbol} of a symbolic index by
1775 ex e = indexed(A, mu_co);
1776 cout << e << " becomes " << e.subs(mu_sym == nu_sym) << endl;
1777 // -> A.mu becomes A.nu
1778 cout << e << " becomes " << e.subs(mu_sym == 0) << endl;
1779 // -> A.mu becomes A.0
1783 As you see, with the second method only the value of the index will get
1784 substituted. Its other properties, including its dimension, remain unchanged.
1785 If you want to change the dimension of an index you have to substitute the
1786 whole index by another one with the new dimension.
1788 Finally, substituting the base expression of an indexed object works as
1793 ex e = indexed(A, mu_co);
1794 cout << e << " becomes " << e.subs(A == A+B) << endl;
1795 // -> A.mu becomes (B+A).mu
1799 @subsection Symmetries
1800 @cindex @code{symmetry} (class)
1801 @cindex @code{sy_none()}
1802 @cindex @code{sy_symm()}
1803 @cindex @code{sy_anti()}
1804 @cindex @code{sy_cycl()}
1806 Indexed objects can have certain symmetry properties with respect to their
1807 indices. Symmetries are specified as a tree of objects of class @code{symmetry}
1808 that is constructed with the helper functions
1811 symmetry sy_none(...);
1812 symmetry sy_symm(...);
1813 symmetry sy_anti(...);
1814 symmetry sy_cycl(...);
1817 @code{sy_none()} stands for no symmetry, @code{sy_symm()} and @code{sy_anti()}
1818 specify fully symmetric or antisymmetric, respectively, and @code{sy_cycl()}
1819 represents a cyclic symmetry. Each of these functions accepts up to four
1820 arguments which can be either symmetry objects themselves or unsigned integer
1821 numbers that represent an index position (counting from 0). A symmetry
1822 specification that consists of only a single @code{sy_symm()}, @code{sy_anti()}
1823 or @code{sy_cycl()} with no arguments specifies the respective symmetry for
1826 Here are some examples of symmetry definitions:
1831 e = indexed(A, i, j);
1832 e = indexed(A, sy_none(), i, j); // equivalent
1833 e = indexed(A, sy_none(0, 1), i, j); // equivalent
1835 // Symmetric in all three indices:
1836 e = indexed(A, sy_symm(), i, j, k);
1837 e = indexed(A, sy_symm(0, 1, 2), i, j, k); // equivalent
1838 e = indexed(A, sy_symm(2, 0, 1), i, j, k); // same symmetry, but yields a
1839 // different canonical order
1841 // Symmetric in the first two indices only:
1842 e = indexed(A, sy_symm(0, 1), i, j, k);
1843 e = indexed(A, sy_none(sy_symm(0, 1), 2), i, j, k); // equivalent
1845 // Antisymmetric in the first and last index only (index ranges need not
1847 e = indexed(A, sy_anti(0, 2), i, j, k);
1848 e = indexed(A, sy_none(sy_anti(0, 2), 1), i, j, k); // equivalent
1850 // An example of a mixed symmetry: antisymmetric in the first two and
1851 // last two indices, symmetric when swapping the first and last index
1852 // pairs (like the Riemann curvature tensor):
1853 e = indexed(A, sy_symm(sy_anti(0, 1), sy_anti(2, 3)), i, j, k, l);
1855 // Cyclic symmetry in all three indices:
1856 e = indexed(A, sy_cycl(), i, j, k);
1857 e = indexed(A, sy_cycl(0, 1, 2), i, j, k); // equivalent
1859 // The following examples are invalid constructions that will throw
1860 // an exception at run time.
1862 // An index may not appear multiple times:
1863 e = indexed(A, sy_symm(0, 0, 1), i, j, k); // ERROR
1864 e = indexed(A, sy_none(sy_symm(0, 1), sy_anti(0, 2)), i, j, k); // ERROR
1866 // Every child of sy_symm(), sy_anti() and sy_cycl() must refer to the
1867 // same number of indices:
1868 e = indexed(A, sy_symm(sy_anti(0, 1), 2), i, j, k); // ERROR
1870 // And of course, you cannot specify indices which are not there:
1871 e = indexed(A, sy_symm(0, 1, 2, 3), i, j, k); // ERROR
1875 If you need to specify more than four indices, you have to use the
1876 @code{.add()} method of the @code{symmetry} class. For example, to specify
1877 full symmetry in the first six indices you would write
1878 @code{sy_symm(0, 1, 2, 3).add(4).add(5)}.
1880 If an indexed object has a symmetry, GiNaC will automatically bring the
1881 indices into a canonical order which allows for some immediate simplifications:
1885 cout << indexed(A, sy_symm(), i, j)
1886 + indexed(A, sy_symm(), j, i) << endl;
1888 cout << indexed(B, sy_anti(), i, j)
1889 + indexed(B, sy_anti(), j, i) << endl;
1891 cout << indexed(B, sy_anti(), i, j, k)
1892 - indexed(B, sy_anti(), j, k, i) << endl;
1897 @cindex @code{get_free_indices()}
1899 @subsection Dummy indices
1901 GiNaC treats certain symbolic index pairs as @dfn{dummy indices} meaning
1902 that a summation over the index range is implied. Symbolic indices which are
1903 not dummy indices are called @dfn{free indices}. Numeric indices are neither
1904 dummy nor free indices.
1906 To be recognized as a dummy index pair, the two indices must be of the same
1907 class and their value must be the same single symbol (an index like
1908 @samp{2*n+1} is never a dummy index). If the indices are of class
1909 @code{varidx} they must also be of opposite variance; if they are of class
1910 @code{spinidx} they must be both dotted or both undotted.
1912 The method @code{.get_free_indices()} returns a vector containing the free
1913 indices of an expression. It also checks that the free indices of the terms
1914 of a sum are consistent:
1918 symbol A("A"), B("B"), C("C");
1920 symbol i_sym("i"), j_sym("j"), k_sym("k"), l_sym("l");
1921 idx i(i_sym, 3), j(j_sym, 3), k(k_sym, 3), l(l_sym, 3);
1923 ex e = indexed(A, i, j) * indexed(B, j, k) + indexed(C, k, l, i, l);
1924 cout << exprseq(e.get_free_indices()) << endl;
1926 // 'j' and 'l' are dummy indices
1928 symbol mu_sym("mu"), nu_sym("nu"), rho_sym("rho"), sigma_sym("sigma");
1929 varidx mu(mu_sym, 4), nu(nu_sym, 4), rho(rho_sym, 4), sigma(sigma_sym, 4);
1931 e = indexed(A, mu, nu) * indexed(B, nu.toggle_variance(), rho)
1932 + indexed(C, mu, sigma, rho, sigma.toggle_variance());
1933 cout << exprseq(e.get_free_indices()) << endl;
1935 // 'nu' is a dummy index, but 'sigma' is not
1937 e = indexed(A, mu, mu);
1938 cout << exprseq(e.get_free_indices()) << endl;
1940 // 'mu' is not a dummy index because it appears twice with the same
1943 e = indexed(A, mu, nu) + 42;
1944 cout << exprseq(e.get_free_indices()) << endl; // ERROR
1945 // this will throw an exception:
1946 // "add::get_free_indices: inconsistent indices in sum"
1950 @cindex @code{simplify_indexed()}
1951 @subsection Simplifying indexed expressions
1953 In addition to the few automatic simplifications that GiNaC performs on
1954 indexed expressions (such as re-ordering the indices of symmetric tensors
1955 and calculating traces and convolutions of matrices and predefined tensors)
1959 ex ex::simplify_indexed(void);
1960 ex ex::simplify_indexed(const scalar_products & sp);
1963 that performs some more expensive operations:
1966 @item it checks the consistency of free indices in sums in the same way
1967 @code{get_free_indices()} does
1968 @item it tries to give dummy indices that appear in different terms of a sum
1969 the same name to allow simplifications like @math{a_i*b_i-a_j*b_j=0}
1970 @item it (symbolically) calculates all possible dummy index summations/contractions
1971 with the predefined tensors (this will be explained in more detail in the
1973 @item it detects contractions that vanish for symmetry reasons, for example
1974 the contraction of a symmetric and a totally antisymmetric tensor
1975 @item as a special case of dummy index summation, it can replace scalar products
1976 of two tensors with a user-defined value
1979 The last point is done with the help of the @code{scalar_products} class
1980 which is used to store scalar products with known values (this is not an
1981 arithmetic class, you just pass it to @code{simplify_indexed()}):
1985 symbol A("A"), B("B"), C("C"), i_sym("i");
1989 sp.add(A, B, 0); // A and B are orthogonal
1990 sp.add(A, C, 0); // A and C are orthogonal
1991 sp.add(A, A, 4); // A^2 = 4 (A has length 2)
1993 e = indexed(A + B, i) * indexed(A + C, i);
1995 // -> (B+A).i*(A+C).i
1997 cout << e.expand(expand_options::expand_indexed).simplify_indexed(sp)
2003 The @code{scalar_products} object @code{sp} acts as a storage for the
2004 scalar products added to it with the @code{.add()} method. This method
2005 takes three arguments: the two expressions of which the scalar product is
2006 taken, and the expression to replace it with. After @code{sp.add(A, B, 0)},
2007 @code{simplify_indexed()} will replace all scalar products of indexed
2008 objects that have the symbols @code{A} and @code{B} as base expressions
2009 with the single value 0. The number, type and dimension of the indices
2010 don't matter; @samp{A~mu~nu*B.mu.nu} would also be replaced by 0.
2012 @cindex @code{expand()}
2013 The example above also illustrates a feature of the @code{expand()} method:
2014 if passed the @code{expand_indexed} option it will distribute indices
2015 over sums, so @samp{(A+B).i} becomes @samp{A.i+B.i}.
2017 @cindex @code{tensor} (class)
2018 @subsection Predefined tensors
2020 Some frequently used special tensors such as the delta, epsilon and metric
2021 tensors are predefined in GiNaC. They have special properties when
2022 contracted with other tensor expressions and some of them have constant
2023 matrix representations (they will evaluate to a number when numeric
2024 indices are specified).
2026 @cindex @code{delta_tensor()}
2027 @subsubsection Delta tensor
2029 The delta tensor takes two indices, is symmetric and has the matrix
2030 representation @code{diag(1, 1, 1, ...)}. It is constructed by the function
2031 @code{delta_tensor()}:
2035 symbol A("A"), B("B");
2037 idx i(symbol("i"), 3), j(symbol("j"), 3),
2038 k(symbol("k"), 3), l(symbol("l"), 3);
2040 ex e = indexed(A, i, j) * indexed(B, k, l)
2041 * delta_tensor(i, k) * delta_tensor(j, l) << endl;
2042 cout << e.simplify_indexed() << endl;
2045 cout << delta_tensor(i, i) << endl;
2050 @cindex @code{metric_tensor()}
2051 @subsubsection General metric tensor
2053 The function @code{metric_tensor()} creates a general symmetric metric
2054 tensor with two indices that can be used to raise/lower tensor indices. The
2055 metric tensor is denoted as @samp{g} in the output and if its indices are of
2056 mixed variance it is automatically replaced by a delta tensor:
2062 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
2064 ex e = metric_tensor(mu, nu) * indexed(A, nu.toggle_variance(), rho);
2065 cout << e.simplify_indexed() << endl;
2068 e = delta_tensor(mu, nu.toggle_variance()) * metric_tensor(nu, rho);
2069 cout << e.simplify_indexed() << endl;
2072 e = metric_tensor(mu.toggle_variance(), nu.toggle_variance())
2073 * metric_tensor(nu, rho);
2074 cout << e.simplify_indexed() << endl;
2077 e = metric_tensor(nu.toggle_variance(), rho.toggle_variance())
2078 * metric_tensor(mu, nu) * (delta_tensor(mu.toggle_variance(), rho)
2079 + indexed(A, mu.toggle_variance(), rho));
2080 cout << e.simplify_indexed() << endl;
2085 @cindex @code{lorentz_g()}
2086 @subsubsection Minkowski metric tensor
2088 The Minkowski metric tensor is a special metric tensor with a constant
2089 matrix representation which is either @code{diag(1, -1, -1, ...)} (negative
2090 signature, the default) or @code{diag(-1, 1, 1, ...)} (positive signature).
2091 It is created with the function @code{lorentz_g()} (although it is output as
2096 varidx mu(symbol("mu"), 4);
2098 e = delta_tensor(varidx(0, 4), mu.toggle_variance())
2099 * lorentz_g(mu, varidx(0, 4)); // negative signature
2100 cout << e.simplify_indexed() << endl;
2103 e = delta_tensor(varidx(0, 4), mu.toggle_variance())
2104 * lorentz_g(mu, varidx(0, 4), true); // positive signature
2105 cout << e.simplify_indexed() << endl;
2110 @cindex @code{spinor_metric()}
2111 @subsubsection Spinor metric tensor
2113 The function @code{spinor_metric()} creates an antisymmetric tensor with
2114 two indices that is used to raise/lower indices of 2-component spinors.
2115 It is output as @samp{eps}:
2121 spinidx A(symbol("A")), B(symbol("B")), C(symbol("C"));
2122 ex A_co = A.toggle_variance(), B_co = B.toggle_variance();
2124 e = spinor_metric(A, B) * indexed(psi, B_co);
2125 cout << e.simplify_indexed() << endl;
2128 e = spinor_metric(A, B) * indexed(psi, A_co);
2129 cout << e.simplify_indexed() << endl;
2132 e = spinor_metric(A_co, B_co) * indexed(psi, B);
2133 cout << e.simplify_indexed() << endl;
2136 e = spinor_metric(A_co, B_co) * indexed(psi, A);
2137 cout << e.simplify_indexed() << endl;
2140 e = spinor_metric(A_co, B_co) * spinor_metric(A, B);
2141 cout << e.simplify_indexed() << endl;
2144 e = spinor_metric(A_co, B_co) * spinor_metric(B, C);
2145 cout << e.simplify_indexed() << endl;
2150 The matrix representation of the spinor metric is @code{[[0, 1], [-1, 0]]}.
2152 @cindex @code{epsilon_tensor()}
2153 @cindex @code{lorentz_eps()}
2154 @subsubsection Epsilon tensor
2156 The epsilon tensor is totally antisymmetric, its number of indices is equal
2157 to the dimension of the index space (the indices must all be of the same
2158 numeric dimension), and @samp{eps.1.2.3...} (resp. @samp{eps~0~1~2...}) is
2159 defined to be 1. Its behavior with indices that have a variance also
2160 depends on the signature of the metric. Epsilon tensors are output as
2163 There are three functions defined to create epsilon tensors in 2, 3 and 4
2167 ex epsilon_tensor(const ex & i1, const ex & i2);
2168 ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
2169 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig = false);
2172 The first two functions create an epsilon tensor in 2 or 3 Euclidean
2173 dimensions, the last function creates an epsilon tensor in a 4-dimensional
2174 Minkowski space (the last @code{bool} argument specifies whether the metric
2175 has negative or positive signature, as in the case of the Minkowski metric
2180 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4),
2181 sig(symbol("sig"), 4), lam(symbol("lam"), 4), bet(symbol("bet"), 4);
2182 e = lorentz_eps(mu, nu, rho, sig) *
2183 lorentz_eps(mu.toggle_variance(), nu.toggle_variance(), lam, bet);
2184 cout << simplify_indexed(e) << endl;
2185 // -> 2*eta~bet~rho*eta~sig~lam-2*eta~sig~bet*eta~rho~lam
2187 idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
2188 symbol A("A"), B("B");
2189 e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(B, k);
2190 cout << simplify_indexed(e) << endl;
2191 // -> -B.k*A.j*eps.i.k.j
2192 e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(A, k);
2193 cout << simplify_indexed(e) << endl;
2198 @subsection Linear algebra
2200 The @code{matrix} class can be used with indices to do some simple linear
2201 algebra (linear combinations and products of vectors and matrices, traces
2202 and scalar products):
2206 idx i(symbol("i"), 2), j(symbol("j"), 2);
2207 symbol x("x"), y("y");
2209 // A is a 2x2 matrix, X is a 2x1 vector
2210 matrix A(2, 2, lst(1, 2, 3, 4)), X(2, 1, lst(x, y));
2212 cout << indexed(A, i, i) << endl;
2215 ex e = indexed(A, i, j) * indexed(X, j);
2216 cout << e.simplify_indexed() << endl;
2217 // -> [[2*y+x],[4*y+3*x]].i
2219 e = indexed(A, i, j) * indexed(X, i) + indexed(X, j) * 2;
2220 cout << e.simplify_indexed() << endl;
2221 // -> [[3*y+3*x,6*y+2*x]].j
2225 You can of course obtain the same results with the @code{matrix::add()},
2226 @code{matrix::mul()} and @code{matrix::trace()} methods (@pxref{Matrices})
2227 but with indices you don't have to worry about transposing matrices.
2229 Matrix indices always start at 0 and their dimension must match the number
2230 of rows/columns of the matrix. Matrices with one row or one column are
2231 vectors and can have one or two indices (it doesn't matter whether it's a
2232 row or a column vector). Other matrices must have two indices.
2234 You should be careful when using indices with variance on matrices. GiNaC
2235 doesn't look at the variance and doesn't know that @samp{F~mu~nu} and
2236 @samp{F.mu.nu} are different matrices. In this case you should use only
2237 one form for @samp{F} and explicitly multiply it with a matrix representation
2238 of the metric tensor.
2241 @node Non-commutative objects, Methods and Functions, Indexed objects, Basic Concepts
2242 @c node-name, next, previous, up
2243 @section Non-commutative objects
2245 GiNaC is equipped to handle certain non-commutative algebras. Three classes of
2246 non-commutative objects are built-in which are mostly of use in high energy
2250 @item Clifford (Dirac) algebra (class @code{clifford})
2251 @item su(3) Lie algebra (class @code{color})
2252 @item Matrices (unindexed) (class @code{matrix})
2255 The @code{clifford} and @code{color} classes are subclasses of
2256 @code{indexed} because the elements of these algebras usually carry
2257 indices. The @code{matrix} class is described in more detail in
2260 Unlike most computer algebra systems, GiNaC does not primarily provide an
2261 operator (often denoted @samp{&*}) for representing inert products of
2262 arbitrary objects. Rather, non-commutativity in GiNaC is a property of the
2263 classes of objects involved, and non-commutative products are formed with
2264 the usual @samp{*} operator, as are ordinary products. GiNaC is capable of
2265 figuring out by itself which objects commute and will group the factors
2266 by their class. Consider this example:
2270 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
2271 idx a(symbol("a"), 8), b(symbol("b"), 8);
2272 ex e = -dirac_gamma(mu) * (2*color_T(a)) * 8 * color_T(b) * dirac_gamma(nu);
2274 // -> -16*(gamma~mu*gamma~nu)*(T.a*T.b)
2278 As can be seen, GiNaC pulls out the overall commutative factor @samp{-16} and
2279 groups the non-commutative factors (the gammas and the su(3) generators)
2280 together while preserving the order of factors within each class (because
2281 Clifford objects commute with color objects). The resulting expression is a
2282 @emph{commutative} product with two factors that are themselves non-commutative
2283 products (@samp{gamma~mu*gamma~nu} and @samp{T.a*T.b}). For clarification,
2284 parentheses are placed around the non-commutative products in the output.
2286 @cindex @code{ncmul} (class)
2287 Non-commutative products are internally represented by objects of the class
2288 @code{ncmul}, as opposed to commutative products which are handled by the
2289 @code{mul} class. You will normally not have to worry about this distinction,
2292 The advantage of this approach is that you never have to worry about using
2293 (or forgetting to use) a special operator when constructing non-commutative
2294 expressions. Also, non-commutative products in GiNaC are more intelligent
2295 than in other computer algebra systems; they can, for example, automatically
2296 canonicalize themselves according to rules specified in the implementation
2297 of the non-commutative classes. The drawback is that to work with other than
2298 the built-in algebras you have to implement new classes yourself. Symbols
2299 always commute and it's not possible to construct non-commutative products
2300 using symbols to represent the algebra elements or generators. User-defined
2301 functions can, however, be specified as being non-commutative.
2303 @cindex @code{return_type()}
2304 @cindex @code{return_type_tinfo()}
2305 Information about the commutativity of an object or expression can be
2306 obtained with the two member functions
2309 unsigned ex::return_type(void) const;
2310 unsigned ex::return_type_tinfo(void) const;
2313 The @code{return_type()} function returns one of three values (defined in
2314 the header file @file{flags.h}), corresponding to three categories of
2315 expressions in GiNaC:
2318 @item @code{return_types::commutative}: Commutes with everything. Most GiNaC
2319 classes are of this kind.
2320 @item @code{return_types::noncommutative}: Non-commutative, belonging to a
2321 certain class of non-commutative objects which can be determined with the
2322 @code{return_type_tinfo()} method. Expressions of this category commute
2323 with everything except @code{noncommutative} expressions of the same
2325 @item @code{return_types::noncommutative_composite}: Non-commutative, composed
2326 of non-commutative objects of different classes. Expressions of this
2327 category don't commute with any other @code{noncommutative} or
2328 @code{noncommutative_composite} expressions.
2331 The value returned by the @code{return_type_tinfo()} method is valid only
2332 when the return type of the expression is @code{noncommutative}. It is a
2333 value that is unique to the class of the object and usually one of the
2334 constants in @file{tinfos.h}, or derived therefrom.
2336 Here are a couple of examples:
2339 @multitable @columnfractions 0.33 0.33 0.34
2340 @item @strong{Expression} @tab @strong{@code{return_type()}} @tab @strong{@code{return_type_tinfo()}}
2341 @item @code{42} @tab @code{commutative} @tab -
2342 @item @code{2*x-y} @tab @code{commutative} @tab -
2343 @item @code{dirac_ONE()} @tab @code{noncommutative} @tab @code{TINFO_clifford}
2344 @item @code{dirac_gamma(mu)*dirac_gamma(nu)} @tab @code{noncommutative} @tab @code{TINFO_clifford}
2345 @item @code{2*color_T(a)} @tab @code{noncommutative} @tab @code{TINFO_color}
2346 @item @code{dirac_ONE()*color_T(a)} @tab @code{noncommutative_composite} @tab -
2350 Note: the @code{return_type_tinfo()} of Clifford objects is only equal to
2351 @code{TINFO_clifford} for objects with a representation label of zero.
2352 Other representation labels yield a different @code{return_type_tinfo()},
2353 but it's the same for any two objects with the same label. This is also true
2356 A last note: With the exception of matrices, positive integer powers of
2357 non-commutative objects are automatically expanded in GiNaC. For example,
2358 @code{pow(a*b, 2)} becomes @samp{a*b*a*b} if @samp{a} and @samp{b} are
2359 non-commutative expressions).
2362 @cindex @code{clifford} (class)
2363 @subsection Clifford algebra
2365 @cindex @code{dirac_gamma()}
2366 Clifford algebra elements (also called Dirac gamma matrices, although GiNaC
2367 doesn't treat them as matrices) are designated as @samp{gamma~mu} and satisfy
2368 @samp{gamma~mu*gamma~nu + gamma~nu*gamma~mu = 2*eta~mu~nu} where @samp{eta~mu~nu}
2369 is the Minkowski metric tensor. Dirac gammas are constructed by the function
2372 ex dirac_gamma(const ex & mu, unsigned char rl = 0);
2375 which takes two arguments: the index and a @dfn{representation label} in the
2376 range 0 to 255 which is used to distinguish elements of different Clifford
2377 algebras (this is also called a @dfn{spin line index}). Gammas with different
2378 labels commute with each other. The dimension of the index can be 4 or (in
2379 the framework of dimensional regularization) any symbolic value. Spinor
2380 indices on Dirac gammas are not supported in GiNaC.
2382 @cindex @code{dirac_ONE()}
2383 The unity element of a Clifford algebra is constructed by
2386 ex dirac_ONE(unsigned char rl = 0);
2389 @strong{Note:} You must always use @code{dirac_ONE()} when referring to
2390 multiples of the unity element, even though it's customary to omit it.
2391 E.g. instead of @code{dirac_gamma(mu)*(dirac_slash(q,4)+m)} you have to
2392 write @code{dirac_gamma(mu)*(dirac_slash(q,4)+m*dirac_ONE())}. Otherwise,
2393 GiNaC may produce incorrect results.
2395 @cindex @code{dirac_gamma5()}
2396 There's a special element @samp{gamma5} that commutes with all other
2397 gammas and in 4 dimensions equals @samp{gamma~0 gamma~1 gamma~2 gamma~3},
2401 ex dirac_gamma5(unsigned char rl = 0);
2404 @cindex @code{dirac_gamma6()}
2405 @cindex @code{dirac_gamma7()}
2406 The two additional functions
2409 ex dirac_gamma6(unsigned char rl = 0);
2410 ex dirac_gamma7(unsigned char rl = 0);
2413 return @code{dirac_ONE(rl) + dirac_gamma5(rl)} and @code{dirac_ONE(rl) - dirac_gamma5(rl)},
2416 @cindex @code{dirac_slash()}
2417 Finally, the function
2420 ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
2423 creates a term that represents a contraction of @samp{e} with the Dirac
2424 Lorentz vector (it behaves like a term of the form @samp{e.mu gamma~mu}
2425 with a unique index whose dimension is given by the @code{dim} argument).
2426 Such slashed expressions are printed with a trailing backslash, e.g. @samp{e\}.
2428 In products of dirac gammas, superfluous unity elements are automatically
2429 removed, squares are replaced by their values and @samp{gamma5} is
2430 anticommuted to the front. The @code{simplify_indexed()} function performs
2431 contractions in gamma strings, for example
2436 symbol a("a"), b("b"), D("D");
2437 varidx mu(symbol("mu"), D);
2438 ex e = dirac_gamma(mu) * dirac_slash(a, D)
2439 * dirac_gamma(mu.toggle_variance());
2441 // -> gamma~mu*a\*gamma.mu
2442 e = e.simplify_indexed();
2445 cout << e.subs(D == 4) << endl;
2451 @cindex @code{dirac_trace()}
2452 To calculate the trace of an expression containing strings of Dirac gammas
2453 you use the function
2456 ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4);
2459 This function takes the trace of all gammas with the specified representation
2460 label; gammas with other labels are left standing. The last argument to
2461 @code{dirac_trace()} is the value to be returned for the trace of the unity
2462 element, which defaults to 4. The @code{dirac_trace()} function is a linear
2463 functional that is equal to the usual trace only in @math{D = 4} dimensions.
2464 In particular, the functional is not cyclic in @math{D != 4} dimensions when
2465 acting on expressions containing @samp{gamma5}, so it's not a proper trace.
2466 This @samp{gamma5} scheme is described in greater detail in
2467 @cite{The Role of gamma5 in Dimensional Regularization}.
2469 The value of the trace itself is also usually different in 4 and in
2470 @math{D != 4} dimensions:
2475 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
2476 ex e = dirac_gamma(mu) * dirac_gamma(nu) *
2477 dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
2478 cout << dirac_trace(e).simplify_indexed() << endl;
2485 varidx mu(symbol("mu"), D), nu(symbol("nu"), D), rho(symbol("rho"), D);
2486 ex e = dirac_gamma(mu) * dirac_gamma(nu) *
2487 dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
2488 cout << dirac_trace(e).simplify_indexed() << endl;
2489 // -> 8*eta~rho~nu-4*eta~rho~nu*D
2493 Here is an example for using @code{dirac_trace()} to compute a value that
2494 appears in the calculation of the one-loop vacuum polarization amplitude in
2499 symbol q("q"), l("l"), m("m"), ldotq("ldotq"), D("D");
2500 varidx mu(symbol("mu"), D), nu(symbol("nu"), D);
2503 sp.add(l, l, pow(l, 2));
2504 sp.add(l, q, ldotq);
2506 ex e = dirac_gamma(mu) *
2507 (dirac_slash(l, D) + dirac_slash(q, D) + m * dirac_ONE()) *
2508 dirac_gamma(mu.toggle_variance()) *
2509 (dirac_slash(l, D) + m * dirac_ONE());
2510 e = dirac_trace(e).simplify_indexed(sp);
2511 e = e.collect(lst(l, ldotq, m));
2513 // -> (8-4*D)*l^2+(8-4*D)*ldotq+4*D*m^2
2517 The @code{canonicalize_clifford()} function reorders all gamma products that
2518 appear in an expression to a canonical (but not necessarily simple) form.
2519 You can use this to compare two expressions or for further simplifications:
2523 varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
2524 ex e = dirac_gamma(mu) * dirac_gamma(nu) + dirac_gamma(nu) * dirac_gamma(mu);
2526 // -> gamma~mu*gamma~nu+gamma~nu*gamma~mu
2528 e = canonicalize_clifford(e);
2535 @cindex @code{color} (class)
2536 @subsection Color algebra
2538 @cindex @code{color_T()}
2539 For computations in quantum chromodynamics, GiNaC implements the base elements
2540 and structure constants of the su(3) Lie algebra (color algebra). The base
2541 elements @math{T_a} are constructed by the function
2544 ex color_T(const ex & a, unsigned char rl = 0);
2547 which takes two arguments: the index and a @dfn{representation label} in the
2548 range 0 to 255 which is used to distinguish elements of different color
2549 algebras. Objects with different labels commute with each other. The
2550 dimension of the index must be exactly 8 and it should be of class @code{idx},
2553 @cindex @code{color_ONE()}
2554 The unity element of a color algebra is constructed by
2557 ex color_ONE(unsigned char rl = 0);
2560 @strong{Note:} You must always use @code{color_ONE()} when referring to
2561 multiples of the unity element, even though it's customary to omit it.
2562 E.g. instead of @code{color_T(a)*(color_T(b)*indexed(X,b)+1)} you have to
2563 write @code{color_T(a)*(color_T(b)*indexed(X,b)+color_ONE())}. Otherwise,
2564 GiNaC may produce incorrect results.
2566 @cindex @code{color_d()}
2567 @cindex @code{color_f()}
2571 ex color_d(const ex & a, const ex & b, const ex & c);
2572 ex color_f(const ex & a, const ex & b, const ex & c);
2575 create the symmetric and antisymmetric structure constants @math{d_abc} and
2576 @math{f_abc} which satisfy @math{@{T_a, T_b@} = 1/3 delta_ab + d_abc T_c}
2577 and @math{[T_a, T_b] = i f_abc T_c}.
2579 @cindex @code{color_h()}
2580 There's an additional function
2583 ex color_h(const ex & a, const ex & b, const ex & c);
2586 which returns the linear combination @samp{color_d(a, b, c)+I*color_f(a, b, c)}.
2588 The function @code{simplify_indexed()} performs some simplifications on
2589 expressions containing color objects:
2594 idx a(symbol("a"), 8), b(symbol("b"), 8), c(symbol("c"), 8),
2595 k(symbol("k"), 8), l(symbol("l"), 8);
2597 e = color_d(a, b, l) * color_f(a, b, k);
2598 cout << e.simplify_indexed() << endl;
2601 e = color_d(a, b, l) * color_d(a, b, k);
2602 cout << e.simplify_indexed() << endl;
2605 e = color_f(l, a, b) * color_f(a, b, k);
2606 cout << e.simplify_indexed() << endl;
2609 e = color_h(a, b, c) * color_h(a, b, c);
2610 cout << e.simplify_indexed() << endl;
2613 e = color_h(a, b, c) * color_T(b) * color_T(c);
2614 cout << e.simplify_indexed() << endl;
2617 e = color_h(a, b, c) * color_T(a) * color_T(b) * color_T(c);
2618 cout << e.simplify_indexed() << endl;
2621 e = color_T(k) * color_T(a) * color_T(b) * color_T(k);
2622 cout << e.simplify_indexed() << endl;
2623 // -> 1/4*delta.b.a*ONE-1/6*T.a*T.b
2627 @cindex @code{color_trace()}
2628 To calculate the trace of an expression containing color objects you use the
2632 ex color_trace(const ex & e, unsigned char rl = 0);
2635 This function takes the trace of all color @samp{T} objects with the
2636 specified representation label; @samp{T}s with other labels are left
2637 standing. For example:
2641 e = color_trace(4 * color_T(a) * color_T(b) * color_T(c));
2643 // -> -I*f.a.c.b+d.a.c.b
2648 @node Methods and Functions, Information About Expressions, Non-commutative objects, Top
2649 @c node-name, next, previous, up
2650 @chapter Methods and Functions
2653 In this chapter the most important algorithms provided by GiNaC will be
2654 described. Some of them are implemented as functions on expressions,
2655 others are implemented as methods provided by expression objects. If
2656 they are methods, there exists a wrapper function around it, so you can
2657 alternatively call it in a functional way as shown in the simple
2662 cout << "As method: " << sin(1).evalf() << endl;
2663 cout << "As function: " << evalf(sin(1)) << endl;
2667 @cindex @code{subs()}
2668 The general rule is that wherever methods accept one or more parameters
2669 (@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
2670 wrapper accepts is the same but preceded by the object to act on
2671 (@var{object}, @var{arg1}, @var{arg2}, @dots{}). This approach is the
2672 most natural one in an OO model but it may lead to confusion for MapleV
2673 users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
2674 would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
2675 @code{A} and @code{x}). On the other hand, since MapleV returns 3 on
2676 @code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
2677 coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
2678 here. Also, users of MuPAD will in most cases feel more comfortable
2679 with GiNaC's convention. All function wrappers are implemented
2680 as simple inline functions which just call the corresponding method and
2681 are only provided for users uncomfortable with OO who are dead set to
2682 avoid method invocations. Generally, nested function wrappers are much
2683 harder to read than a sequence of methods and should therefore be
2684 avoided if possible. On the other hand, not everything in GiNaC is a
2685 method on class @code{ex} and sometimes calling a function cannot be
2689 * Information About Expressions::
2690 * Substituting Expressions::
2691 * Pattern Matching and Advanced Substitutions::
2692 * Applying a Function on Subexpressions::
2693 * Polynomial Arithmetic:: Working with polynomials.
2694 * Rational Expressions:: Working with rational functions.
2695 * Symbolic Differentiation::
2696 * Series Expansion:: Taylor and Laurent expansion.
2698 * Built-in Functions:: List of predefined mathematical functions.
2699 * Input/Output:: Input and output of expressions.
2703 @node Information About Expressions, Substituting Expressions, Methods and Functions, Methods and Functions
2704 @c node-name, next, previous, up
2705 @section Getting information about expressions
2707 @subsection Checking expression types
2708 @cindex @code{is_a<@dots{}>()}
2709 @cindex @code{is_exactly_a<@dots{}>()}
2710 @cindex @code{ex_to<@dots{}>()}
2711 @cindex Converting @code{ex} to other classes
2712 @cindex @code{info()}
2713 @cindex @code{return_type()}
2714 @cindex @code{return_type_tinfo()}
2716 Sometimes it's useful to check whether a given expression is a plain number,
2717 a sum, a polynomial with integer coefficients, or of some other specific type.
2718 GiNaC provides a couple of functions for this:
2721 bool is_a<T>(const ex & e);
2722 bool is_exactly_a<T>(const ex & e);
2723 bool ex::info(unsigned flag);
2724 unsigned ex::return_type(void) const;
2725 unsigned ex::return_type_tinfo(void) const;
2728 When the test made by @code{is_a<T>()} returns true, it is safe to call
2729 one of the functions @code{ex_to<T>()}, where @code{T} is one of the
2730 class names (@xref{The Class Hierarchy}, for a list of all classes). For
2731 example, assuming @code{e} is an @code{ex}:
2736 if (is_a<numeric>(e))
2737 numeric n = ex_to<numeric>(e);
2742 @code{is_a<T>(e)} allows you to check whether the top-level object of
2743 an expression @samp{e} is an instance of the GiNaC class @samp{T}
2744 (@xref{The Class Hierarchy}, for a list of all classes). This is most useful,
2745 e.g., for checking whether an expression is a number, a sum, or a product:
2752 is_a<numeric>(e1); // true
2753 is_a<numeric>(e2); // false
2754 is_a<add>(e1); // false
2755 is_a<add>(e2); // true
2756 is_a<mul>(e1); // false
2757 is_a<mul>(e2); // false
2761 In contrast, @code{is_exactly_a<T>(e)} allows you to check whether the
2762 top-level object of an expression @samp{e} is an instance of the GiNaC
2763 class @samp{T}, not including parent classes.
2765 The @code{info()} method is used for checking certain attributes of
2766 expressions. The possible values for the @code{flag} argument are defined
2767 in @file{ginac/flags.h}, the most important being explained in the following
2771 @multitable @columnfractions .30 .70
2772 @item @strong{Flag} @tab @strong{Returns true if the object is@dots{}}
2773 @item @code{numeric}
2774 @tab @dots{}a number (same as @code{is_<numeric>(...)})
2776 @tab @dots{}a real integer, rational or float (i.e. is not complex)
2777 @item @code{rational}
2778 @tab @dots{}an exact rational number (integers are rational, too)
2779 @item @code{integer}
2780 @tab @dots{}a (non-complex) integer
2781 @item @code{crational}
2782 @tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
2783 @item @code{cinteger}
2784 @tab @dots{}a (complex) integer (such as @math{2-3*I})
2785 @item @code{positive}
2786 @tab @dots{}not complex and greater than 0
2787 @item @code{negative}
2788 @tab @dots{}not complex and less than 0
2789 @item @code{nonnegative}
2790 @tab @dots{}not complex and greater than or equal to 0
2792 @tab @dots{}an integer greater than 0
2794 @tab @dots{}an integer less than 0
2795 @item @code{nonnegint}
2796 @tab @dots{}an integer greater than or equal to 0
2798 @tab @dots{}an even integer
2800 @tab @dots{}an odd integer
2802 @tab @dots{}a prime integer (probabilistic primality test)
2803 @item @code{relation}
2804 @tab @dots{}a relation (same as @code{is_a<relational>(...)})
2805 @item @code{relation_equal}
2806 @tab @dots{}a @code{==} relation
2807 @item @code{relation_not_equal}
2808 @tab @dots{}a @code{!=} relation
2809 @item @code{relation_less}
2810 @tab @dots{}a @code{<} relation
2811 @item @code{relation_less_or_equal}
2812 @tab @dots{}a @code{<=} relation
2813 @item @code{relation_greater}
2814 @tab @dots{}a @code{>} relation
2815 @item @code{relation_greater_or_equal}
2816 @tab @dots{}a @code{>=} relation
2818 @tab @dots{}a symbol (same as @code{is_a<symbol>(...)})
2820 @tab @dots{}a list (same as @code{is_a<lst>(...)})
2821 @item @code{polynomial}
2822 @tab @dots{}a polynomial (i.e. only consists of sums and products of numbers and symbols with positive integer powers)
2823 @item @code{integer_polynomial}
2824 @tab @dots{}a polynomial with (non-complex) integer coefficients
2825 @item @code{cinteger_polynomial}
2826 @tab @dots{}a polynomial with (possibly complex) integer coefficients (such as @math{2-3*I})
2827 @item @code{rational_polynomial}
2828 @tab @dots{}a polynomial with (non-complex) rational coefficients
2829 @item @code{crational_polynomial}
2830 @tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I})
2831 @item @code{rational_function}
2832 @tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)})
2833 @item @code{algebraic}
2834 @tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1})
2838 To determine whether an expression is commutative or non-commutative and if
2839 so, with which other expressions it would commute, you use the methods
2840 @code{return_type()} and @code{return_type_tinfo()}. @xref{Non-commutative objects},
2841 for an explanation of these.
2844 @subsection Accessing subexpressions
2845 @cindex @code{nops()}
2848 @cindex @code{relational} (class)
2850 GiNaC provides the two methods
2853 unsigned ex::nops();
2854 ex ex::op(unsigned i);
2857 for accessing the subexpressions in the container-like GiNaC classes like
2858 @code{add}, @code{mul}, @code{lst}, and @code{function}. @code{nops()}
2859 determines the number of subexpressions (@samp{operands}) contained, while
2860 @code{op()} returns the @code{i}-th (0..@code{nops()-1}) subexpression.
2861 In the case of a @code{power} object, @code{op(0)} will return the basis
2862 and @code{op(1)} the exponent. For @code{indexed} objects, @code{op(0)}
2863 is the base expression and @code{op(i)}, @math{i>0} are the indices.
2865 The left-hand and right-hand side expressions of objects of class
2866 @code{relational} (and only of these) can also be accessed with the methods
2874 @subsection Comparing expressions
2875 @cindex @code{is_equal()}
2876 @cindex @code{is_zero()}
2878 Expressions can be compared with the usual C++ relational operators like
2879 @code{==}, @code{>}, and @code{<} but if the expressions contain symbols,
2880 the result is usually not determinable and the result will be @code{false},
2881 except in the case of the @code{!=} operator. You should also be aware that
2882 GiNaC will only do the most trivial test for equality (subtracting both
2883 expressions), so something like @code{(pow(x,2)+x)/x==x+1} will return
2886 Actually, if you construct an expression like @code{a == b}, this will be
2887 represented by an object of the @code{relational} class (@pxref{Relations})
2888 which is not evaluated until (explicitly or implicitly) cast to a @code{bool}.
2890 There are also two methods
2893 bool ex::is_equal(const ex & other);
2897 for checking whether one expression is equal to another, or equal to zero,
2900 @strong{Warning:} You will also find an @code{ex::compare()} method in the
2901 GiNaC header files. This method is however only to be used internally by
2902 GiNaC to establish a canonical sort order for terms, and using it to compare
2903 expressions will give very surprising results.
2906 @node Substituting Expressions, Pattern Matching and Advanced Substitutions, Information About Expressions, Methods and Functions
2907 @c node-name, next, previous, up
2908 @section Substituting expressions
2909 @cindex @code{subs()}
2911 Algebraic objects inside expressions can be replaced with arbitrary
2912 expressions via the @code{.subs()} method:
2915 ex ex::subs(const ex & e);
2916 ex ex::subs(const lst & syms, const lst & repls);
2919 In the first form, @code{subs()} accepts a relational of the form
2920 @samp{object == expression} or a @code{lst} of such relationals:
2924 symbol x("x"), y("y");
2926 ex e1 = 2*x^2-4*x+3;
2927 cout << "e1(7) = " << e1.subs(x == 7) << endl;
2931 cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
2936 If you specify multiple substitutions, they are performed in parallel, so e.g.
2937 @code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
2939 The second form of @code{subs()} takes two lists, one for the objects to be
2940 replaced and one for the expressions to be substituted (both lists must
2941 contain the same number of elements). Using this form, you would write
2942 @code{subs(lst(x, y), lst(y, x))} to exchange @samp{x} and @samp{y}.
2944 @code{subs()} performs syntactic substitution of any complete algebraic
2945 object; it does not try to match sub-expressions as is demonstrated by the
2950 symbol x("x"), y("y"), z("z");
2952 ex e1 = pow(x+y, 2);
2953 cout << e1.subs(x+y == 4) << endl;
2956 ex e2 = sin(x)*sin(y)*cos(x);
2957 cout << e2.subs(sin(x) == cos(x)) << endl;
2958 // -> cos(x)^2*sin(y)
2961 cout << e3.subs(x+y == 4) << endl;
2963 // (and not 4+z as one might expect)
2967 A more powerful form of substitution using wildcards is described in the
2971 @node Pattern Matching and Advanced Substitutions, Applying a Function on Subexpressions, Substituting Expressions, Methods and Functions
2972 @c node-name, next, previous, up
2973 @section Pattern matching and advanced substitutions
2974 @cindex @code{wildcard} (class)
2975 @cindex Pattern matching
2977 GiNaC allows the use of patterns for checking whether an expression is of a
2978 certain form or contains subexpressions of a certain form, and for
2979 substituting expressions in a more general way.
2981 A @dfn{pattern} is an algebraic expression that optionally contains wildcards.
2982 A @dfn{wildcard} is a special kind of object (of class @code{wildcard}) that
2983 represents an arbitrary expression. Every wildcard has a @dfn{label} which is
2984 an unsigned integer number to allow having multiple different wildcards in a
2985 pattern. Wildcards are printed as @samp{$label} (this is also the way they
2986 are specified in @command{ginsh}). In C++ code, wildcard objects are created
2990 ex wild(unsigned label = 0);
2993 which is simply a wrapper for the @code{wildcard()} constructor with a shorter
2996 Some examples for patterns:
2998 @multitable @columnfractions .5 .5
2999 @item @strong{Constructed as} @tab @strong{Output as}
3000 @item @code{wild()} @tab @samp{$0}
3001 @item @code{pow(x,wild())} @tab @samp{x^$0}
3002 @item @code{atan2(wild(1),wild(2))} @tab @samp{atan2($1,$2)}
3003 @item @code{indexed(A,idx(wild(),3))} @tab @samp{A.$0}
3009 @item Wildcards behave like symbols and are subject to the same algebraic
3010 rules. E.g., @samp{$0+2*$0} is automatically transformed to @samp{3*$0}.
3011 @item As shown in the last example, to use wildcards for indices you have to
3012 use them as the value of an @code{idx} object. This is because indices must
3013 always be of class @code{idx} (or a subclass).
3014 @item Wildcards only represent expressions or subexpressions. It is not
3015 possible to use them as placeholders for other properties like index
3016 dimension or variance, representation labels, symmetry of indexed objects
3018 @item Because wildcards are commutative, it is not possible to use wildcards
3019 as part of noncommutative products.
3020 @item A pattern does not have to contain wildcards. @samp{x} and @samp{x+y}
3021 are also valid patterns.
3024 @cindex @code{match()}
3025 The most basic application of patterns is to check whether an expression
3026 matches a given pattern. This is done by the function
3029 bool ex::match(const ex & pattern);
3030 bool ex::match(const ex & pattern, lst & repls);
3033 This function returns @code{true} when the expression matches the pattern
3034 and @code{false} if it doesn't. If used in the second form, the actual
3035 subexpressions matched by the wildcards get returned in the @code{repls}
3036 object as a list of relations of the form @samp{wildcard == expression}.
3037 If @code{match()} returns false, the state of @code{repls} is undefined.
3038 For reproducible results, the list should be empty when passed to
3039 @code{match()}, but it is also possible to find similarities in multiple
3040 expressions by passing in the result of a previous match.
3042 The matching algorithm works as follows:
3045 @item A single wildcard matches any expression. If one wildcard appears
3046 multiple times in a pattern, it must match the same expression in all
3047 places (e.g. @samp{$0} matches anything, and @samp{$0*($0+1)} matches
3048 @samp{x*(x+1)} but not @samp{x*(y+1)}).
3049 @item If the expression is not of the same class as the pattern, the match
3050 fails (i.e. a sum only matches a sum, a function only matches a function,
3052 @item If the pattern is a function, it only matches the same function
3053 (i.e. @samp{sin($0)} matches @samp{sin(x)} but doesn't match @samp{exp(x)}).
3054 @item Except for sums and products, the match fails if the number of
3055 subexpressions (@code{nops()}) is not equal to the number of subexpressions
3057 @item If there are no subexpressions, the expressions and the pattern must
3058 be equal (in the sense of @code{is_equal()}).
3059 @item Except for sums and products, each subexpression (@code{op()}) must
3060 match the corresponding subexpression of the pattern.
3063 Sums (@code{add}) and products (@code{mul}) are treated in a special way to
3064 account for their commutativity and associativity:
3067 @item If the pattern contains a term or factor that is a single wildcard,
3068 this one is used as the @dfn{global wildcard}. If there is more than one
3069 such wildcard, one of them is chosen as the global wildcard in a random
3071 @item Every term/factor of the pattern, except the global wildcard, is
3072 matched against every term of the expression in sequence. If no match is
3073 found, the whole match fails. Terms that did match are not considered in
3075 @item If there are no unmatched terms left, the match succeeds. Otherwise
3076 the match fails unless there is a global wildcard in the pattern, in
3077 which case this wildcard matches the remaining terms.
3080 In general, having more than one single wildcard as a term of a sum or a
3081 factor of a product (such as @samp{a+$0+$1}) will lead to unpredictable or
3084 Here are some examples in @command{ginsh} to demonstrate how it works (the
3085 @code{match()} function in @command{ginsh} returns @samp{FAIL} if the
3086 match fails, and the list of wildcard replacements otherwise):
3089 > match((x+y)^a,(x+y)^a);
3091 > match((x+y)^a,(x+y)^b);
3093 > match((x+y)^a,$1^$2);
3095 > match((x+y)^a,$1^$1);
3097 > match((x+y)^(x+y),$1^$1);
3099 > match((x+y)^(x+y),$1^$2);
3101 > match((a+b)*(a+c),($1+b)*($1+c));
3103 > match((a+b)*(a+c),(a+$1)*(a+$2));
3105 (Unpredictable. The result might also be [$1==c,$2==b].)
3106 > match((a+b)*(a+c),($1+$2)*($1+$3));
3107 (The result is undefined. Due to the sequential nature of the algorithm
3108 and the re-ordering of terms in GiNaC, the match for the first factor
3109 may be @{$1==a,$2==b@} in which case the match for the second factor
3110 succeeds, or it may be @{$1==b,$2==a@} which causes the second match to
3112 > match(a*(x+y)+a*z+b,a*$1+$2);
3113 (This is also ambiguous and may return either @{$1==z,$2==a*(x+y)+b@} or
3114 @{$1=x+y,$2=a*z+b@}.)
3115 > match(a+b+c+d+e+f,c);
3117 > match(a+b+c+d+e+f,c+$0);
3119 > match(a+b+c+d+e+f,c+e+$0);
3121 > match(a+b,a+b+$0);
3123 > match(a*b^2,a^$1*b^$2);
3125 (The matching is syntactic, not algebraic, and "a" doesn't match "a^$1"
3126 even though a==a^1.)
3127 > match(x*atan2(x,x^2),$0*atan2($0,$0^2));
3129 > match(atan2(y,x^2),atan2(y,$0));
3133 @cindex @code{has()}
3134 A more general way to look for patterns in expressions is provided by the
3138 bool ex::has(const ex & pattern);
3141 This function checks whether a pattern is matched by an expression itself or
3142 by any of its subexpressions.
3144 Again some examples in @command{ginsh} for illustration (in @command{ginsh},
3145 @code{has()} returns @samp{1} for @code{true} and @samp{0} for @code{false}):
3148 > has(x*sin(x+y+2*a),y);
3150 > has(x*sin(x+y+2*a),x+y);
3152 (This is because in GiNaC, "x+y" is not a subexpression of "x+y+2*a" (which
3153 has the subexpressions "x", "y" and "2*a".)
3154 > has(x*sin(x+y+2*a),x+y+$1);
3156 (But this is possible.)
3157 > has(x*sin(2*(x+y)+2*a),x+y);
3159 (This fails because "2*(x+y)" automatically gets converted to "2*x+2*y" of
3160 which "x+y" is not a subexpression.)
3163 (Although x^1==x and x^0==1, neither "x" nor "1" are actually of the form
3165 > has(4*x^2-x+3,$1*x);
3167 > has(4*x^2+x+3,$1*x);
3169 (Another possible pitfall. The first expression matches because the term
3170 "-x" has the form "(-1)*x" in GiNaC. To check whether a polynomial
3171 contains a linear term you should use the coeff() function instead.)
3174 @cindex @code{find()}
3178 bool ex::find(const ex & pattern, lst & found);
3181 works a bit like @code{has()} but it doesn't stop upon finding the first
3182 match. Instead, it appends all found matches to the specified list. If there
3183 are multiple occurrences of the same expression, it is entered only once to
3184 the list. @code{find()} returns false if no matches were found (in
3185 @command{ginsh}, it returns an empty list):
3188 > find(1+x+x^2+x^3,x);
3190 > find(1+x+x^2+x^3,y);
3192 > find(1+x+x^2+x^3,x^$1);
3194 (Note the absence of "x".)
3195 > expand((sin(x)+sin(y))*(a+b));
3196 sin(y)*a+sin(x)*b+sin(x)*a+sin(y)*b
3201 @cindex @code{subs()}
3202 Probably the most useful application of patterns is to use them for
3203 substituting expressions with the @code{subs()} method. Wildcards can be
3204 used in the search patterns as well as in the replacement expressions, where
3205 they get replaced by the expressions matched by them. @code{subs()} doesn't
3206 know anything about algebra; it performs purely syntactic substitutions.
3211 > subs(a^2+b^2+(x+y)^2,$1^2==$1^3);
3213 > subs(a^4+b^4+(x+y)^4,$1^2==$1^3);
3215 > subs((a+b+c)^2,a+b==x);
3217 > subs((a+b+c)^2,a+b+$1==x+$1);
3219 > subs(a+2*b,a+b==x);
3221 > subs(4*x^3-2*x^2+5*x-1,x==a);
3223 > subs(4*x^3-2*x^2+5*x-1,x^$0==a^$0);
3225 > subs(sin(1+sin(x)),sin($1)==cos($1));
3227 > expand(subs(a*sin(x+y)^2+a*cos(x+y)^2+b,cos($1)^2==1-sin($1)^2));
3231 The last example would be written in C++ in this way:
3235 symbol a("a"), b("b"), x("x"), y("y");
3236 e = a*pow(sin(x+y), 2) + a*pow(cos(x+y), 2) + b;
3237 e = e.subs(pow(cos(wild()), 2) == 1-pow(sin(wild()), 2));
3238 cout << e.expand() << endl;
3244 @node Applying a Function on Subexpressions, Polynomial Arithmetic, Pattern Matching and Advanced Substitutions, Methods and Functions
3245 @c node-name, next, previous, up
3246 @section Applying a Function on Subexpressions
3247 @cindex Tree traversal
3248 @cindex @code{map()}
3250 Sometimes you may want to perform an operation on specific parts of an
3251 expression while leaving the general structure of it intact. An example
3252 of this would be a matrix trace operation: the trace of a sum is the sum
3253 of the traces of the individual terms. That is, the trace should @dfn{map}
3254 on the sum, by applying itself to each of the sum's operands. It is possible
3255 to do this manually which usually results in code like this:
3260 if (is_a<matrix>(e))
3261 return ex_to<matrix>(e).trace();
3262 else if (is_a<add>(e)) @{
3264 for (unsigned i=0; i<e.nops(); i++)
3265 sum += calc_trace(e.op(i));
3267 @} else if (is_a<mul>)(e)) @{
3275 This is, however, slightly inefficient (if the sum is very large it can take
3276 a long time to add the terms one-by-one), and its applicability is limited to
3277 a rather small class of expressions. If @code{calc_trace()} is called with
3278 a relation or a list as its argument, you will probably want the trace to
3279 be taken on both sides of the relation or of all elements of the list.
3281 GiNaC offers the @code{map()} method to aid in the implementation of such
3285 ex ex::map(map_function & f) const;
3286 ex ex::map(ex (*f)(const ex & e)) const;
3289 In the first (preferred) form, @code{map()} takes a function object that
3290 is subclassed from the @code{map_function} class. In the second form, it
3291 takes a pointer to a function that accepts and returns an expression.
3292 @code{map()} constructs a new expression of the same type, applying the
3293 specified function on all subexpressions (in the sense of @code{op()}),
3296 The use of a function object makes it possible to supply more arguments to
3297 the function that is being mapped, or to keep local state information.
3298 The @code{map_function} class declares a virtual function call operator
3299 that you can overload. Here is a sample implementation of @code{calc_trace()}
3300 that uses @code{map()} in a recursive fashion:
3303 struct calc_trace : public map_function @{
3304 ex operator()(const ex &e)
3306 if (is_a<matrix>(e))
3307 return ex_to<matrix>(e).trace();
3308 else if (is_a<mul>(e)) @{
3311 return e.map(*this);
3316 This function object could then be used like this:
3320 ex M = ... // expression with matrices
3321 calc_trace do_trace;
3322 ex tr = do_trace(M);
3326 Here is another example for you to meditate over. It removes quadratic
3327 terms in a variable from an expanded polynomial:
3330 struct map_rem_quad : public map_function @{
3332 map_rem_quad(const ex & var_) : var(var_) @{@}
3334 ex operator()(const ex & e)
3336 if (is_a<add>(e) || is_a<mul>(e))
3337 return e.map(*this);
3338 else if (is_a<power>(e) &&
3339 e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
3349 symbol x("x"), y("y");
3352 for (int i=0; i<8; i++)
3353 e += pow(x, i) * pow(y, 8-i) * (i+1);
3355 // -> 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
3357 map_rem_quad rem_quad(x);
3358 cout << rem_quad(e) << endl;
3359 // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
3363 @command{ginsh} offers a slightly different implementation of @code{map()}
3364 that allows applying algebraic functions to operands. The second argument
3365 to @code{map()} is an expression containing the wildcard @samp{$0} which
3366 acts as the placeholder for the operands:
3371 > map(a+2*b,sin($0));
3373 > map(@{a,b,c@},$0^2+$0);
3374 @{a^2+a,b^2+b,c^2+c@}
3377 Note that it is only possible to use algebraic functions in the second
3378 argument. You can not use functions like @samp{diff()}, @samp{op()},
3379 @samp{subs()} etc. because these are evaluated immediately:
3382 > map(@{a,b,c@},diff($0,a));
3384 This is because "diff($0,a)" evaluates to "0", so the command is equivalent
3385 to "map(@{a,b,c@},0)".
3389 @node Polynomial Arithmetic, Rational Expressions, Applying a Function on Subexpressions, Methods and Functions
3390 @c node-name, next, previous, up
3391 @section Polynomial arithmetic
3393 @subsection Expanding and collecting
3394 @cindex @code{expand()}
3395 @cindex @code{collect()}
3397 A polynomial in one or more variables has many equivalent
3398 representations. Some useful ones serve a specific purpose. Consider
3399 for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
3400 21*y*z + 4*z^2} (written down here in output-style). It is equivalent
3401 to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
3402 representations are the recursive ones where one collects for exponents
3403 in one of the three variable. Since the factors are themselves
3404 polynomials in the remaining two variables the procedure can be
3405 repeated. In our example, two possibilities would be @math{(4*y + z)*x
3406 + 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
3409 To bring an expression into expanded form, its method
3415 may be called. In our example above, this corresponds to @math{4*x*y +
3416 x*z + 20*y^2 + 21*y*z + 4*z^2}. Again, since the canonical form in
3417 GiNaC is not easily guessable you should be prepared to see different
3418 orderings of terms in such sums!
3420 Another useful representation of multivariate polynomials is as a
3421 univariate polynomial in one of the variables with the coefficients
3422 being polynomials in the remaining variables. The method
3423 @code{collect()} accomplishes this task:
3426 ex ex::collect(const ex & s, bool distributed = false);
3429 The first argument to @code{collect()} can also be a list of objects in which
3430 case the result is either a recursively collected polynomial, or a polynomial
3431 in a distributed form with terms like @math{c*x1^e1*...*xn^en}, as specified
3432 by the @code{distributed} flag.
3434 Note that the original polynomial needs to be in expanded form (for the
3435 variables concerned) in order for @code{collect()} to be able to find the
3436 coefficients properly.
3438 The following @command{ginsh} transcript shows an application of @code{collect()}
3439 together with @code{find()}:
3442 > a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
3443 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
3444 > collect(a,@{p,q@});
3445 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)
3446 > collect(a,find(a,sin($1)));
3447 (1+q+d+q*d+d*p+p)*sin(y)+(1+q+d+q*d+d*p+p)*sin(x)
3448 > collect(a,@{find(a,sin($1)),p,q@});
3449 (1+(1+d)*p+d+q*(1+d))*sin(x)+(1+(1+d)*p+d+q*(1+d))*sin(y)
3450 > collect(a,@{find(a,sin($1)),d@});
3451 (1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
3454 @subsection Degree and coefficients
3455 @cindex @code{degree()}
3456 @cindex @code{ldegree()}
3457 @cindex @code{coeff()}
3459 The degree and low degree of a polynomial can be obtained using the two
3463 int ex::degree(const ex & s);
3464 int ex::ldegree(const ex & s);
3467 which also work reliably on non-expanded input polynomials (they even work
3468 on rational functions, returning the asymptotic degree). To extract
3469 a coefficient with a certain power from an expanded polynomial you use
3472 ex ex::coeff(const ex & s, int n);
3475 You can also obtain the leading and trailing coefficients with the methods
3478 ex ex::lcoeff(const ex & s);
3479 ex ex::tcoeff(const ex & s);
3482 which are equivalent to @code{coeff(s, degree(s))} and @code{coeff(s, ldegree(s))},
3485 An application is illustrated in the next example, where a multivariate
3486 polynomial is analyzed:
3490 symbol x("x"), y("y");
3491 ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
3492 - pow(x+y,2) + 2*pow(y+2,2) - 8;
3493 ex Poly = PolyInp.expand();
3495 for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) @{
3496 cout << "The x^" << i << "-coefficient is "
3497 << Poly.coeff(x,i) << endl;
3499 cout << "As polynomial in y: "
3500 << Poly.collect(y) << endl;
3504 When run, it returns an output in the following fashion:
3507 The x^0-coefficient is y^2+11*y
3508 The x^1-coefficient is 5*y^2-2*y
3509 The x^2-coefficient is -1
3510 The x^3-coefficient is 4*y
3511 As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
3514 As always, the exact output may vary between different versions of GiNaC
3515 or even from run to run since the internal canonical ordering is not
3516 within the user's sphere of influence.
3518 @code{degree()}, @code{ldegree()}, @code{coeff()}, @code{lcoeff()},
3519 @code{tcoeff()} and @code{collect()} can also be used to a certain degree
3520 with non-polynomial expressions as they not only work with symbols but with
3521 constants, functions and indexed objects as well:
3525 symbol a("a"), b("b"), c("c");
3526 idx i(symbol("i"), 3);
3528 ex e = pow(sin(x) - cos(x), 4);
3529 cout << e.degree(cos(x)) << endl;
3531 cout << e.expand().coeff(sin(x), 3) << endl;
3534 e = indexed(a+b, i) * indexed(b+c, i);
3535 e = e.expand(expand_options::expand_indexed);
3536 cout << e.collect(indexed(b, i)) << endl;
3537 // -> a.i*c.i+(a.i+c.i)*b.i+b.i^2
3542 @subsection Polynomial division
3543 @cindex polynomial division
3546 @cindex pseudo-remainder
3547 @cindex @code{quo()}
3548 @cindex @code{rem()}
3549 @cindex @code{prem()}
3550 @cindex @code{divide()}
3555 ex quo(const ex & a, const ex & b, const symbol & x);
3556 ex rem(const ex & a, const ex & b, const symbol & x);
3559 compute the quotient and remainder of univariate polynomials in the variable
3560 @samp{x}. The results satisfy @math{a = b*quo(a, b, x) + rem(a, b, x)}.
3562 The additional function
3565 ex prem(const ex & a, const ex & b, const symbol & x);
3568 computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
3569 @math{c*a = b*q + prem(a, b, x)}, where @math{c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1)}.
3571 Exact division of multivariate polynomials is performed by the function
3574 bool divide(const ex & a, const ex & b, ex & q);
3577 If @samp{b} divides @samp{a} over the rationals, this function returns @code{true}
3578 and returns the quotient in the variable @code{q}. Otherwise it returns @code{false}
3579 in which case the value of @code{q} is undefined.
3582 @subsection Unit, content and primitive part
3583 @cindex @code{unit()}
3584 @cindex @code{content()}
3585 @cindex @code{primpart()}
3590 ex ex::unit(const symbol & x);
3591 ex ex::content(const symbol & x);
3592 ex ex::primpart(const symbol & x);
3595 return the unit part, content part, and primitive polynomial of a multivariate
3596 polynomial with respect to the variable @samp{x} (the unit part being the sign
3597 of the leading coefficient, the content part being the GCD of the coefficients,
3598 and the primitive polynomial being the input polynomial divided by the unit and
3599 content parts). The product of unit, content, and primitive part is the
3600 original polynomial.
3603 @subsection GCD and LCM
3606 @cindex @code{gcd()}
3607 @cindex @code{lcm()}
3609 The functions for polynomial greatest common divisor and least common
3610 multiple have the synopsis
3613 ex gcd(const ex & a, const ex & b);
3614 ex lcm(const ex & a, const ex & b);
3617 The functions @code{gcd()} and @code{lcm()} accept two expressions
3618 @code{a} and @code{b} as arguments and return a new expression, their
3619 greatest common divisor or least common multiple, respectively. If the
3620 polynomials @code{a} and @code{b} are coprime @code{gcd(a,b)} returns 1
3621 and @code{lcm(a,b)} returns the product of @code{a} and @code{b}.
3624 #include <ginac/ginac.h>
3625 using namespace GiNaC;
3629 symbol x("x"), y("y"), z("z");
3630 ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
3631 ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
3633 ex P_gcd = gcd(P_a, P_b);
3635 ex P_lcm = lcm(P_a, P_b);
3636 // 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
3641 @subsection Square-free decomposition
3642 @cindex square-free decomposition
3643 @cindex factorization
3644 @cindex @code{sqrfree()}
3646 GiNaC still lacks proper factorization support. Some form of
3647 factorization is, however, easily implemented by noting that factors
3648 appearing in a polynomial with power two or more also appear in the
3649 derivative and hence can easily be found by computing the GCD of the
3650 original polynomial and its derivatives. Any decent system has an
3651 interface for this so called square-free factorization. So we provide
3654 ex sqrfree(const ex & a, const lst & l = lst());
3656 Here is an example that by the way illustrates how the exact form of the
3657 result may slightly depend on the order of differentiation, calling for
3658 some care with subsequent processing of the result:
3661 symbol x("x"), y("y");
3662 ex BiVarPol = expand(pow(2-2*y,3) * pow(1+x*y,2) * pow(x-2*y,2) * (x+y));
3664 cout << sqrfree(BiVarPol, lst(x,y)) << endl;
3665 // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)
3667 cout << sqrfree(BiVarPol, lst(y,x)) << endl;
3668 // -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x)
3670 cout << sqrfree(BiVarPol) << endl;
3671 // -> depending on luck, any of the above
3674 Note also, how factors with the same exponents are not fully factorized
3678 @node Rational Expressions, Symbolic Differentiation, Polynomial Arithmetic, Methods and Functions
3679 @c node-name, next, previous, up
3680 @section Rational expressions
3682 @subsection The @code{normal} method
3683 @cindex @code{normal()}
3684 @cindex simplification
3685 @cindex temporary replacement
3687 Some basic form of simplification of expressions is called for frequently.
3688 GiNaC provides the method @code{.normal()}, which converts a rational function
3689 into an equivalent rational function of the form @samp{numerator/denominator}
3690 where numerator and denominator are coprime. If the input expression is already
3691 a fraction, it just finds the GCD of numerator and denominator and cancels it,
3692 otherwise it performs fraction addition and multiplication.
3694 @code{.normal()} can also be used on expressions which are not rational functions
3695 as it will replace all non-rational objects (like functions or non-integer
3696 powers) by temporary symbols to bring the expression to the domain of rational
3697 functions before performing the normalization, and re-substituting these
3698 symbols afterwards. This algorithm is also available as a separate method
3699 @code{.to_rational()}, described below.
3701 This means that both expressions @code{t1} and @code{t2} are indeed
3702 simplified in this little code snippet:
3707 ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
3708 ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
3709 std::cout << "t1 is " << t1.normal() << std::endl;
3710 std::cout << "t2 is " << t2.normal() << std::endl;
3714 Of course this works for multivariate polynomials too, so the ratio of
3715 the sample-polynomials from the section about GCD and LCM above would be
3716 normalized to @code{P_a/P_b} = @code{(4*y+z)/(y+3*z)}.
3719 @subsection Numerator and denominator
3722 @cindex @code{numer()}
3723 @cindex @code{denom()}
3724 @cindex @code{numer_denom()}
3726 The numerator and denominator of an expression can be obtained with
3731 ex ex::numer_denom();
3734 These functions will first normalize the expression as described above and
3735 then return the numerator, denominator, or both as a list, respectively.
3736 If you need both numerator and denominator, calling @code{numer_denom()} is
3737 faster than using @code{numer()} and @code{denom()} separately.
3740 @subsection Converting to a rational expression
3741 @cindex @code{to_rational()}
3743 Some of the methods described so far only work on polynomials or rational
3744 functions. GiNaC provides a way to extend the domain of these functions to
3745 general expressions by using the temporary replacement algorithm described
3746 above. You do this by calling
3749 ex ex::to_rational(lst &l);
3752 on the expression to be converted. The supplied @code{lst} will be filled
3753 with the generated temporary symbols and their replacement expressions in
3754 a format that can be used directly for the @code{subs()} method. It can also
3755 already contain a list of replacements from an earlier application of
3756 @code{.to_rational()}, so it's possible to use it on multiple expressions
3757 and get consistent results.
3764 ex a = pow(sin(x), 2) - pow(cos(x), 2);
3765 ex b = sin(x) + cos(x);
3768 divide(a.to_rational(l), b.to_rational(l), q);
3769 cout << q.subs(l) << endl;
3773 will print @samp{sin(x)-cos(x)}.
3776 @node Symbolic Differentiation, Series Expansion, Rational Expressions, Methods and Functions
3777 @c node-name, next, previous, up
3778 @section Symbolic differentiation
3779 @cindex differentiation
3780 @cindex @code{diff()}
3782 @cindex product rule
3784 GiNaC's objects know how to differentiate themselves. Thus, a
3785 polynomial (class @code{add}) knows that its derivative is the sum of
3786 the derivatives of all the monomials:
3790 symbol x("x"), y("y"), z("z");
3791 ex P = pow(x, 5) + pow(x, 2) + y;
3793 cout << P.diff(x,2) << endl;
3795 cout << P.diff(y) << endl; // 1
3797 cout << P.diff(z) << endl; // 0
3802 If a second integer parameter @var{n} is given, the @code{diff} method
3803 returns the @var{n}th derivative.
3805 If @emph{every} object and every function is told what its derivative
3806 is, all derivatives of composed objects can be calculated using the
3807 chain rule and the product rule. Consider, for instance the expression
3808 @code{1/cosh(x)}. Since the derivative of @code{cosh(x)} is
3809 @code{sinh(x)} and the derivative of @code{pow(x,-1)} is
3810 @code{-pow(x,-2)}, GiNaC can readily compute the composition. It turns
3811 out that the composition is the generating function for Euler Numbers,
3812 i.e. the so called @var{n}th Euler number is the coefficient of
3813 @code{x^n/n!} in the expansion of @code{1/cosh(x)}. We may use this
3814 identity to code a function that generates Euler numbers in just three
3817 @cindex Euler numbers
3819 #include <ginac/ginac.h>
3820 using namespace GiNaC;
3822 ex EulerNumber(unsigned n)
3825 const ex generator = pow(cosh(x),-1);
3826 return generator.diff(x,n).subs(x==0);
3831 for (unsigned i=0; i<11; i+=2)
3832 std::cout << EulerNumber(i) << std::endl;
3837 When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
3838 @code{-61}, @code{1385}, @code{-50521}. We increment the loop variable
3839 @code{i} by two since all odd Euler numbers vanish anyways.
3842 @node Series Expansion, Symmetrization, Symbolic Differentiation, Methods and Functions
3843 @c node-name, next, previous, up
3844 @section Series expansion
3845 @cindex @code{series()}
3846 @cindex Taylor expansion
3847 @cindex Laurent expansion
3848 @cindex @code{pseries} (class)
3849 @cindex @code{Order()}
3851 Expressions know how to expand themselves as a Taylor series or (more
3852 generally) a Laurent series. As in most conventional Computer Algebra
3853 Systems, no distinction is made between those two. There is a class of
3854 its own for storing such series (@code{class pseries}) and a built-in
3855 function (called @code{Order}) for storing the order term of the series.
3856 As a consequence, if you want to work with series, i.e. multiply two
3857 series, you need to call the method @code{ex::series} again to convert
3858 it to a series object with the usual structure (expansion plus order
3859 term). A sample application from special relativity could read:
3862 #include <ginac/ginac.h>
3863 using namespace std;
3864 using namespace GiNaC;
3868 symbol v("v"), c("c");
3870 ex gamma = 1/sqrt(1 - pow(v/c,2));
3871 ex mass_nonrel = gamma.series(v==0, 10);
3873 cout << "the relativistic mass increase with v is " << endl
3874 << mass_nonrel << endl;
3876 cout << "the inverse square of this series is " << endl
3877 << pow(mass_nonrel,-2).series(v==0, 10) << endl;
3881 Only calling the series method makes the last output simplify to
3882 @math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
3883 series raised to the power @math{-2}.
3885 @cindex M@'echain's formula
3886 As another instructive application, let us calculate the numerical
3887 value of Archimedes' constant
3891 (for which there already exists the built-in constant @code{Pi})
3892 using M@'echain's amazing formula
3894 $\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
3897 @math{Pi==16*atan(1/5)-4*atan(1/239)}.
3899 We may expand the arcus tangent around @code{0} and insert the fractions
3900 @code{1/5} and @code{1/239}. But, as we have seen, a series in GiNaC
3901 carries an order term with it and the question arises what the system is
3902 supposed to do when the fractions are plugged into that order term. The
3903 solution is to use the function @code{series_to_poly()} to simply strip
3907 #include <ginac/ginac.h>
3908 using namespace GiNaC;
3910 ex mechain_pi(int degr)
3913 ex pi_expansion = series_to_poly(atan(x).series(x,degr));
3914 ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
3915 -4*pi_expansion.subs(x==numeric(1,239));
3921 using std::cout; // just for fun, another way of...
3922 using std::endl; // ...dealing with this namespace std.
3924 for (int i=2; i<12; i+=2) @{
3925 pi_frac = mechain_pi(i);
3926 cout << i << ":\t" << pi_frac << endl
3927 << "\t" << pi_frac.evalf() << endl;
3933 Note how we just called @code{.series(x,degr)} instead of
3934 @code{.series(x==0,degr)}. This is a simple shortcut for @code{ex}'s
3935 method @code{series()}: if the first argument is a symbol the expression
3936 is expanded in that symbol around point @code{0}. When you run this
3937 program, it will type out:
3941 3.1832635983263598326
3942 4: 5359397032/1706489875
3943 3.1405970293260603143
3944 6: 38279241713339684/12184551018734375
3945 3.141621029325034425
3946 8: 76528487109180192540976/24359780855939418203125
3947 3.141591772182177295
3948 10: 327853873402258685803048818236/104359128170408663038552734375
3949 3.1415926824043995174
3953 @node Symmetrization, Built-in Functions, Series Expansion, Methods and Functions
3954 @c node-name, next, previous, up
3955 @section Symmetrization
3956 @cindex @code{symmetrize()}
3957 @cindex @code{antisymmetrize()}
3958 @cindex @code{symmetrize_cyclic()}
3963 ex ex::symmetrize(const lst & l);
3964 ex ex::antisymmetrize(const lst & l);
3965 ex ex::symmetrize_cyclic(const lst & l);
3968 symmetrize an expression by returning the sum over all symmetric,
3969 antisymmetric or cyclic permutations of the specified list of objects,
3970 weighted by the number of permutations.
3972 The three additional methods
3975 ex ex::symmetrize();
3976 ex ex::antisymmetrize();
3977 ex ex::symmetrize_cyclic();
3980 symmetrize or antisymmetrize an expression over its free indices.
3982 Symmetrization is most useful with indexed expressions but can be used with
3983 almost any kind of object (anything that is @code{subs()}able):
3987 idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
3988 symbol A("A"), B("B"), a("a"), b("b"), c("c");
3990 cout << indexed(A, i, j).symmetrize() << endl;
3991 // -> 1/2*A.j.i+1/2*A.i.j
3992 cout << indexed(A, i, j, k).antisymmetrize(lst(i, j)) << endl;
3993 // -> -1/2*A.j.i.k+1/2*A.i.j.k
3994 cout << lst(a, b, c).symmetrize_cyclic(lst(a, b, c)) << endl;
3995 // -> 1/3*@{a,b,c@}+1/3*@{b,c,a@}+1/3*@{c,a,b@}
4000 @node Built-in Functions, Input/Output, Symmetrization, Methods and Functions
4001 @c node-name, next, previous, up
4002 @section Predefined mathematical functions
4004 GiNaC contains the following predefined mathematical functions:
4007 @multitable @columnfractions .30 .70
4008 @item @strong{Name} @tab @strong{Function}
4011 @cindex @code{abs()}
4012 @item @code{csgn(x)}
4014 @cindex @code{csgn()}
4015 @item @code{sqrt(x)}
4016 @tab square root (not a GiNaC function, rather an alias for @code{pow(x, numeric(1, 2))})
4017 @cindex @code{sqrt()}
4020 @cindex @code{sin()}
4023 @cindex @code{cos()}
4026 @cindex @code{tan()}
4027 @item @code{asin(x)}
4029 @cindex @code{asin()}
4030 @item @code{acos(x)}
4032 @cindex @code{acos()}
4033 @item @code{atan(x)}
4034 @tab inverse tangent
4035 @cindex @code{atan()}
4036 @item @code{atan2(y, x)}
4037 @tab inverse tangent with two arguments
4038 @item @code{sinh(x)}
4039 @tab hyperbolic sine
4040 @cindex @code{sinh()}
4041 @item @code{cosh(x)}
4042 @tab hyperbolic cosine
4043 @cindex @code{cosh()}
4044 @item @code{tanh(x)}
4045 @tab hyperbolic tangent
4046 @cindex @code{tanh()}
4047 @item @code{asinh(x)}
4048 @tab inverse hyperbolic sine
4049 @cindex @code{asinh()}
4050 @item @code{acosh(x)}
4051 @tab inverse hyperbolic cosine
4052 @cindex @code{acosh()}
4053 @item @code{atanh(x)}
4054 @tab inverse hyperbolic tangent
4055 @cindex @code{atanh()}
4057 @tab exponential function
4058 @cindex @code{exp()}
4060 @tab natural logarithm
4061 @cindex @code{log()}
4064 @cindex @code{Li2()}
4065 @item @code{zeta(x)}
4066 @tab Riemann's zeta function
4067 @cindex @code{zeta()}
4068 @item @code{zeta(n, x)}
4069 @tab derivatives of Riemann's zeta function
4070 @item @code{tgamma(x)}
4072 @cindex @code{tgamma()}
4073 @cindex Gamma function
4074 @item @code{lgamma(x)}
4075 @tab logarithm of Gamma function
4076 @cindex @code{lgamma()}
4077 @item @code{beta(x, y)}
4078 @tab Beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
4079 @cindex @code{beta()}
4081 @tab psi (digamma) function
4082 @cindex @code{psi()}
4083 @item @code{psi(n, x)}
4084 @tab derivatives of psi function (polygamma functions)
4085 @item @code{factorial(n)}
4086 @tab factorial function
4087 @cindex @code{factorial()}
4088 @item @code{binomial(n, m)}
4089 @tab binomial coefficients
4090 @cindex @code{binomial()}
4091 @item @code{Order(x)}
4092 @tab order term function in truncated power series
4093 @cindex @code{Order()}
4098 For functions that have a branch cut in the complex plane GiNaC follows
4099 the conventions for C++ as defined in the ANSI standard as far as
4100 possible. In particular: the natural logarithm (@code{log}) and the
4101 square root (@code{sqrt}) both have their branch cuts running along the
4102 negative real axis where the points on the axis itself belong to the
4103 upper part (i.e. continuous with quadrant II). The inverse
4104 trigonometric and hyperbolic functions are not defined for complex
4105 arguments by the C++ standard, however. In GiNaC we follow the
4106 conventions used by CLN, which in turn follow the carefully designed
4107 definitions in the Common Lisp standard. It should be noted that this
4108 convention is identical to the one used by the C99 standard and by most
4109 serious CAS. It is to be expected that future revisions of the C++
4110 standard incorporate these functions in the complex domain in a manner
4111 compatible with C99.
4114 @node Input/Output, Extending GiNaC, Built-in Functions, Methods and Functions
4115 @c node-name, next, previous, up
4116 @section Input and output of expressions
4119 @subsection Expression output
4121 @cindex output of expressions
4123 The easiest way to print an expression is to write it to a stream:
4128 ex e = 4.5+pow(x,2)*3/2;
4129 cout << e << endl; // prints '(4.5)+3/2*x^2'
4133 The output format is identical to the @command{ginsh} input syntax and
4134 to that used by most computer algebra systems, but not directly pastable
4135 into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
4136 is printed as @samp{x^2}).
4138 It is possible to print expressions in a number of different formats with
4142 void ex::print(const print_context & c, unsigned level = 0);
4145 @cindex @code{print_context} (class)
4146 The type of @code{print_context} object passed in determines the format
4147 of the output. The possible types are defined in @file{ginac/print.h}.
4148 All constructors of @code{print_context} and derived classes take an
4149 @code{ostream &} as their first argument.
4151 To print an expression in a way that can be directly used in a C or C++
4152 program, you pass a @code{print_csrc} object like this:
4156 cout << "float f = ";
4157 e.print(print_csrc_float(cout));
4160 cout << "double d = ";
4161 e.print(print_csrc_double(cout));
4164 cout << "cl_N n = ";
4165 e.print(print_csrc_cl_N(cout));
4170 The three possible types mostly affect the way in which floating point
4171 numbers are written.
4173 The above example will produce (note the @code{x^2} being converted to @code{x*x}):
4176 float f = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
4177 double d = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
4178 cl_N n = (cln::cl_F("3.0")/cln::cl_F("2.0"))*(x*x)+cln::cl_F("4.5");
4181 The @code{print_context} type @code{print_tree} provides a dump of the
4182 internal structure of an expression for debugging purposes:
4186 e.print(print_tree(cout));
4193 add, hash=0x0, flags=0x3, nops=2
4194 power, hash=0x9, flags=0x3, nops=2
4195 x (symbol), serial=3, hash=0x44a113a6, flags=0xf
4196 2 (numeric), hash=0x80000042, flags=0xf
4197 3/2 (numeric), hash=0x80000061, flags=0xf
4200 4.5L0 (numeric), hash=0x8000004b, flags=0xf
4204 This kind of output is also available in @command{ginsh} as the @code{print()}
4207 Another useful output format is for LaTeX parsing in mathematical mode.
4208 It is rather similar to the default @code{print_context} but provides
4209 some braces needed by LaTeX for delimiting boxes and also converts some
4210 common objects to conventional LaTeX names. It is possible to give symbols
4211 a special name for LaTeX output by supplying it as a second argument to
4212 the @code{symbol} constructor.
4214 For example, the code snippet
4219 ex foo = lgamma(x).series(x==0,3);
4220 foo.print(print_latex(std::cout));
4226 @{(-\ln(x))@}+@{(-\gamma_E)@} x+@{(1/12 \pi^2)@} x^@{2@}+\mathcal@{O@}(x^3)
4229 @cindex Tree traversal
4230 If you need any fancy special output format, e.g. for interfacing GiNaC
4231 with other algebra systems or for producing code for different
4232 programming languages, you can always traverse the expression tree yourself:
4235 static void my_print(const ex & e)
4237 if (is_a<function>(e))
4238 cout << ex_to<function>(e).get_name();
4240 cout << e.bp->class_name();
4242 unsigned n = e.nops();
4244 for (unsigned i=0; i<n; i++) @{
4256 my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
4264 add(power(numeric(3),symbol(x)),mul(sin(mul(power(constant(Pi),numeric(-1)),
4265 symbol(y))),numeric(-2)))
4268 If you need an output format that makes it possible to accurately
4269 reconstruct an expression by feeding the output to a suitable parser or
4270 object factory, you should consider storing the expression in an
4271 @code{archive} object and reading the object properties from there.
4272 See the section on archiving for more information.
4275 @subsection Expression input
4276 @cindex input of expressions
4278 GiNaC provides no way to directly read an expression from a stream because
4279 you will usually want the user to be able to enter something like @samp{2*x+sin(y)}
4280 and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
4281 @code{y} you defined in your program and there is no way to specify the
4282 desired symbols to the @code{>>} stream input operator.
4284 Instead, GiNaC lets you construct an expression from a string, specifying the
4285 list of symbols to be used:
4289 symbol x("x"), y("y");
4290 ex e("2*x+sin(y)", lst(x, y));
4294 The input syntax is the same as that used by @command{ginsh} and the stream
4295 output operator @code{<<}. The symbols in the string are matched by name to
4296 the symbols in the list and if GiNaC encounters a symbol not specified in
4297 the list it will throw an exception.
4299 With this constructor, it's also easy to implement interactive GiNaC programs:
4304 #include <stdexcept>
4305 #include <ginac/ginac.h>
4306 using namespace std;
4307 using namespace GiNaC;
4314 cout << "Enter an expression containing 'x': ";
4319 cout << "The derivative of " << e << " with respect to x is ";
4320 cout << e.diff(x) << ".\n";
4321 @} catch (exception &p) @{
4322 cerr << p.what() << endl;
4328 @subsection Archiving
4329 @cindex @code{archive} (class)
4332 GiNaC allows creating @dfn{archives} of expressions which can be stored
4333 to or retrieved from files. To create an archive, you declare an object
4334 of class @code{archive} and archive expressions in it, giving each
4335 expression a unique name:
4339 using namespace std;
4340 #include <ginac/ginac.h>
4341 using namespace GiNaC;
4345 symbol x("x"), y("y"), z("z");
4347 ex foo = sin(x + 2*y) + 3*z + 41;
4351 a.archive_ex(foo, "foo");
4352 a.archive_ex(bar, "the second one");
4356 The archive can then be written to a file:
4360 ofstream out("foobar.gar");
4366 The file @file{foobar.gar} contains all information that is needed to
4367 reconstruct the expressions @code{foo} and @code{bar}.
4369 @cindex @command{viewgar}
4370 The tool @command{viewgar} that comes with GiNaC can be used to view
4371 the contents of GiNaC archive files:
4374 $ viewgar foobar.gar
4375 foo = 41+sin(x+2*y)+3*z
4376 the second one = 42+sin(x+2*y)+3*z
4379 The point of writing archive files is of course that they can later be
4385 ifstream in("foobar.gar");
4390 And the stored expressions can be retrieved by their name:
4396 ex ex1 = a2.unarchive_ex(syms, "foo");
4397 ex ex2 = a2.unarchive_ex(syms, "the second one");
4399 cout << ex1 << endl; // prints "41+sin(x+2*y)+3*z"
4400 cout << ex2 << endl; // prints "42+sin(x+2*y)+3*z"
4401 cout << ex1.subs(x == 2) << endl; // prints "41+sin(2+2*y)+3*z"
4405 Note that you have to supply a list of the symbols which are to be inserted
4406 in the expressions. Symbols in archives are stored by their name only and
4407 if you don't specify which symbols you have, unarchiving the expression will
4408 create new symbols with that name. E.g. if you hadn't included @code{x} in
4409 the @code{syms} list above, the @code{ex1.subs(x == 2)} statement would
4410 have had no effect because the @code{x} in @code{ex1} would have been a
4411 different symbol than the @code{x} which was defined at the beginning of
4412 the program, although both would appear as @samp{x} when printed.
4414 You can also use the information stored in an @code{archive} object to
4415 output expressions in a format suitable for exact reconstruction. The
4416 @code{archive} and @code{archive_node} classes have a couple of member
4417 functions that let you access the stored properties:
4420 static void my_print2(const archive_node & n)
4423 n.find_string("class", class_name);
4424 cout << class_name << "(";
4426 archive_node::propinfovector p;
4427 n.get_properties(p);
4429 unsigned num = p.size();
4430 for (unsigned i=0; i<num; i++) @{
4431 const string &name = p[i].name;
4432 if (name == "class")
4434 cout << name << "=";
4436 unsigned count = p[i].count;
4440 for (unsigned j=0; j<count; j++) @{
4441 switch (p[i].type) @{
4442 case archive_node::PTYPE_BOOL: @{
4444 n.find_bool(name, x, j);
4445 cout << (x ? "true" : "false");
4448 case archive_node::PTYPE_UNSIGNED: @{
4450 n.find_unsigned(name, x, j);
4454 case archive_node::PTYPE_STRING: @{
4456 n.find_string(name, x, j);
4457 cout << '\"' << x << '\"';
4460 case archive_node::PTYPE_NODE: @{
4461 const archive_node &x = n.find_ex_node(name, j);
4483 ex e = pow(2, x) - y;
4485 my_print2(ar.get_top_node(0)); cout << endl;
4493 add(rest=@{power(basis=numeric(number="2"),exponent=symbol(name="x")),
4494 symbol(name="y")@},coeff=@{numeric(number="1"),numeric(number="-1")@},
4495 overall_coeff=numeric(number="0"))
4498 Be warned, however, that the set of properties and their meaning for each
4499 class may change between GiNaC versions.
4502 @node Extending GiNaC, What does not belong into GiNaC, Input/Output, Top
4503 @c node-name, next, previous, up
4504 @chapter Extending GiNaC
4506 By reading so far you should have gotten a fairly good understanding of
4507 GiNaC's design-patterns. From here on you should start reading the
4508 sources. All we can do now is issue some recommendations how to tackle
4509 GiNaC's many loose ends in order to fulfill everybody's dreams. If you
4510 develop some useful extension please don't hesitate to contact the GiNaC
4511 authors---they will happily incorporate them into future versions.
4514 * What does not belong into GiNaC:: What to avoid.
4515 * Symbolic functions:: Implementing symbolic functions.
4516 * Adding classes:: Defining new algebraic classes.
4520 @node What does not belong into GiNaC, Symbolic functions, Extending GiNaC, Extending GiNaC
4521 @c node-name, next, previous, up
4522 @section What doesn't belong into GiNaC
4524 @cindex @command{ginsh}
4525 First of all, GiNaC's name must be read literally. It is designed to be
4526 a library for use within C++. The tiny @command{ginsh} accompanying
4527 GiNaC makes this even more clear: it doesn't even attempt to provide a
4528 language. There are no loops or conditional expressions in
4529 @command{ginsh}, it is merely a window into the library for the
4530 programmer to test stuff (or to show off). Still, the design of a
4531 complete CAS with a language of its own, graphical capabilities and all
4532 this on top of GiNaC is possible and is without doubt a nice project for
4535 There are many built-in functions in GiNaC that do not know how to
4536 evaluate themselves numerically to a precision declared at runtime
4537 (using @code{Digits}). Some may be evaluated at certain points, but not
4538 generally. This ought to be fixed. However, doing numerical
4539 computations with GiNaC's quite abstract classes is doomed to be
4540 inefficient. For this purpose, the underlying foundation classes
4541 provided by CLN are much better suited.
4544 @node Symbolic functions, Adding classes, What does not belong into GiNaC, Extending GiNaC
4545 @c node-name, next, previous, up
4546 @section Symbolic functions
4548 The easiest and most instructive way to start with is probably to
4549 implement your own function. GiNaC's functions are objects of class
4550 @code{function}. The preprocessor is then used to convert the function
4551 names to objects with a corresponding serial number that is used
4552 internally to identify them. You usually need not worry about this
4553 number. New functions may be inserted into the system via a kind of
4554 `registry'. It is your responsibility to care for some functions that
4555 are called when the user invokes certain methods. These are usual
4556 C++-functions accepting a number of @code{ex} as arguments and returning
4557 one @code{ex}. As an example, if we have a look at a simplified
4558 implementation of the cosine trigonometric function, we first need a
4559 function that is called when one wishes to @code{eval} it. It could
4560 look something like this:
4563 static ex cos_eval_method(const ex & x)
4565 // if (!x%(2*Pi)) return 1
4566 // if (!x%Pi) return -1
4567 // if (!x%Pi/2) return 0
4568 // care for other cases...
4569 return cos(x).hold();
4573 @cindex @code{hold()}
4575 The last line returns @code{cos(x)} if we don't know what else to do and
4576 stops a potential recursive evaluation by saying @code{.hold()}, which
4577 sets a flag to the expression signaling that it has been evaluated. We
4578 should also implement a method for numerical evaluation and since we are
4579 lazy we sweep the problem under the rug by calling someone else's
4580 function that does so, in this case the one in class @code{numeric}:
4583 static ex cos_evalf(const ex & x)
4585 if (is_a<numeric>(x))
4586 return cos(ex_to<numeric>(x));
4588 return cos(x).hold();
4592 Differentiation will surely turn up and so we need to tell @code{cos}
4593 what the first derivative is (higher derivatives (@code{.diff(x,3)} for
4594 instance are then handled automatically by @code{basic::diff} and
4598 static ex cos_deriv(const ex & x, unsigned diff_param)
4604 @cindex product rule
4605 The second parameter is obligatory but uninteresting at this point. It
4606 specifies which parameter to differentiate in a partial derivative in
4607 case the function has more than one parameter and its main application
4608 is for correct handling of the chain rule. For Taylor expansion, it is
4609 enough to know how to differentiate. But if the function you want to
4610 implement does have a pole somewhere in the complex plane, you need to
4611 write another method for Laurent expansion around that point.
4613 Now that all the ingredients for @code{cos} have been set up, we need
4614 to tell the system about it. This is done by a macro and we are not
4615 going to describe how it expands, please consult your preprocessor if you
4619 REGISTER_FUNCTION(cos, eval_func(cos_eval).
4620 evalf_func(cos_evalf).
4621 derivative_func(cos_deriv));
4624 The first argument is the function's name used for calling it and for
4625 output. The second binds the corresponding methods as options to this
4626 object. Options are separated by a dot and can be given in an arbitrary
4627 order. GiNaC functions understand several more options which are always
4628 specified as @code{.option(params)}, for example a method for series
4629 expansion @code{.series_func(cos_series)}. Again, if no series
4630 expansion method is given, GiNaC defaults to simple Taylor expansion,
4631 which is correct if there are no poles involved as is the case for the
4632 @code{cos} function. The way GiNaC handles poles in case there are any
4633 is best understood by studying one of the examples, like the Gamma
4634 (@code{tgamma}) function for instance. (In essence the function first
4635 checks if there is a pole at the evaluation point and falls back to
4636 Taylor expansion if there isn't. Then, the pole is regularized by some
4637 suitable transformation.) Also, the new function needs to be declared
4638 somewhere. This may also be done by a convenient preprocessor macro:
4641 DECLARE_FUNCTION_1P(cos)
4644 The suffix @code{_1P} stands for @emph{one parameter}. Of course, this
4645 implementation of @code{cos} is very incomplete and lacks several safety
4646 mechanisms. Please, have a look at the real implementation in GiNaC.
4647 (By the way: in case you are worrying about all the macros above we can
4648 assure you that functions are GiNaC's most macro-intense classes. We
4649 have done our best to avoid macros where we can.)
4652 @node Adding classes, A Comparison With Other CAS, Symbolic functions, Extending GiNaC
4653 @c node-name, next, previous, up
4654 @section Adding classes
4656 If you are doing some very specialized things with GiNaC you may find that
4657 you have to implement your own algebraic classes to fit your needs. This
4658 section will explain how to do this by giving the example of a simple
4659 'string' class. After reading this section you will know how to properly
4660 declare a GiNaC class and what the minimum required member functions are
4661 that you have to implement. We only cover the implementation of a 'leaf'
4662 class here (i.e. one that doesn't contain subexpressions). Creating a
4663 container class like, for example, a class representing tensor products is
4664 more involved but this section should give you enough information so you can
4665 consult the source to GiNaC's predefined classes if you want to implement
4666 something more complicated.
4668 @subsection GiNaC's run-time type information system
4670 @cindex hierarchy of classes
4672 All algebraic classes (that is, all classes that can appear in expressions)
4673 in GiNaC are direct or indirect subclasses of the class @code{basic}. So a
4674 @code{basic *} (which is essentially what an @code{ex} is) represents a
4675 generic pointer to an algebraic class. Occasionally it is necessary to find
4676 out what the class of an object pointed to by a @code{basic *} really is.
4677 Also, for the unarchiving of expressions it must be possible to find the
4678 @code{unarchive()} function of a class given the class name (as a string). A
4679 system that provides this kind of information is called a run-time type
4680 information (RTTI) system. The C++ language provides such a thing (see the
4681 standard header file @file{<typeinfo>}) but for efficiency reasons GiNaC
4682 implements its own, simpler RTTI.
4684 The RTTI in GiNaC is based on two mechanisms:
4689 The @code{basic} class declares a member variable @code{tinfo_key} which
4690 holds an unsigned integer that identifies the object's class. These numbers
4691 are defined in the @file{tinfos.h} header file for the built-in GiNaC
4692 classes. They all start with @code{TINFO_}.
4695 By means of some clever tricks with static members, GiNaC maintains a list
4696 of information for all classes derived from @code{basic}. The information
4697 available includes the class names, the @code{tinfo_key}s, and pointers
4698 to the unarchiving functions. This class registry is defined in the
4699 @file{registrar.h} header file.
4703 The disadvantage of this proprietary RTTI implementation is that there's
4704 a little more to do when implementing new classes (C++'s RTTI works more
4705 or less automatic) but don't worry, most of the work is simplified by
4708 @subsection A minimalistic example
4710 Now we will start implementing a new class @code{mystring} that allows
4711 placing character strings in algebraic expressions (this is not very useful,
4712 but it's just an example). This class will be a direct subclass of
4713 @code{basic}. You can use this sample implementation as a starting point
4714 for your own classes.
4716 The code snippets given here assume that you have included some header files
4722 #include <stdexcept>
4723 using namespace std;
4725 #include <ginac/ginac.h>
4726 using namespace GiNaC;
4729 The first thing we have to do is to define a @code{tinfo_key} for our new
4730 class. This can be any arbitrary unsigned number that is not already taken
4731 by one of the existing classes but it's better to come up with something
4732 that is unlikely to clash with keys that might be added in the future. The
4733 numbers in @file{tinfos.h} are modeled somewhat after the class hierarchy
4734 which is not a requirement but we are going to stick with this scheme:
4737 const unsigned TINFO_mystring = 0x42420001U;
4740 Now we can write down the class declaration. The class stores a C++
4741 @code{string} and the user shall be able to construct a @code{mystring}
4742 object from a C or C++ string:
4745 class mystring : public basic
4747 GINAC_DECLARE_REGISTERED_CLASS(mystring, basic)
4750 mystring(const string &s);
4751 mystring(const char *s);
4757 GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
4760 The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
4761 macros are defined in @file{registrar.h}. They take the name of the class
4762 and its direct superclass as arguments and insert all required declarations
4763 for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
4764 the first line after the opening brace of the class definition. The
4765 @code{GINAC_IMPLEMENT_REGISTERED_CLASS} may appear anywhere else in the
4766 source (at global scope, of course, not inside a function).
4768 @code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
4769 declarations of the default and copy constructor, the destructor, the
4770 assignment operator and a couple of other functions that are required. It
4771 also defines a type @code{inherited} which refers to the superclass so you
4772 don't have to modify your code every time you shuffle around the class
4773 hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} implements the copy
4774 constructor, the destructor and the assignment operator.
4776 Now there are nine member functions we have to implement to get a working
4782 @code{mystring()}, the default constructor.
4785 @code{void destroy(bool call_parent)}, which is used in the destructor and the
4786 assignment operator to free dynamically allocated members. The @code{call_parent}
4787 specifies whether the @code{destroy()} function of the superclass is to be
4791 @code{void copy(const mystring &other)}, which is used in the copy constructor
4792 and assignment operator to copy the member variables over from another
4793 object of the same class.
4796 @code{void archive(archive_node &n)}, the archiving function. This stores all
4797 information needed to reconstruct an object of this class inside an
4798 @code{archive_node}.
4801 @code{mystring(const archive_node &n, const lst &sym_lst)}, the unarchiving
4802 constructor. This constructs an instance of the class from the information
4803 found in an @code{archive_node}.
4806 @code{ex unarchive(const archive_node &n, const lst &sym_lst)}, the static
4807 unarchiving function. It constructs a new instance by calling the unarchiving
4811 @code{int compare_same_type(const basic &other)}, which is used internally
4812 by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
4813 -1, depending on the relative order of this object and the @code{other}
4814 object. If it returns 0, the objects are considered equal.
4815 @strong{Note:} This has nothing to do with the (numeric) ordering
4816 relationship expressed by @code{<}, @code{>=} etc (which cannot be defined
4817 for non-numeric classes). For example, @code{numeric(1).compare_same_type(numeric(2))}
4818 may return +1 even though 1 is clearly smaller than 2. Every GiNaC class
4819 must provide a @code{compare_same_type()} function, even those representing
4820 objects for which no reasonable algebraic ordering relationship can be
4824 And, of course, @code{mystring(const string &s)} and @code{mystring(const char *s)}
4825 which are the two constructors we declared.
4829 Let's proceed step-by-step. The default constructor looks like this:
4832 mystring::mystring() : inherited(TINFO_mystring)
4834 // dynamically allocate resources here if required
4838 The golden rule is that in all constructors you have to set the
4839 @code{tinfo_key} member to the @code{TINFO_*} value of your class. Otherwise
4840 it will be set by the constructor of the superclass and all hell will break
4841 loose in the RTTI. For your convenience, the @code{basic} class provides
4842 a constructor that takes a @code{tinfo_key} value, which we are using here
4843 (remember that in our case @code{inherited = basic}). If the superclass
4844 didn't have such a constructor, we would have to set the @code{tinfo_key}
4845 to the right value manually.
4847 In the default constructor you should set all other member variables to
4848 reasonable default values (we don't need that here since our @code{str}
4849 member gets set to an empty string automatically). The constructor(s) are of
4850 course also the right place to allocate any dynamic resources you require.
4852 Next, the @code{destroy()} function:
4855 void mystring::destroy(bool call_parent)
4857 // free dynamically allocated resources here if required
4859 inherited::destroy(call_parent);
4863 This function is where we free all dynamically allocated resources. We
4864 don't have any so we're not doing anything here, but if we had, for
4865 example, used a C-style @code{char *} to store our string, this would be
4866 the place to @code{delete[]} the string storage. If @code{call_parent}
4867 is true, we have to call the @code{destroy()} function of the superclass
4868 after we're done (to mimic C++'s automatic invocation of superclass
4869 destructors where @code{destroy()} is called from outside a destructor).
4871 The @code{copy()} function just copies over the member variables from
4875 void mystring::copy(const mystring &other)
4877 inherited::copy(other);
4882 We can simply overwrite the member variables here. There's no need to worry
4883 about dynamically allocated storage. The assignment operator (which is
4884 automatically defined by @code{GINAC_IMPLEMENT_REGISTERED_CLASS}, as you
4885 recall) calls @code{destroy()} before it calls @code{copy()}. You have to
4886 explicitly call the @code{copy()} function of the superclass here so
4887 all the member variables will get copied.
4889 Next are the three functions for archiving. You have to implement them even
4890 if you don't plan to use archives, but the minimum required implementation
4891 is really simple. First, the archiving function:
4894 void mystring::archive(archive_node &n) const
4896 inherited::archive(n);
4897 n.add_string("string", str);
4901 The only thing that is really required is calling the @code{archive()}
4902 function of the superclass. Optionally, you can store all information you
4903 deem necessary for representing the object into the passed
4904 @code{archive_node}. We are just storing our string here. For more
4905 information on how the archiving works, consult the @file{archive.h} header
4908 The unarchiving constructor is basically the inverse of the archiving
4912 mystring::mystring(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
4914 n.find_string("string", str);
4918 If you don't need archiving, just leave this function empty (but you must
4919 invoke the unarchiving constructor of the superclass). Note that we don't
4920 have to set the @code{tinfo_key} here because it is done automatically
4921 by the unarchiving constructor of the @code{basic} class.
4923 Finally, the unarchiving function:
4926 ex mystring::unarchive(const archive_node &n, const lst &sym_lst)
4928 return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
4932 You don't have to understand how exactly this works. Just copy these
4933 four lines into your code literally (replacing the class name, of
4934 course). It calls the unarchiving constructor of the class and unless
4935 you are doing something very special (like matching @code{archive_node}s
4936 to global objects) you don't need a different implementation. For those
4937 who are interested: setting the @code{dynallocated} flag puts the object
4938 under the control of GiNaC's garbage collection. It will get deleted
4939 automatically once it is no longer referenced.
4941 Our @code{compare_same_type()} function uses a provided function to compare
4945 int mystring::compare_same_type(const basic &other) const
4947 const mystring &o = static_cast<const mystring &>(other);
4948 int cmpval = str.compare(o.str);
4951 else if (cmpval < 0)
4958 Although this function takes a @code{basic &}, it will always be a reference
4959 to an object of exactly the same class (objects of different classes are not
4960 comparable), so the cast is safe. If this function returns 0, the two objects
4961 are considered equal (in the sense that @math{A-B=0}), so you should compare
4962 all relevant member variables.
4964 Now the only thing missing is our two new constructors:
4967 mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s)
4969 // dynamically allocate resources here if required
4972 mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s)
4974 // dynamically allocate resources here if required
4978 No surprises here. We set the @code{str} member from the argument and
4979 remember to pass the right @code{tinfo_key} to the @code{basic} constructor.
4981 That's it! We now have a minimal working GiNaC class that can store
4982 strings in algebraic expressions. Let's confirm that the RTTI works:
4985 ex e = mystring("Hello, world!");
4986 cout << is_a<mystring>(e) << endl;
4989 cout << e.bp->class_name() << endl;
4993 Obviously it does. Let's see what the expression @code{e} looks like:
4997 // -> [mystring object]
5000 Hm, not exactly what we expect, but of course the @code{mystring} class
5001 doesn't yet know how to print itself. This is done in the @code{print()}
5002 member function. Let's say that we wanted to print the string surrounded
5006 class mystring : public basic
5010 void print(const print_context &c, unsigned level = 0) const;
5014 void mystring::print(const print_context &c, unsigned level) const
5016 // print_context::s is a reference to an ostream
5017 c.s << '\"' << str << '\"';
5021 The @code{level} argument is only required for container classes to
5022 correctly parenthesize the output. Let's try again to print the expression:
5026 // -> "Hello, world!"
5029 Much better. The @code{mystring} class can be used in arbitrary expressions:
5032 e += mystring("GiNaC rulez");
5034 // -> "GiNaC rulez"+"Hello, world!"
5037 (GiNaC's automatic term reordering is in effect here), or even
5040 e = pow(mystring("One string"), 2*sin(Pi-mystring("Another string")));
5042 // -> "One string"^(2*sin(-"Another string"+Pi))
5045 Whether this makes sense is debatable but remember that this is only an
5046 example. At least it allows you to implement your own symbolic algorithms
5049 Note that GiNaC's algebraic rules remain unchanged:
5052 e = mystring("Wow") * mystring("Wow");
5056 e = pow(mystring("First")-mystring("Second"), 2);
5057 cout << e.expand() << endl;
5058 // -> -2*"First"*"Second"+"First"^2+"Second"^2
5061 There's no way to, for example, make GiNaC's @code{add} class perform string
5062 concatenation. You would have to implement this yourself.
5064 @subsection Automatic evaluation
5066 @cindex @code{hold()}
5067 @cindex @code{eval()}
5069 When dealing with objects that are just a little more complicated than the
5070 simple string objects we have implemented, chances are that you will want to
5071 have some automatic simplifications or canonicalizations performed on them.
5072 This is done in the evaluation member function @code{eval()}. Let's say that
5073 we wanted all strings automatically converted to lowercase with
5074 non-alphabetic characters stripped, and empty strings removed:
5077 class mystring : public basic
5081 ex eval(int level = 0) const;
5085 ex mystring::eval(int level) const
5088 for (int i=0; i<str.length(); i++) @{
5090 if (c >= 'A' && c <= 'Z')
5091 new_str += tolower(c);
5092 else if (c >= 'a' && c <= 'z')
5096 if (new_str.length() == 0)
5099 return mystring(new_str).hold();
5103 The @code{level} argument is used to limit the recursion depth of the
5104 evaluation. We don't have any subexpressions in the @code{mystring}
5105 class so we are not concerned with this. If we had, we would call the
5106 @code{eval()} functions of the subexpressions with @code{level - 1} as
5107 the argument if @code{level != 1}. The @code{hold()} member function
5108 sets a flag in the object that prevents further evaluation. Otherwise
5109 we might end up in an endless loop. When you want to return the object
5110 unmodified, use @code{return this->hold();}.
5112 Let's confirm that it works:
5115 ex e = mystring("Hello, world!") + mystring("!?#");
5119 e = mystring("Wow!") + mystring("WOW") + mystring(" W ** o ** W");
5124 @subsection Other member functions
5126 We have implemented only a small set of member functions to make the class
5127 work in the GiNaC framework. For a real algebraic class, there are probably
5128 some more functions that you will want to re-implement, such as
5129 @code{evalf()}, @code{series()} or @code{op()}. Have a look at @file{basic.h}
5130 or the header file of the class you want to make a subclass of to see
5131 what's there. One member function that you will most likely want to
5132 implement for terminal classes like the described string class is
5133 @code{calcchash()} that returns an @code{unsigned} hash value for the object
5134 which will allow GiNaC to compare and canonicalize expressions much more
5137 You can, of course, also add your own new member functions. Remember,
5138 that the RTTI may be used to get information about what kinds of objects
5139 you are dealing with (the position in the class hierarchy) and that you
5140 can always extract the bare object from an @code{ex} by stripping the
5141 @code{ex} off using the @code{ex_to<mystring>(e)} function when that
5142 should become a need.
5144 That's it. May the source be with you!
5147 @node A Comparison With Other CAS, Advantages, Adding classes, Top
5148 @c node-name, next, previous, up
5149 @chapter A Comparison With Other CAS
5152 This chapter will give you some information on how GiNaC compares to
5153 other, traditional Computer Algebra Systems, like @emph{Maple},
5154 @emph{Mathematica} or @emph{Reduce}, where it has advantages and
5155 disadvantages over these systems.
5158 * Advantages:: Strengths of the GiNaC approach.
5159 * Disadvantages:: Weaknesses of the GiNaC approach.
5160 * Why C++?:: Attractiveness of C++.
5163 @node Advantages, Disadvantages, A Comparison With Other CAS, A Comparison With Other CAS
5164 @c node-name, next, previous, up
5167 GiNaC has several advantages over traditional Computer
5168 Algebra Systems, like
5173 familiar language: all common CAS implement their own proprietary
5174 grammar which you have to learn first (and maybe learn again when your
5175 vendor decides to `enhance' it). With GiNaC you can write your program
5176 in common C++, which is standardized.
5180 structured data types: you can build up structured data types using
5181 @code{struct}s or @code{class}es together with STL features instead of
5182 using unnamed lists of lists of lists.
5185 strongly typed: in CAS, you usually have only one kind of variables
5186 which can hold contents of an arbitrary type. This 4GL like feature is
5187 nice for novice programmers, but dangerous.
5190 development tools: powerful development tools exist for C++, like fancy
5191 editors (e.g. with automatic indentation and syntax highlighting),
5192 debuggers, visualization tools, documentation generators@dots{}
5195 modularization: C++ programs can easily be split into modules by
5196 separating interface and implementation.
5199 price: GiNaC is distributed under the GNU Public License which means
5200 that it is free and available with source code. And there are excellent
5201 C++-compilers for free, too.
5204 extendable: you can add your own classes to GiNaC, thus extending it on
5205 a very low level. Compare this to a traditional CAS that you can
5206 usually only extend on a high level by writing in the language defined
5207 by the parser. In particular, it turns out to be almost impossible to
5208 fix bugs in a traditional system.
5211 multiple interfaces: Though real GiNaC programs have to be written in
5212 some editor, then be compiled, linked and executed, there are more ways
5213 to work with the GiNaC engine. Many people want to play with
5214 expressions interactively, as in traditional CASs. Currently, two such
5215 windows into GiNaC have been implemented and many more are possible: the
5216 tiny @command{ginsh} that is part of the distribution exposes GiNaC's
5217 types to a command line and second, as a more consistent approach, an
5218 interactive interface to the Cint C++ interpreter has been put together
5219 (called GiNaC-cint) that allows an interactive scripting interface
5220 consistent with the C++ language. It is available from the usual GiNaC
5224 seamless integration: it is somewhere between difficult and impossible
5225 to call CAS functions from within a program written in C++ or any other
5226 programming language and vice versa. With GiNaC, your symbolic routines
5227 are part of your program. You can easily call third party libraries,
5228 e.g. for numerical evaluation or graphical interaction. All other
5229 approaches are much more cumbersome: they range from simply ignoring the
5230 problem (i.e. @emph{Maple}) to providing a method for `embedding' the
5231 system (i.e. @emph{Yacas}).
5234 efficiency: often large parts of a program do not need symbolic
5235 calculations at all. Why use large integers for loop variables or
5236 arbitrary precision arithmetics where @code{int} and @code{double} are
5237 sufficient? For pure symbolic applications, GiNaC is comparable in
5238 speed with other CAS.
5243 @node Disadvantages, Why C++?, Advantages, A Comparison With Other CAS
5244 @c node-name, next, previous, up
5245 @section Disadvantages
5247 Of course it also has some disadvantages:
5252 advanced features: GiNaC cannot compete with a program like
5253 @emph{Reduce} which exists for more than 30 years now or @emph{Maple}
5254 which grows since 1981 by the work of dozens of programmers, with
5255 respect to mathematical features. Integration, factorization,
5256 non-trivial simplifications, limits etc. are missing in GiNaC (and are
5257 not planned for the near future).
5260 portability: While the GiNaC library itself is designed to avoid any
5261 platform dependent features (it should compile on any ANSI compliant C++
5262 compiler), the currently used version of the CLN library (fast large
5263 integer and arbitrary precision arithmetics) can only by compiled
5264 without hassle on systems with the C++ compiler from the GNU Compiler
5265 Collection (GCC).@footnote{This is because CLN uses PROVIDE/REQUIRE like
5266 macros to let the compiler gather all static initializations, which
5267 works for GNU C++ only. Feel free to contact the authors in case you
5268 really believe that you need to use a different compiler. We have
5269 occasionally used other compilers and may be able to give you advice.}
5270 GiNaC uses recent language features like explicit constructors, mutable
5271 members, RTTI, @code{dynamic_cast}s and STL, so ANSI compliance is meant
5272 literally. Recent GCC versions starting at 2.95.3, although itself not
5273 yet ANSI compliant, support all needed features.
5278 @node Why C++?, Internal Structures, Disadvantages, A Comparison With Other CAS
5279 @c node-name, next, previous, up
5282 Why did we choose to implement GiNaC in C++ instead of Java or any other
5283 language? C++ is not perfect: type checking is not strict (casting is
5284 possible), separation between interface and implementation is not
5285 complete, object oriented design is not enforced. The main reason is
5286 the often scolded feature of operator overloading in C++. While it may
5287 be true that operating on classes with a @code{+} operator is rarely
5288 meaningful, it is perfectly suited for algebraic expressions. Writing
5289 @math{3x+5y} as @code{3*x+5*y} instead of
5290 @code{x.times(3).plus(y.times(5))} looks much more natural.
5291 Furthermore, the main developers are more familiar with C++ than with
5292 any other programming language.
5295 @node Internal Structures, Expressions are reference counted, Why C++? , Top
5296 @c node-name, next, previous, up
5297 @appendix Internal Structures
5300 * Expressions are reference counted::
5301 * Internal representation of products and sums::
5304 @node Expressions are reference counted, Internal representation of products and sums, Internal Structures, Internal Structures
5305 @c node-name, next, previous, up
5306 @appendixsection Expressions are reference counted
5308 @cindex reference counting
5309 @cindex copy-on-write
5310 @cindex garbage collection
5311 An expression is extremely light-weight since internally it works like a
5312 handle to the actual representation and really holds nothing more than a
5313 pointer to some other object. What this means in practice is that
5314 whenever you create two @code{ex} and set the second equal to the first
5315 no copying process is involved. Instead, the copying takes place as soon
5316 as you try to change the second. Consider the simple sequence of code:
5320 #include <ginac/ginac.h>
5321 using namespace std;
5322 using namespace GiNaC;
5326 symbol x("x"), y("y"), z("z");
5329 e1 = sin(x + 2*y) + 3*z + 41;
5330 e2 = e1; // e2 points to same object as e1
5331 cout << e2 << endl; // prints sin(x+2*y)+3*z+41
5332 e2 += 1; // e2 is copied into a new object
5333 cout << e2 << endl; // prints sin(x+2*y)+3*z+42
5337 The line @code{e2 = e1;} creates a second expression pointing to the
5338 object held already by @code{e1}. The time involved for this operation
5339 is therefore constant, no matter how large @code{e1} was. Actual
5340 copying, however, must take place in the line @code{e2 += 1;} because
5341 @code{e1} and @code{e2} are not handles for the same object any more.
5342 This concept is called @dfn{copy-on-write semantics}. It increases
5343 performance considerably whenever one object occurs multiple times and
5344 represents a simple garbage collection scheme because when an @code{ex}
5345 runs out of scope its destructor checks whether other expressions handle
5346 the object it points to too and deletes the object from memory if that
5347 turns out not to be the case. A slightly less trivial example of
5348 differentiation using the chain-rule should make clear how powerful this
5353 symbol x("x"), y("y");
5357 ex e3 = diff(sin(e2), x); // first derivative of sin(e2) by x
5358 cout << e1 << endl // prints x+3*y
5359 << e2 << endl // prints (x+3*y)^3
5360 << e3 << endl; // prints 3*(x+3*y)^2*cos((x+3*y)^3)
5364 Here, @code{e1} will actually be referenced three times while @code{e2}
5365 will be referenced two times. When the power of an expression is built,
5366 that expression needs not be copied. Likewise, since the derivative of
5367 a power of an expression can be easily expressed in terms of that
5368 expression, no copying of @code{e1} is involved when @code{e3} is
5369 constructed. So, when @code{e3} is constructed it will print as
5370 @code{3*(x+3*y)^2*cos((x+3*y)^3)} but the argument of @code{cos()} only
5371 holds a reference to @code{e2} and the factor in front is just
5374 As a user of GiNaC, you cannot see this mechanism of copy-on-write
5375 semantics. When you insert an expression into a second expression, the
5376 result behaves exactly as if the contents of the first expression were
5377 inserted. But it may be useful to remember that this is not what
5378 happens. Knowing this will enable you to write much more efficient
5379 code. If you still have an uncertain feeling with copy-on-write
5380 semantics, we recommend you have a look at the
5381 @uref{http://www.cerfnet.com/~mpcline/c++-faq-lite/, C++-FAQ lite} by
5382 Marshall Cline. Chapter 16 covers this issue and presents an
5383 implementation which is pretty close to the one in GiNaC.
5386 @node Internal representation of products and sums, Package Tools, Expressions are reference counted, Internal Structures
5387 @c node-name, next, previous, up
5388 @appendixsection Internal representation of products and sums
5390 @cindex representation
5393 @cindex @code{power}
5394 Although it should be completely transparent for the user of
5395 GiNaC a short discussion of this topic helps to understand the sources
5396 and also explain performance to a large degree. Consider the
5397 unexpanded symbolic expression
5399 $2d^3 \left( 4a + 5b - 3 \right)$
5402 @math{2*d^3*(4*a+5*b-3)}
5404 which could naively be represented by a tree of linear containers for
5405 addition and multiplication, one container for exponentiation with base
5406 and exponent and some atomic leaves of symbols and numbers in this
5411 @cindex pair-wise representation
5412 However, doing so results in a rather deeply nested tree which will
5413 quickly become inefficient to manipulate. We can improve on this by
5414 representing the sum as a sequence of terms, each one being a pair of a
5415 purely numeric multiplicative coefficient and its rest. In the same
5416 spirit we can store the multiplication as a sequence of terms, each
5417 having a numeric exponent and a possibly complicated base, the tree
5418 becomes much more flat:
5422 The number @code{3} above the symbol @code{d} shows that @code{mul}
5423 objects are treated similarly where the coefficients are interpreted as
5424 @emph{exponents} now. Addition of sums of terms or multiplication of
5425 products with numerical exponents can be coded to be very efficient with
5426 such a pair-wise representation. Internally, this handling is performed
5427 by most CAS in this way. It typically speeds up manipulations by an
5428 order of magnitude. The overall multiplicative factor @code{2} and the
5429 additive term @code{-3} look somewhat out of place in this
5430 representation, however, since they are still carrying a trivial
5431 exponent and multiplicative factor @code{1} respectively. Within GiNaC,
5432 this is avoided by adding a field that carries an overall numeric
5433 coefficient. This results in the realistic picture of internal
5436 $2d^3 \left( 4a + 5b - 3 \right)$:
5439 @math{2*d^3*(4*a+5*b-3)}:
5445 This also allows for a better handling of numeric radicals, since
5446 @code{sqrt(2)} can now be carried along calculations. Now it should be
5447 clear, why both classes @code{add} and @code{mul} are derived from the
5448 same abstract class: the data representation is the same, only the
5449 semantics differs. In the class hierarchy, methods for polynomial
5450 expansion and the like are reimplemented for @code{add} and @code{mul},
5451 but the data structure is inherited from @code{expairseq}.
5454 @node Package Tools, ginac-config, Internal representation of products and sums, Top
5455 @c node-name, next, previous, up
5456 @appendix Package Tools
5458 If you are creating a software package that uses the GiNaC library,
5459 setting the correct command line options for the compiler and linker
5460 can be difficult. GiNaC includes two tools to make this process easier.
5463 * ginac-config:: A shell script to detect compiler and linker flags.
5464 * AM_PATH_GINAC:: Macro for GNU automake.
5468 @node ginac-config, AM_PATH_GINAC, Package Tools, Package Tools
5469 @c node-name, next, previous, up
5470 @section @command{ginac-config}
5471 @cindex ginac-config
5473 @command{ginac-config} is a shell script that you can use to determine
5474 the compiler and linker command line options required to compile and
5475 link a program with the GiNaC library.
5477 @command{ginac-config} takes the following flags:
5481 Prints out the version of GiNaC installed.
5483 Prints '-I' flags pointing to the installed header files.
5485 Prints out the linker flags necessary to link a program against GiNaC.
5486 @item --prefix[=@var{PREFIX}]
5487 If @var{PREFIX} is specified, overrides the configured value of @env{$prefix}.
5488 (And of exec-prefix, unless @code{--exec-prefix} is also specified)
5489 Otherwise, prints out the configured value of @env{$prefix}.
5490 @item --exec-prefix[=@var{PREFIX}]
5491 If @var{PREFIX} is specified, overrides the configured value of @env{$exec_prefix}.
5492 Otherwise, prints out the configured value of @env{$exec_prefix}.
5495 Typically, @command{ginac-config} will be used within a configure
5496 script, as described below. It, however, can also be used directly from
5497 the command line using backquotes to compile a simple program. For
5501 c++ -o simple `ginac-config --cppflags` simple.cpp `ginac-config --libs`
5504 This command line might expand to (for example):
5507 cc -o simple -I/usr/local/include simple.cpp -L/usr/local/lib \
5508 -lginac -lcln -lstdc++
5511 Not only is the form using @command{ginac-config} easier to type, it will
5512 work on any system, no matter how GiNaC was configured.
5515 @node AM_PATH_GINAC, Configure script options, ginac-config, Package Tools
5516 @c node-name, next, previous, up
5517 @section @samp{AM_PATH_GINAC}
5518 @cindex AM_PATH_GINAC
5520 For packages configured using GNU automake, GiNaC also provides
5521 a macro to automate the process of checking for GiNaC.
5524 AM_PATH_GINAC([@var{MINIMUM-VERSION}, [@var{ACTION-IF-FOUND} [, @var{ACTION-IF-NOT-FOUND}]]])
5532 Determines the location of GiNaC using @command{ginac-config}, which is
5533 either found in the user's path, or from the environment variable
5534 @env{GINACLIB_CONFIG}.
5537 Tests the installed libraries to make sure that their version
5538 is later than @var{MINIMUM-VERSION}. (A default version will be used
5542 If the required version was found, sets the @env{GINACLIB_CPPFLAGS} variable
5543 to the output of @command{ginac-config --cppflags} and the @env{GINACLIB_LIBS}
5544 variable to the output of @command{ginac-config --libs}, and calls
5545 @samp{AC_SUBST()} for these variables so they can be used in generated
5546 makefiles, and then executes @var{ACTION-IF-FOUND}.
5549 If the required version was not found, sets @env{GINACLIB_CPPFLAGS} and
5550 @env{GINACLIB_LIBS} to empty strings, and executes @var{ACTION-IF-NOT-FOUND}.
5554 This macro is in file @file{ginac.m4} which is installed in
5555 @file{$datadir/aclocal}. Note that if automake was installed with a
5556 different @samp{--prefix} than GiNaC, you will either have to manually
5557 move @file{ginac.m4} to automake's @file{$datadir/aclocal}, or give
5558 aclocal the @samp{-I} option when running it.
5561 * Configure script options:: Configuring a package that uses AM_PATH_GINAC.
5562 * Example package:: Example of a package using AM_PATH_GINAC.
5566 @node Configure script options, Example package, AM_PATH_GINAC, AM_PATH_GINAC
5567 @c node-name, next, previous, up
5568 @subsection Configuring a package that uses @samp{AM_PATH_GINAC}
5570 Simply make sure that @command{ginac-config} is in your path, and run
5571 the configure script.
5578 The directory where the GiNaC libraries are installed needs
5579 to be found by your system's dynamic linker.
5581 This is generally done by
5584 editing @file{/etc/ld.so.conf} and running @command{ldconfig}
5590 setting the environment variable @env{LD_LIBRARY_PATH},
5593 or, as a last resort,
5596 giving a @samp{-R} or @samp{-rpath} flag (depending on your linker) when
5597 running configure, for instance:
5600 LDFLAGS=-R/home/cbauer/lib ./configure
5605 You can also specify a @command{ginac-config} not in your path by
5606 setting the @env{GINACLIB_CONFIG} environment variable to the
5607 name of the executable
5610 If you move the GiNaC package from its installed location,
5611 you will either need to modify @command{ginac-config} script
5612 manually to point to the new location or rebuild GiNaC.
5623 --with-ginac-prefix=@var{PREFIX}
5624 --with-ginac-exec-prefix=@var{PREFIX}
5627 are provided to override the prefix and exec-prefix that were stored
5628 in the @command{ginac-config} shell script by GiNaC's configure. You are
5629 generally better off configuring GiNaC with the right path to begin with.
5633 @node Example package, Bibliography, Configure script options, AM_PATH_GINAC
5634 @c node-name, next, previous, up
5635 @subsection Example of a package using @samp{AM_PATH_GINAC}
5637 The following shows how to build a simple package using automake
5638 and the @samp{AM_PATH_GINAC} macro. The program used here is @file{simple.cpp}:
5641 #include <ginac/ginac.h>
5645 GiNaC::symbol x("x");
5646 GiNaC::ex a = GiNaC::sin(x);
5647 std::cout << "Derivative of " << a
5648 << " is " << a.diff(x) << std::endl;
5653 You should first read the introductory portions of the automake
5654 Manual, if you are not already familiar with it.
5656 Two files are needed, @file{configure.in}, which is used to build the
5660 dnl Process this file with autoconf to produce a configure script.
5662 AM_INIT_AUTOMAKE(simple.cpp, 1.0.0)
5668 AM_PATH_GINAC(0.9.0, [
5669 LIBS="$LIBS $GINACLIB_LIBS"
5670 CPPFLAGS="$CPPFLAGS $GINACLIB_CPPFLAGS"
5671 ], AC_MSG_ERROR([need to have GiNaC installed]))
5676 The only command in this which is not standard for automake
5677 is the @samp{AM_PATH_GINAC} macro.
5679 That command does the following: If a GiNaC version greater or equal
5680 than 0.7.0 is found, then it adds @env{$GINACLIB_LIBS} to @env{$LIBS}
5681 and @env{$GINACLIB_CPPFLAGS} to @env{$CPPFLAGS}. Otherwise, it dies with
5682 the error message `need to have GiNaC installed'
5684 And the @file{Makefile.am}, which will be used to build the Makefile.
5687 ## Process this file with automake to produce Makefile.in
5688 bin_PROGRAMS = simple
5689 simple_SOURCES = simple.cpp
5692 This @file{Makefile.am}, says that we are building a single executable,
5693 from a single sourcefile @file{simple.cpp}. Since every program
5694 we are building uses GiNaC we simply added the GiNaC options
5695 to @env{$LIBS} and @env{$CPPFLAGS}, but in other circumstances, we might
5696 want to specify them on a per-program basis: for instance by
5700 simple_LDADD = $(GINACLIB_LIBS)
5701 INCLUDES = $(GINACLIB_CPPFLAGS)
5704 to the @file{Makefile.am}.
5706 To try this example out, create a new directory and add the three
5709 Now execute the following commands:
5712 $ automake --add-missing
5717 You now have a package that can be built in the normal fashion
5726 @node Bibliography, Concept Index, Example package, Top
5727 @c node-name, next, previous, up
5728 @appendix Bibliography
5733 @cite{ISO/IEC 14882:1998: Programming Languages: C++}
5736 @cite{CLN: A Class Library for Numbers}, @email{haible@@ilog.fr, Bruno Haible}
5739 @cite{The C++ Programming Language}, Bjarne Stroustrup, 3rd Edition, ISBN 0-201-88954-4, Addison Wesley
5742 @cite{C++ FAQs}, Marshall Cline, ISBN 0-201-58958-3, 1995, Addison Wesley
5745 @cite{Algorithms for Computer Algebra}, Keith O. Geddes, Stephen R. Czapor,
5746 and George Labahn, ISBN 0-7923-9259-0, 1992, Kluwer Academic Publishers, Norwell, Massachusetts
5749 @cite{Computer Algebra: Systems and Algorithms for Algebraic Computation},
5750 James H. Davenport, Yvon Siret, and Evelyne Tournier, ISBN 0-12-204230-1, 1988,
5751 Academic Press, London
5754 @cite{Computer Algebra Systems - A Practical Guide},
5755 Michael J. Wester (editor), ISBN 0-471-98353-5, 1999, Wiley, Chichester
5758 @cite{The Art of Computer Programming, Vol 2: Seminumerical Algorithms},
5759 Donald E. Knuth, ISBN 0-201-89684-2, 1998, Addison Wesley
5762 @cite{The Role of gamma5 in Dimensional Regularization}, Dirk Kreimer, hep-ph/9401354
5767 @node Concept Index, , Bibliography, Top
5768 @c node-name, next, previous, up
5769 @unnumbered Concept Index