This is a tutorial that documents GiNaC @value{VERSION}, an open
framework for symbolic computation within the C++ programming language.
-Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
Permission is granted to make and distribute verbatim copies of
this manual provided the copyright notice and this permission notice
@page
@vskip 0pt plus 1filll
-Copyright @copyright{} 1999-2001 Johannes Gutenberg University Mainz, Germany
+Copyright @copyright{} 1999-2003 Johannes Gutenberg University Mainz, Germany
@sp 2
Permission is granted to make and distribute verbatim copies of
this manual provided the copyright notice and this permission notice
present day computer algebra systems (CAS) are linguistically and
semantically impoverished. Although they are quite powerful tools for
learning math and solving particular problems they lack modern
-linguistical structures that allow for the creation of large-scale
+linguistic structures that allow for the creation of large-scale
projects. GiNaC is an attempt to overcome this situation by extending a
well established and standardized computer language (C++) by some
fundamental symbolic capabilities, thus allowing for integrated systems
@section License
The GiNaC framework for symbolic computation within the C++ programming
-language is Copyright @copyright{} 1999-2001 Johannes Gutenberg
+language is Copyright @copyright{} 1999-2003 Johannes Gutenberg
University Mainz, Germany.
This program is free software; you can redistribute it and/or
pointless) bivariate polynomial with some large coefficients:
@example
+#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
generates Hermite polynomials in a specified free variable.
@example
+#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
[[1,1],[2,-1]]
> A+2*M;
[[1,1],[2,-1]]+2*[[1,3],[-3,2]]
-> evalm(");
+> evalm(%);
[[3,7],[-4,3]]
+> B = [ [0, 0, a], [b, 1, -b], [-1/a, 0, 0] ];
+> evalm(B^(2^12345));
+[[1,0,0],[0,1,0],[0,0,1]]
@end example
Multivariate polynomials and rational functions may be expanded,
> series(tgamma(x),x==0,3);
x^(-1)-Euler+(1/12*Pi^2+1/2*Euler^2)*x+
(-1/3*zeta(3)-1/12*Pi^2*Euler-1/6*Euler^3)*x^2+Order(x^3)
-> evalf(");
+> evalf(%);
x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
-(0.90747907608088628905)*x^2+Order(x^3)
> series(tgamma(2*sin(x)-2),x==Pi/2,6);
-Euler-1/12+Order((x-1/2*Pi)^3)
@end example
-Here we have made use of the @command{ginsh}-command @code{"} to pop the
+Here we have made use of the @command{ginsh}-command @code{%} to pop the
previously evaluated element from @command{ginsh}'s internal stack.
If you ever wanted to convert units in C or C++ and found this is
In order to install GiNaC on your system, some prerequisites need to be
met. First of all, you need to have a C++-compiler adhering to the
-ANSI-standard @cite{ISO/IEC 14882:1998(E)}. We used @acronym{GCC} for
-development so if you have a different compiler you are on your own.
-For the configuration to succeed you need a Posix compliant shell
-installed in @file{/bin/sh}, GNU @command{bash} is fine. Perl is needed
-by the built process as well, since some of the source files are
-automatically generated by Perl scripts. Last but not least, Bruno
-Haible's library @acronym{CLN} is extensively used and needs to be
-installed on your system. Please get it either from
-@uref{ftp://ftp.santafe.edu/pub/gnu/}, from
+ANSI-standard @cite{ISO/IEC 14882:1998(E)}. We used GCC for development
+so if you have a different compiler you are on your own. For the
+configuration to succeed you need a Posix compliant shell installed in
+@file{/bin/sh}, GNU @command{bash} is fine. Perl is needed by the built
+process as well, since some of the source files are automatically
+generated by Perl scripts. Last but not least, Bruno Haible's library
+CLN is extensively used and needs to be installed on your system.
+Please get it either from @uref{ftp://ftp.santafe.edu/pub/gnu/}, from
@uref{ftp://ftpthep.physik.uni-mainz.de/pub/gnu/, GiNaC's FTP site} or
from @uref{ftp://ftp.ilog.fr/pub/Users/haible/gnu/, Bruno Haible's FTP
site} (it is covered by GPL) and install it prior to trying to install
@end itemize
-In addition, you may specify some environment variables.
-@env{CXX} holds the path and the name of the C++ compiler
-in case you want to override the default in your path. (The
-@command{configure} script searches your path for @command{c++},
-@command{g++}, @command{gcc}, @command{CC}, @command{cxx}
-and @command{cc++} in that order.) It may be very useful to
-define some compiler flags with the @env{CXXFLAGS} environment
-variable, like optimization, debugging information and warning
-levels. If omitted, it defaults to @option{-g -O2}.
+In addition, you may specify some environment variables. @env{CXX}
+holds the path and the name of the C++ compiler in case you want to
+override the default in your path. (The @command{configure} script
+searches your path for @command{c++}, @command{g++}, @command{gcc},
+@command{CC}, @command{cxx} and @command{cc++} in that order.) It may
+be very useful to define some compiler flags with the @env{CXXFLAGS}
+environment variable, like optimization, debugging information and
+warning levels. If omitted, it defaults to @option{-g
+-O2}.@footnote{The @command{configure} script is itself generated from
+the file @file{configure.ac}. It is only distributed in packaged
+releases of GiNaC. If you got the naked sources, e.g. from CVS, you
+must generate @command{configure} along with the various
+@file{Makefile.in} by using the @command{autogen.sh} script. This will
+require a fair amount of support from your local toolchain, though.}
The whole process is illustrated in the following two
examples. (Substitute @command{setenv @var{VARIABLE} @var{value}} for
@end example
And here is a configuration for a private static GiNaC library with
-several components sitting in custom places (site-wide @acronym{GCC} and
-private @acronym{CLN}). The compiler is pursuaded to be picky and full
-assertions and debugging information are switched on:
+several components sitting in custom places (site-wide GCC and private
+CLN). The compiler is persuaded to be picky and full assertions and
+debugging information are switched on:
@example
$ export CXX=/usr/local/gnu/bin/c++
$ export CPPFLAGS="$(CPPFLAGS) -I$(HOME)/include"
-$ export CXXFLAGS="$(CXXFLAGS) -DDO_GINAC_ASSERT -ggdb -Wall -ansi -pedantic"
+$ export CXXFLAGS="$(CXXFLAGS) -DDO_GINAC_ASSERT -ggdb -Wall -pedantic"
$ export LDFLAGS="$(LDFLAGS) -L$(HOME)/lib"
$ ./configure --disable-shared --prefix=$(HOME)
@end example
to fiddle around with optimization.
Generally, the top-level Makefile runs recursively to the
-subdirectories. It is therfore safe to go into any subdirectory
+subdirectories. It is therefore safe to go into any subdirectory
(@code{doc/}, @code{ginsh/}, @dots{}) and simply type @code{make}
@var{target} there in case something went wrong.
@menu
* Expressions:: The fundamental GiNaC class.
+* Automatic evaluation:: Evaluation and canonicalization.
+* Error handling:: How the library reports errors.
* The Class Hierarchy:: Overview of GiNaC's classes.
* Symbols:: Symbolic objects.
* Numbers:: Numerical objects.
* Constants:: Pre-defined constants.
-* Fundamental containers:: The power, add and mul classes.
+* Fundamental containers:: Sums, products and powers.
* Lists:: Lists of expressions.
* Mathematical functions:: Mathematical functions.
* Relations:: Equality, Inequality and all that.
@end menu
-@node Expressions, The Class Hierarchy, Basic Concepts, Basic Concepts
+@node Expressions, Automatic evaluation, Basic Concepts, Basic Concepts
@c node-name, next, previous, up
@section Expressions
@cindex expression (class @code{ex})
@code{ex}.
-@node The Class Hierarchy, Symbols, Expressions, Basic Concepts
+@node Automatic evaluation, Error handling, Expressions, Basic Concepts
+@c node-name, next, previous, up
+@section Automatic evaluation and canonicalization of expressions
+@cindex evaluation
+
+GiNaC performs some automatic transformations on expressions, to simplify
+them and put them into a canonical form. Some examples:
+
+@example
+ex MyEx1 = 2*x - 1 + x; // 3*x-1
+ex MyEx2 = x - x; // 0
+ex MyEx3 = cos(2*Pi); // 1
+ex MyEx4 = x*y/x; // y
+@end example
+
+This behavior is usually referred to as @dfn{automatic} or @dfn{anonymous
+evaluation}. GiNaC only performs transformations that are
+
+@itemize @bullet
+@item
+at most of complexity @math{O(n log n)}
+@item
+algebraically correct, possibly except for a set of measure zero (e.g.
+@math{x/x} is transformed to @math{1} although this is incorrect for @math{x=0})
+@end itemize
+
+There are two types of automatic transformations in GiNaC that may not
+behave in an entirely obvious way at first glance:
+
+@itemize
+@item
+The terms of sums and products (and some other things like the arguments of
+symmetric functions, the indices of symmetric tensors etc.) are re-ordered
+into a canonical form that is deterministic, but not lexicographical or in
+any other way easily guessable (it almost always depends on the number and
+order of the symbols you define). However, constructing the same expression
+twice, either implicitly or explicitly, will always result in the same
+canonical form.
+@item
+Expressions of the form 'number times sum' are automatically expanded (this
+has to do with GiNaC's internal representation of sums and products). For
+example
+@example
+ex MyEx5 = 2*(x + y); // 2*x+2*y
+ex MyEx6 = z*(x + y); // z*(x+y)
+@end example
+@end itemize
+
+The general rule is that when you construct expressions, GiNaC automatically
+creates them in canonical form, which might differ from the form you typed in
+your program. This may create some awkward looking output (@samp{-y+x} instead
+of @samp{x-y}) but allows for more efficient operation and usually yields
+some immediate simplifications.
+
+@cindex @code{eval()}
+Internally, the anonymous evaluator in GiNaC is implemented by the methods
+
+@example
+ex ex::eval(int level = 0) const;
+ex basic::eval(int level = 0) const;
+@end example
+
+but unless you are extending GiNaC with your own classes or functions, there
+should never be any reason to call them explicitly. All GiNaC methods that
+transform expressions, like @code{subs()} or @code{normal()}, automatically
+re-evaluate their results.
+
+
+@node Error handling, The Class Hierarchy, Automatic evaluation, Basic Concepts
+@c node-name, next, previous, up
+@section Error handling
+@cindex exceptions
+@cindex @code{pole_error} (class)
+
+GiNaC reports run-time errors by throwing C++ exceptions. All exceptions
+generated by GiNaC are subclassed from the standard @code{exception} class
+defined in the @file{<stdexcept>} header. In addition to the predefined
+@code{logic_error}, @code{domain_error}, @code{out_of_range},
+@code{invalid_argument}, @code{runtime_error}, @code{range_error} and
+@code{overflow_error} types, GiNaC also defines a @code{pole_error}
+exception that gets thrown when trying to evaluate a mathematical function
+at a singularity.
+
+The @code{pole_error} class has a member function
+
+@example
+int pole_error::degree() const;
+@end example
+
+that returns the order of the singularity (or 0 when the pole is
+logarithmic or the order is undefined).
+
+When using GiNaC it is useful to arrange for exceptions to be catched in
+the main program even if you don't want to do any special error handling.
+Otherwise whenever an error occurs in GiNaC, it will be delegated to the
+default exception handler of your C++ compiler's run-time system which
+usually only aborts the program without giving any information what went
+wrong.
+
+Here is an example for a @code{main()} function that catches and prints
+exceptions generated by GiNaC:
+
+@example
+#include <iostream>
+#include <stdexcept>
+#include <ginac/ginac.h>
+using namespace std;
+using namespace GiNaC;
+
+int main()
+@{
+ try @{
+ ...
+ // code using GiNaC
+ ...
+ @} catch (exception &p) @{
+ cerr << p.what() << endl;
+ return 1;
+ @}
+ return 0;
+@}
+@end example
+
+
+@node The Class Hierarchy, Symbols, Error handling, Basic Concepts
@c node-name, next, previous, up
@section The Class Hierarchy
@cindex container
@cindex atom
-To get an idea about what kinds of symbolic composits may be built we
+To get an idea about what kinds of symbolic composites may be built we
have a look at the most important classes in the class hierarchy and
some of the relations among the classes:
@item @code{varidx} @tab Index with variance
@item @code{spinidx} @tab Index with variance and dot (used in Weyl-van-der-Waerden spinor formalism)
@item @code{wildcard} @tab Wildcard for pattern matching
+@item @code{structure} @tab Template for user-defined classes
@end multitable
@end cartouche
+
@node Symbols, Numbers, The Class Hierarchy, Basic Concepts
@c node-name, next, previous, up
@section Symbols
@cindex CLN
@cindex rational
@cindex fraction
-For storing numerical things, GiNaC uses Bruno Haible's library
-@acronym{CLN}. The classes therein serve as foundation classes for
-GiNaC. @acronym{CLN} stands for Class Library for Numbers or
-alternatively for Common Lisp Numbers. In order to find out more about
-@acronym{CLN}'s internals the reader is refered to the documentation of
-that library. @inforef{Introduction, , cln}, for more
-information. Suffice to say that it is by itself build on top of another
-library, the GNU Multiple Precision library @acronym{GMP}, which is an
+For storing numerical things, GiNaC uses Bruno Haible's library CLN.
+The classes therein serve as foundation classes for GiNaC. CLN stands
+for Class Library for Numbers or alternatively for Common Lisp Numbers.
+In order to find out more about CLN's internals, the reader is referred to
+the documentation of that library. @inforef{Introduction, , cln}, for
+more information. Suffice to say that it is by itself build on top of
+another library, the GNU Multiple Precision library GMP, which is an
extremely fast library for arbitrary long integers and rationals as well
as arbitrary precision floating point numbers. It is very commonly used
-by several popular cryptographic applications. @acronym{CLN} extends
-@acronym{GMP} by several useful things: First, it introduces the complex
-number field over either reals (i.e. floating point numbers with
-arbitrary precision) or rationals. Second, it automatically converts
-rationals to integers if the denominator is unity and complex numbers to
-real numbers if the imaginary part vanishes and also correctly treats
-algebraic functions. Third it provides good implementations of
-state-of-the-art algorithms for all trigonometric and hyperbolic
-functions as well as for calculation of some useful constants.
+by several popular cryptographic applications. CLN extends GMP by
+several useful things: First, it introduces the complex number field
+over either reals (i.e. floating point numbers with arbitrary precision)
+or rationals. Second, it automatically converts rationals to integers
+if the denominator is unity and complex numbers to real numbers if the
+imaginary part vanishes and also correctly treats algebraic functions.
+Third it provides good implementations of state-of-the-art algorithms
+for all trigonometric and hyperbolic functions as well as for
+calculation of some useful constants.
The user can construct an object of class @code{numeric} in several
ways. The following example shows the four most important constructors.
integers, construction from C-float and construction from a string:
@example
+#include <iostream>
#include <ginac/ginac.h>
using namespace GiNaC;
numeric trott("1.0841015122311136151E-2");
std::cout << two*p << std::endl; // floating point 6.283...
+ ...
+@end example
+
+@cindex @code{I}
+@cindex complex numbers
+The imaginary unit in GiNaC is a predefined @code{numeric} object with the
+name @code{I}:
+
+@example
+ ...
+ numeric z1 = 2-3*I; // exact complex number 2-3i
+ numeric z2 = 5.9+1.6*I; // complex floating point number
@}
@end example
-It may be tempting to construct numbers writing @code{numeric r(3/2)}.
+It may be tempting to construct fractions by writing @code{numeric r(3/2)}.
This would, however, call C's built-in operator @code{/} for integers
first and result in a numeric holding a plain integer 1. @strong{Never
use the operator @code{/} on integers} unless you know exactly what you
digits:
@example
+#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
@example
in 17 digits:
-0.333333333333333333
-3.14159265358979324
+0.33333333333333333334
+3.1415926535897932385
in 60 digits:
-0.333333333333333333333333333333333333333333333333333333333333333333
-3.14159265358979323846264338327950288419716939937510582097494459231
+0.33333333333333333333333333333333333333333333333333333333333333333334
+3.1415926535897932384626433832795028841971693993751058209749445923078
@end example
+@cindex rounding
+Note that the last number is not necessarily rounded as you would
+naively expect it to be rounded in the decimal system. But note also,
+that in both cases you got a couple of extra digits. This is because
+numbers are internally stored by CLN as chunks of binary digits in order
+to match your machine's word size and to not waste precision. Thus, on
+architectures with different word size, the above output might even
+differ with regard to actually computed digits.
+
It should be clear that objects of class @code{numeric} should be used
for constructing numbers or for doing arithmetic with them. The objects
one deals with most of the time are the polymorphic expressions @code{ex}.
some multiple of its denominator and test what comes out:
@example
+#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
holds a rational number represented as integer numerator and integer
denominator. When multiplied by 10, the denominator becomes unity and
the result is automatically converted to a pure integer again.
-Internally, the underlying @acronym{CLN} is responsible for this
-behaviour and we refer the reader to @acronym{CLN}'s documentation.
-Suffice to say that the same behaviour applies to complex numbers as
-well as return values of certain functions. Complex numbers are
-automatically converted to real numbers if the imaginary part becomes
-zero. The full set of tests that can be applied is listed in the
-following table.
+Internally, the underlying CLN is responsible for this behavior and we
+refer the reader to CLN's documentation. Suffice to say that
+the same behavior applies to complex numbers as well as return values of
+certain functions. Complex numbers are automatically converted to real
+numbers if the imaginary part becomes zero. The full set of tests that
+can be applied is listed in the following table.
@cartouche
@multitable @columnfractions .30 .70
@end multitable
@end cartouche
+@subsection Converting numbers
+
+Sometimes it is desirable to convert a @code{numeric} object back to a
+built-in arithmetic type (@code{int}, @code{double}, etc.). The @code{numeric}
+class provides a couple of methods for this purpose:
+
+@cindex @code{to_int()}
+@cindex @code{to_long()}
+@cindex @code{to_double()}
+@cindex @code{to_cl_N()}
+@example
+int numeric::to_int() const;
+long numeric::to_long() const;
+double numeric::to_double() const;
+cln::cl_N numeric::to_cl_N() const;
+@end example
+
+@code{to_int()} and @code{to_long()} only work when the number they are
+applied on is an exact integer. Otherwise the program will halt with a
+message like @samp{Not a 32-bit integer}. @code{to_double()} applied on a
+rational number will return a floating-point approximation. Both
+@code{to_int()/to_long()} and @code{to_double()} discard the imaginary
+part of complex numbers.
+
@node Constants, Fundamental containers, Numbers, Basic Concepts
@c node-name, next, previous, up
@node Fundamental containers, Lists, Constants, Basic Concepts
@c node-name, next, previous, up
-@section Fundamental containers: the @code{power}, @code{add} and @code{mul} classes
+@section Sums, products and powers
@cindex polynomial
@cindex @code{add}
@cindex @code{mul}
@cindex @code{power}
-Simple polynomial expressions are written down in GiNaC pretty much like
+Simple rational expressions are written down in GiNaC pretty much like
in other CAS or like expressions involving numerical variables in C.
The necessary operators @code{+}, @code{-}, @code{*} and @code{/} have
been overloaded to achieve this goal. When you run the following
Also, expressions involving integer exponents are very frequently used,
which makes it even more dangerous to overload @code{^} since it is then
hard to distinguish between the semantics as exponentiation and the one
-for exclusive or. (It would be embarassing to return @code{1} where one
+for exclusive or. (It would be embarrassing to return @code{1} where one
has requested @code{2^3}.)
@end itemize
and safe simplifications are carried out like transforming
@code{3*x+4-x} to @code{2*x+4}.
-The general rule is that when you construct such objects, GiNaC
-automatically creates them in canonical form, which might differ from
-the form you typed in your program. This allows for rapid comparison of
-expressions, since after all @code{a-a} is simply zero. Note, that the
-canonical form is not necessarily lexicographical ordering or in any way
-easily guessable. It is only guaranteed that constructing the same
-expression twice, either implicitly or explicitly, results in the same
-canonical form.
-
@node Lists, Mathematical functions, Fundamental containers, Basic Concepts
@c node-name, next, previous, up
@cindex @code{prepend()}
@cindex @code{remove_first()}
@cindex @code{remove_last()}
+@cindex @code{remove_all()}
The GiNaC class @code{lst} serves for holding a @dfn{list} of arbitrary
-expressions. These are sometimes used to supply a variable number of
-arguments of the same type to GiNaC methods such as @code{subs()} and
-@code{to_rational()}, so you should have a basic understanding about them.
+expressions. They are not as ubiquitous as in many other computer algebra
+packages, but are sometimes used to supply a variable number of arguments of
+the same type to GiNaC methods such as @code{subs()} and @code{to_rational()},
+so you should have a basic understanding of them.
Lists of up to 16 expressions can be directly constructed from single
expressions:
symbol x("x"), y("y");
lst l(x, 2, y, x+y);
// now, l is a list holding the expressions 'x', '2', 'y', and 'x+y'
- // ...
+ ...
@end example
Use the @code{nops()} method to determine the size (number of expressions) of
-a list and the @code{op()} method to access individual elements:
+a list and the @code{op()} method or the @code{[]} operator to access
+individual elements:
@example
- // ...
- cout << l.nops() << endl; // prints '4'
- cout << l.op(2) << " " << l.op(0) << endl; // prints 'y x'
- // ...
+ ...
+ cout << l.nops() << endl; // prints '4'
+ cout << l.op(2) << " " << l[0] << endl; // prints 'y x'
+ ...
+@end example
+
+As with the standard @code{list<T>} container, accessing random elements of a
+@code{lst} is generally an operation of order @math{O(N)}. Faster read-only
+sequential access to the elements of a list is possible with the
+iterator types provided by the @code{lst} class:
+
+@example
+typedef ... lst::const_iterator;
+typedef ... lst::const_reverse_iterator;
+lst::const_iterator lst::begin() const;
+lst::const_iterator lst::end() const;
+lst::const_reverse_iterator lst::rbegin() const;
+lst::const_reverse_iterator lst::rend() const;
+@end example
+
+For example, to print the elements of a list individually you can use:
+
+@example
+ ...
+ // O(N)
+ for (lst::const_iterator i = l.begin(); i != l.end(); ++i)
+ cout << *i << endl;
+ ...
+@end example
+
+which is one order faster than
+
+@example
+ ...
+ // O(N^2)
+ for (size_t i = 0; i < l.nops(); ++i)
+ cout << l.op(i) << endl;
+ ...
+@end example
+
+These iterators also allow you to use some of the algorithms provided by
+the C++ standard library:
+
+@example
+ ...
+ // print the elements of the list (requires #include <iterator>)
+ copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));
+
+ // sum up the elements of the list (requires #include <numeric>)
+ ex sum = accumulate(l.begin(), l.end(), ex(0));
+ cout << sum << endl; // prints '2+2*x+2*y'
+ ...
+@end example
+
+@code{lst} is one of the few GiNaC classes that allow in-place modifications
+(the only other one is @code{matrix}). You can modify single elements:
+
+@example
+ ...
+ l[1] = 42; // l is now @{x, 42, y, x+y@}
+ l.let_op(1) = 7; // l is now @{x, 7, y, x+y@}
+ ...
@end example
You can append or prepend an expression to a list with the @code{append()}
and @code{prepend()} methods:
@example
- // ...
- l.append(4*x); // l is now @{x, 2, y, x+y, 4*x@}
- l.prepend(0); // l is now @{0, x, 2, y, x+y, 4*x@}
- // ...
+ ...
+ l.append(4*x); // l is now @{x, 7, y, x+y, 4*x@}
+ l.prepend(0); // l is now @{0, x, 7, y, x+y, 4*x@}
+ ...
@end example
-Finally you can remove the first or last element of a list with
-@code{remove_first()} and @code{remove_last()}:
+You can remove the first or last element of a list with @code{remove_first()}
+and @code{remove_last()}:
@example
- // ...
- l.remove_first(); // l is now @{x, 2, y, x+y, 4*x@}
- l.remove_last(); // l is now @{x, 2, y, x+y@}
+ ...
+ l.remove_first(); // l is now @{x, 7, y, x+y, 4*x@}
+ l.remove_last(); // l is now @{x, 7, y, x+y@}
+ ...
+@end example
+
+You can remove all the elements of a list with @code{remove_all()}:
+
+@example
+ ...
+ l.remove_all(); // l is now empty
+ ...
+@end example
+
+You can bring the elements of a list into a canonical order with @code{sort()}:
+
+@example
+ ...
+ lst l1(x, 2, y, x+y);
+ lst l2(2, x+y, x, y);
+ l1.sort();
+ l2.sort();
+ // l1 and l2 are now equal
+ ...
+@end example
+
+Finally, you can remove all but the first element of consecutive groups of
+elements with @code{unique()}:
+
+@example
+ ...
+ lst l3(x, 2, 2, 2, y, x+y, y+x);
+ l3.unique(); // l3 is now @{x, 2, y, x+y@}
@}
@end example
functions, where the argument list is templated. This means that
whenever you call @code{GiNaC::sin(1)} it is equivalent to
@code{sin(ex(1))} and will therefore not result in a floating point
-numeber. Unless of course the function prototype is explicitly
+number. Unless of course the function prototype is explicitly
overridden -- which is the case for arguments of type @code{numeric}
(not wrapped inside an @code{ex}). Hence, in order to obtain a floating
point number of class @code{numeric} you should call
There are a couple of ways to construct matrices, with or without preset
elements:
+@cindex @code{lst_to_matrix()}
+@cindex @code{diag_matrix()}
+@cindex @code{unit_matrix()}
+@cindex @code{symbolic_matrix()}
@example
matrix::matrix(unsigned r, unsigned c);
matrix::matrix(unsigned r, unsigned c, const lst & l);
ex lst_to_matrix(const lst & l);
ex diag_matrix(const lst & l);
+ex unit_matrix(unsigned x);
+ex unit_matrix(unsigned r, unsigned c);
+ex symbolic_matrix(unsigned r, unsigned c, const string & base_name);
+ex symbolic_matrix(unsigned r, unsigned c, const string & base_name, const string & tex_base_name);
@end example
The first two functions are @code{matrix} constructors which create a matrix
with @samp{r} rows and @samp{c} columns. The matrix elements can be
initialized from a (flat) list of expressions @samp{l}. Otherwise they are
all set to zero. The @code{lst_to_matrix()} function constructs a matrix
-from a list of lists, each list representing a matrix row. Finally,
-@code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
-elements. Note that the last two functions return expressions, not matrix
-objects.
+from a list of lists, each list representing a matrix row. @code{diag_matrix()}
+constructs a diagonal matrix given the list of diagonal elements.
+@code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r} by @samp{c})
+unit matrix. And finally, @code{symbolic_matrix} constructs a matrix filled
+with newly generated symbols made of the specified base name and the
+position of each element in the matrix.
Matrix elements can be accessed and set using the parenthesis (function call)
operator:
the @code{op()} method. But C++-style subscripting with square brackets
@samp{[]} is not available.
-Here are a couple of examples that all construct the same 2x2 diagonal
-matrix:
+Here are a couple of examples of constructing matrices:
@example
@{
symbol a("a"), b("b");
- ex e;
matrix M(2, 2);
M(0, 0) = a;
M(1, 1) = b;
- e = M;
-
- e = matrix(2, 2, lst(a, 0, 0, b));
+ cout << M << endl;
+ // -> [[a,0],[0,b]]
- e = lst_to_matrix(lst(lst(a, 0), lst(0, b)));
+ cout << matrix(2, 2, lst(a, 0, 0, b)) << endl;
+ // -> [[a,0],[0,b]]
- e = diag_matrix(lst(a, b));
+ cout << lst_to_matrix(lst(lst(a, 0), lst(0, b))) << endl;
+ // -> [[a,0],[0,b]]
- cout << e << endl;
+ cout << diag_matrix(lst(a, b)) << endl;
// -> [[a,0],[0,b]]
+
+ cout << unit_matrix(3) << endl;
+ // -> [[1,0,0],[0,1,0],[0,0,1]]
+
+ cout << symbolic_matrix(2, 3, "x") << endl;
+ // -> [[x00,x01,x02],[x10,x11,x12]]
@}
@end example
@cindex @code{transpose()}
-@cindex @code{inverse()}
There are three ways to do arithmetic with matrices. The first (and most
-efficient one) is to use the methods provided by the @code{matrix} class:
+direct one) is to use the methods provided by the @code{matrix} class:
@example
matrix matrix::add(const matrix & other) const;
matrix matrix::mul(const matrix & other) const;
matrix matrix::mul_scalar(const ex & other) const;
matrix matrix::pow(const ex & expn) const;
-matrix matrix::transpose(void) const;
-matrix matrix::inverse(void) const;
+matrix matrix::transpose() const;
@end example
All of these methods return the result as a new matrix object. Here is an
The @code{matrix} class provides a couple of additional methods for
computing determinants, traces, and characteristic polynomials:
+@cindex @code{determinant()}
+@cindex @code{trace()}
+@cindex @code{charpoly()}
+@example
+ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
+ex matrix::trace() const;
+ex matrix::charpoly(const ex & lambda) const;
+@end example
+
+The @samp{algo} argument of @code{determinant()} allows to select
+between different algorithms for calculating the determinant. The
+asymptotic speed (as parametrized by the matrix size) can greatly differ
+between those algorithms, depending on the nature of the matrix'
+entries. The possible values are defined in the @file{flags.h} header
+file. By default, GiNaC uses a heuristic to automatically select an
+algorithm that is likely (but not guaranteed) to give the result most
+quickly.
+
+@cindex @code{inverse()}
+@cindex @code{solve()}
+Matrices may also be inverted using the @code{ex matrix::inverse()}
+method and linear systems may be solved with:
+
@example
-ex matrix::determinant(unsigned algo = determinant_algo::automatic) const;
-ex matrix::trace(void) const;
-ex matrix::charpoly(const symbol & lambda) const;
+matrix matrix::solve(const matrix & vars, const matrix & rhs, unsigned algo=solve_algo::automatic) const;
@end example
-The @samp{algo} argument of @code{determinant()} allows to select between
-different algorithms for calculating the determinant. The possible values
-are defined in the @file{flags.h} header file. By default, GiNaC uses a
-heuristic to automatically select an algorithm that is likely to give the
-result most quickly.
+Assuming the matrix object this method is applied on is an @code{m}
+times @code{n} matrix, then @code{vars} must be a @code{n} times
+@code{p} matrix of symbolic indeterminates and @code{rhs} a @code{m}
+times @code{p} matrix. The returned matrix then has dimension @code{n}
+times @code{p} and in the case of an underdetermined system will still
+contain some of the indeterminates from @code{vars}. If the system is
+overdetermined, an exception is thrown.
@node Indexed objects, Non-commutative objects, Matrices, Basic Concepts
A simple example shall illustrate the concepts:
@example
+#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
symbol A("A");
cout << indexed(A, i, j) << endl;
// -> A.i.j
+ cout << index_dimensions << indexed(A, i, j) << endl;
+ // -> A.i[3].j[3]
+ cout << dflt; // reset cout to default output format (dimensions hidden)
...
@end example
the symbol @code{A} as its base expression and the two indices @code{i} and
@code{j}.
+The dimensions of indices are normally not visible in the output, but one
+can request them to be printed with the @code{index_dimensions} manipulator,
+as shown above.
+
Note the difference between the indices @code{i} and @code{j} which are of
-class @code{idx}, and the index values which are the sybols @code{i_sym}
+class @code{idx}, and the index values which are the symbols @code{i_sym}
and @code{j_sym}. The indices of indexed objects cannot directly be symbols
or numbers but must be index objects. For example, the following is not
correct and will raise an exception:
The methods
@example
-ex idx::get_value(void);
-ex idx::get_dimension(void);
+ex idx::get_value();
+ex idx::get_dimension();
@end example
return the value and dimension of an @code{idx} object. If you have an index
There are also the methods
@example
-bool idx::is_numeric(void);
-bool idx::is_symbolic(void);
-bool idx::is_dim_numeric(void);
-bool idx::is_dim_symbolic(void);
+bool idx::is_numeric();
+bool idx::is_symbolic();
+bool idx::is_dim_numeric();
+bool idx::is_dim_symbolic();
@end example
for checking whether the value and dimension are numeric or symbolic
constructor. The two methods
@example
-bool varidx::is_covariant(void);
-bool varidx::is_contravariant(void);
+bool varidx::is_covariant();
+bool varidx::is_contravariant();
@end example
allow you to check the variance of a @code{varidx} object (use @code{ex_to<varidx>()}
method
@example
-ex varidx::toggle_variance(void);
+ex varidx::toggle_variance();
@end example
which makes a new index with the same value and dimension but the opposite
methods
@example
-bool spinidx::is_dotted(void);
-bool spinidx::is_undotted(void);
+bool spinidx::is_dotted();
+bool spinidx::is_undotted();
@end example
allow you to check whether or not a @code{spinidx} object is dotted (use
Finally, the two methods
@example
-ex spinidx::toggle_dot(void);
-ex spinidx::toggle_variance_dot(void);
+ex spinidx::toggle_dot();
+ex spinidx::toggle_variance_dot();
@end example
create a new index with the same value and dimension but opposite dottedness
// -> 2*A.j.i
cout << indexed(B, sy_anti(), i, j)
+ indexed(B, sy_anti(), j, i) << endl;
- // -> -B.j.i
+ // -> 0
cout << indexed(B, sy_anti(), i, j, k)
- + indexed(B, sy_anti(), j, i, k) << endl;
+ - indexed(B, sy_anti(), j, k, i) << endl;
// -> 0
...
@end example
@cindex @code{get_free_indices()}
-@cindex Dummy index
+@cindex dummy index
@subsection Dummy indices
GiNaC treats certain symbolic index pairs as @dfn{dummy indices} meaning
dummy nor free indices.
To be recognized as a dummy index pair, the two indices must be of the same
-class and dimension and their value must be the same single symbol (an index
-like @samp{2*n+1} is never a dummy index). If the indices are of class
+class and their value must be the same single symbol (an index like
+@samp{2*n+1} is never a dummy index). If the indices are of class
@code{varidx} they must also be of opposite variance; if they are of class
@code{spinidx} they must be both dotted or both undotted.
there is the method
@example
-ex ex::simplify_indexed(void);
+ex ex::simplify_indexed();
ex ex::simplify_indexed(const scalar_products & sp);
@end example
@itemize
@item it checks the consistency of free indices in sums in the same way
@code{get_free_indices()} does
-@item it tries to give dumy indices that appear in different terms of a sum
+@item it tries to give dummy indices that appear in different terms of a sum
the same name to allow simplifications like @math{a_i*b_i-a_j*b_j=0}
@item it (symbolically) calculates all possible dummy index summations/contractions
with the predefined tensors (this will be explained in more detail in the
next section)
+@item it detects contractions that vanish for symmetry reasons, for example
+ the contraction of a symmetric and a totally antisymmetric tensor
@item as a special case of dummy index summation, it can replace scalar products
of two tensors with a user-defined value
@end itemize
The epsilon tensor is totally antisymmetric, its number of indices is equal
to the dimension of the index space (the indices must all be of the same
numeric dimension), and @samp{eps.1.2.3...} (resp. @samp{eps~0~1~2...}) is
-defined to be 1. Its behaviour with indices that have a variance also
+defined to be 1. Its behavior with indices that have a variance also
depends on the signature of the metric. Epsilon tensors are output as
@samp{eps}.
dimensions, the last function creates an epsilon tensor in a 4-dimensional
Minkowski space (the last @code{bool} argument specifies whether the metric
has negative or positive signature, as in the case of the Minkowski metric
-tensor).
+tensor):
+
+@example
+@{
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4),
+ sig(symbol("sig"), 4), lam(symbol("lam"), 4), bet(symbol("bet"), 4);
+ e = lorentz_eps(mu, nu, rho, sig) *
+ lorentz_eps(mu.toggle_variance(), nu.toggle_variance(), lam, bet);
+ cout << simplify_indexed(e) << endl;
+ // -> 2*eta~bet~rho*eta~sig~lam-2*eta~sig~bet*eta~rho~lam
+
+ idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
+ symbol A("A"), B("B");
+ e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(B, k);
+ cout << simplify_indexed(e) << endl;
+ // -> -B.k*A.j*eps.i.k.j
+ e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(A, k);
+ cout << simplify_indexed(e) << endl;
+ // -> 0
+@}
+@end example
@subsection Linear algebra
@end itemize
The @code{clifford} and @code{color} classes are subclasses of
-@code{indexed} because the elements of these algebras ususally carry
+@code{indexed} because the elements of these algebras usually carry
indices. The @code{matrix} class is described in more detail in
@ref{Matrices}.
obtained with the two member functions
@example
-unsigned ex::return_type(void) const;
-unsigned ex::return_type_tinfo(void) const;
+unsigned ex::return_type() const;
+unsigned ex::return_type_tinfo() const;
@end example
The @code{return_type()} function returns one of three values (defined in
ex dirac_ONE(unsigned char rl = 0);
@end example
+@strong{Note:} You must always use @code{dirac_ONE()} when referring to
+multiples of the unity element, even though it's customary to omit it.
+E.g. instead of @code{dirac_gamma(mu)*(dirac_slash(q,4)+m)} you have to
+write @code{dirac_gamma(mu)*(dirac_slash(q,4)+m*dirac_ONE())}. Otherwise,
+GiNaC will complain and/or produce incorrect results.
+
@cindex @code{dirac_gamma5()}
-and there's a special element @samp{gamma5} that commutes with all other
-gammas and in 4 dimensions equals @samp{gamma~0 gamma~1 gamma~2 gamma~3},
-provided by
+There is a special element @samp{gamma5} that commutes with all other
+gammas, has a unit square, and in 4 dimensions equals
+@samp{gamma~0 gamma~1 gamma~2 gamma~3}, provided by
@example
ex dirac_gamma5(unsigned char rl = 0);
@end example
-@cindex @code{dirac_gamma6()}
-@cindex @code{dirac_gamma7()}
-The two additional functions
+@cindex @code{dirac_gammaL()}
+@cindex @code{dirac_gammaR()}
+The chiral projectors @samp{(1+/-gamma5)/2} are also available as proper
+objects, constructed by
@example
-ex dirac_gamma6(unsigned char rl = 0);
-ex dirac_gamma7(unsigned char rl = 0);
+ex dirac_gammaL(unsigned char rl = 0);
+ex dirac_gammaR(unsigned char rl = 0);
@end example
-return @code{dirac_ONE(rl) + dirac_gamma5(rl)} and @code{dirac_ONE(rl) - dirac_gamma5(rl)},
-respectively.
+They observe the relations @samp{gammaL^2 = gammaL}, @samp{gammaR^2 = gammaR},
+and @samp{gammaL gammaR = gammaR gammaL = 0}.
@cindex @code{dirac_slash()}
Finally, the function
ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
@end example
-creates a term of the form @samp{e.mu gamma~mu} with a new and unique index
-whose dimension is given by the @code{dim} argument.
+creates a term that represents a contraction of @samp{e} with the Dirac
+Lorentz vector (it behaves like a term of the form @samp{e.mu gamma~mu}
+with a unique index whose dimension is given by the @code{dim} argument).
+Such slashed expressions are printed with a trailing backslash, e.g. @samp{e\}.
In products of dirac gammas, superfluous unity elements are automatically
-removed, squares are replaced by their values and @samp{gamma5} is
-anticommuted to the front. The @code{simplify_indexed()} function performs
-contractions in gamma strings, for example
+removed, squares are replaced by their values, and @samp{gamma5}, @samp{gammaL}
+and @samp{gammaR} are moved to the front.
+
+The @code{simplify_indexed()} function performs contractions in gamma strings,
+for example
@example
@{
ex e = dirac_gamma(mu) * dirac_slash(a, D)
* dirac_gamma(mu.toggle_variance());
cout << e << endl;
- // -> (gamma~mu*gamma~symbol10*gamma.mu)*a.symbol10
+ // -> gamma~mu*a\*gamma.mu
e = e.simplify_indexed();
cout << e << endl;
- // -> -gamma~symbol10*a.symbol10*D+2*gamma~symbol10*a.symbol10
+ // -> -D*a\+2*a\
cout << e.subs(D == 4) << endl;
- // -> -2*gamma~symbol10*a.symbol10
- // [ == -2 * dirac_slash(a, D) ]
+ // -> -2*a\
...
@}
@end example
ex color_ONE(unsigned char rl = 0);
@end example
+@strong{Note:} You must always use @code{color_ONE()} when referring to
+multiples of the unity element, even though it's customary to omit it.
+E.g. instead of @code{color_T(a)*(color_T(b)*indexed(X,b)+1)} you have to
+write @code{color_T(a)*(color_T(b)*indexed(X,b)+color_ONE())}. Otherwise,
+GiNaC may produce incorrect results.
+
@cindex @code{color_d()}
@cindex @code{color_f()}
-and the functions
+The functions
@example
ex color_d(const ex & a, const ex & b, const ex & c);
@menu
* Information About Expressions::
+* Numerical Evaluation::
* Substituting Expressions::
* Pattern Matching and Advanced Substitutions::
* Applying a Function on Subexpressions::
+* Visitors and Tree Traversal::
* Polynomial Arithmetic:: Working with polynomials.
* Rational Expressions:: Working with rational functions.
* Symbolic Differentiation::
* Series Expansion:: Taylor and Laurent expansion.
* Symmetrization::
* Built-in Functions:: List of predefined mathematical functions.
+* Solving Linear Systems of Equations::
* Input/Output:: Input and output of expressions.
@end menu
-@node Information About Expressions, Substituting Expressions, Methods and Functions, Methods and Functions
+@node Information About Expressions, Numerical Evaluation, Methods and Functions, Methods and Functions
@c node-name, next, previous, up
@section Getting information about expressions
bool is_a<T>(const ex & e);
bool is_exactly_a<T>(const ex & e);
bool ex::info(unsigned flag);
-unsigned ex::return_type(void) const;
-unsigned ex::return_type_tinfo(void) const;
+unsigned ex::return_type() const;
+unsigned ex::return_type_tinfo() const;
@end example
When the test made by @code{is_a<T>()} returns true, it is safe to call
GiNaC provides the two methods
@example
-unsigned ex::nops();
-ex ex::op(unsigned i);
+size_t ex::nops();
+ex ex::op(size_t i);
@end example
for accessing the subexpressions in the container-like GiNaC classes like
Actually, if you construct an expression like @code{a == b}, this will be
represented by an object of the @code{relational} class (@pxref{Relations})
-which is not evaluated until (explicitly or implicitely) cast to a @code{bool}.
+which is not evaluated until (explicitly or implicitly) cast to a @code{bool}.
There are also two methods
for checking whether one expression is equal to another, or equal to zero,
respectively.
-@strong{Warning:} You will also find an @code{ex::compare()} method in the
-GiNaC header files. This method is however only to be used internally by
-GiNaC to establish a canonical sort order for terms, and using it to compare
-expressions will give very surprising results.
+
+@subsection Ordering expressions
+@cindex @code{ex_is_less} (class)
+@cindex @code{ex_is_equal} (class)
+@cindex @code{compare()}
+
+Sometimes it is necessary to establish a mathematically well-defined ordering
+on a set of arbitrary expressions, for example to use expressions as keys
+in a @code{std::map<>} container, or to bring a vector of expressions into
+a canonical order (which is done internally by GiNaC for sums and products).
+
+The operators @code{<}, @code{>} etc. described in the last section cannot
+be used for this, as they don't implement an ordering relation in the
+mathematical sense. In particular, they are not guaranteed to be
+antisymmetric: if @samp{a} and @samp{b} are different expressions, and
+@code{a < b} yields @code{false}, then @code{b < a} doesn't necessarily
+yield @code{true}.
+
+By default, STL classes and algorithms use the @code{<} and @code{==}
+operators to compare objects, which are unsuitable for expressions, but GiNaC
+provides two functors that can be supplied as proper binary comparison
+predicates to the STL:
+
+@example
+class ex_is_less : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+
+class ex_is_equal : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+@end example
+
+For example, to define a @code{map} that maps expressions to strings you
+have to use
+
+@example
+std::map<ex, std::string, ex_is_less> myMap;
+@end example
+
+Omitting the @code{ex_is_less} template parameter will introduce spurious
+bugs because the map operates improperly.
+
+Other examples for the use of the functors:
+
+@example
+std::vector<ex> v;
+// fill vector
+...
+
+// sort vector
+std::sort(v.begin(), v.end(), ex_is_less());
+
+// count the number of expressions equal to '1'
+unsigned num_ones = std::count_if(v.begin(), v.end(),
+ std::bind2nd(ex_is_equal(), 1));
+@end example
+
+The implementation of @code{ex_is_less} uses the member function
+
+@example
+int ex::compare(const ex & other) const;
+@end example
+
+which returns @math{0} if @code{*this} and @code{other} are equal, @math{-1}
+if @code{*this} sorts before @code{other}, and @math{1} if @code{*this} sorts
+after @code{other}.
+
+
+@node Numerical Evaluation, Substituting Expressions, Information About Expressions, Methods and Functions
+@c node-name, next, previous, up
+@section Numercial Evaluation
+@cindex @code{evalf()}
+
+GiNaC keeps algebraic expressions, numbers and constants in their exact form.
+To evaluate them using floating-point arithmetic you need to call
+
+@example
+ex ex::evalf(int level = 0) const;
+@end example
+
+@cindex @code{Digits}
+The accuracy of the evaluation is controlled by the global object @code{Digits}
+which can be assigned an integer value. The default value of @code{Digits}
+is 17. @xref{Numbers}, for more information and examples.
+
+To evaluate an expression to a @code{double} floating-point number you can
+call @code{evalf()} followed by @code{numeric::to_double()}, like this:
+
+@example
+@{
+ // Approximate sin(x/Pi)
+ symbol x("x");
+ ex e = series(sin(x/Pi), x == 0, 6);
+
+ // Evaluate numerically at x=0.1
+ ex f = evalf(e.subs(x == 0.1));
+
+ // ex_to<numeric> is an unsafe cast, so check the type first
+ if (is_a<numeric>(f)) @{
+ double d = ex_to<numeric>(f).to_double();
+ cout << d << endl;
+ // -> 0.0318256
+ @} else
+ // error
+@}
+@end example
-@node Substituting Expressions, Pattern Matching and Advanced Substitutions, Information About Expressions, Methods and Functions
+@node Substituting Expressions, Pattern Matching and Advanced Substitutions, Numerical Evaluation, Methods and Functions
@c node-name, next, previous, up
@section Substituting expressions
@cindex @code{subs()}
expressions via the @code{.subs()} method:
@example
-ex ex::subs(const ex & e);
-ex ex::subs(const lst & syms, const lst & repls);
+ex ex::subs(const ex & e, unsigned options = 0);
+ex ex::subs(const lst & syms, const lst & repls, unsigned options = 0);
@end example
In the first form, @code{subs()} accepts a relational of the form
contain the same number of elements). Using this form, you would write
@code{subs(lst(x, y), lst(y, x))} to exchange @samp{x} and @samp{y}.
+The optional last argument to @code{subs()} is a combination of
+@code{subs_options} flags. There are two options available:
+@code{subs_options::no_pattern} disables pattern matching, which makes
+large @code{subs()} operations significantly faster if you are not using
+patterns. The second option, @code{subs_options::algebraic} enables
+algebraic substitutions in products and powers.
+@ref{Pattern Matching and Advanced Substitutions}, for more information
+about patterns and algebraic substitutions.
+
@code{subs()} performs syntactic substitution of any complete algebraic
object; it does not try to match sub-expressions as is demonstrated by the
following example:
are also valid patterns.
@end itemize
+@subsection Matching expressions
@cindex @code{match()}
The most basic application of patterns is to check whether an expression
matches a given pattern. This is done by the function
In general, having more than one single wildcard as a term of a sum or a
factor of a product (such as @samp{a+$0+$1}) will lead to unpredictable or
-amgiguous results.
+ambiguous results.
Here are some examples in @command{ginsh} to demonstrate how it works (the
@code{match()} function in @command{ginsh} returns @samp{FAIL} if the
@{$0==x^2@}
@end example
+@subsection Matching parts of expressions
@cindex @code{has()}
A more general way to look for patterns in expressions is provided by the
member function
(Note the absence of "x".)
> expand((sin(x)+sin(y))*(a+b));
sin(y)*a+sin(x)*b+sin(x)*a+sin(y)*b
-> find(",sin($1));
+> find(%,sin($1));
@{sin(y),sin(x)@}
@end example
+@subsection Substituting expressions
@cindex @code{subs()}
Probably the most useful application of patterns is to use them for
substituting expressions with the @code{subs()} method. Wildcards can be
b^3+a^3+(x+y)^3
> subs(a^4+b^4+(x+y)^4,$1^2==$1^3);
b^4+a^4+(x+y)^4
-> subs((a+b+c)^2,a+b=x);
+> subs((a+b+c)^2,a+b==x);
(a+b+c)^2
> subs((a+b+c)^2,a+b+$1==x+$1);
(x+c)^2
-> subs(a+2*b,a+b=x);
+> subs(a+2*b,a+b==x);
a+2*b
> subs(4*x^3-2*x^2+5*x-1,x==a);
-1+5*a-2*a^2+4*a^3
@}
@end example
+@subsection Algebraic substitutions
+Supplying the @code{subs_options::algebraic} option to @code{subs()}
+enables smarter, algebraic substitutions in products and powers. If you want
+to substitute some factors of a product, you only need to list these factors
+in your pattern. Furthermore, if an (integer) power of some expression occurs
+in your pattern and in the expression that you want the substitution to occur
+in, it can be substituted as many times as possible, without getting negative
+powers.
+
+An example clarifies it all (hopefully):
+
+@example
+cout << (a*a*a*a+b*b*b*b+pow(x+y,4)).subs(wild()*wild()==pow(wild(),3),
+ subs_options::algebraic) << endl;
+// --> (y+x)^6+b^6+a^6
+
+cout << ((a+b+c)*(a+b+c)).subs(a+b==x,subs_options::algebraic) << endl;
+// --> (c+b+a)^2
+// Powers and products are smart, but addition is just the same.
+
+cout << ((a+b+c)*(a+b+c)).subs(a+b+wild()==x+wild(), subs_options::algebraic)
+ << endl;
+// --> (x+c)^2
+// As I said: addition is just the same.
+
+cout << (pow(a,5)*pow(b,7)+2*b).subs(b*b*a==x,subs_options::algebraic) << endl;
+// --> x^3*b*a^2+2*b
+
+cout << (pow(a,-5)*pow(b,-7)+2*b).subs(1/(b*b*a)==x,subs_options::algebraic)
+ << endl;
+// --> 2*b+x^3*b^(-1)*a^(-2)
+
+cout << (4*x*x*x-2*x*x+5*x-1).subs(x==a,subs_options::algebraic) << endl;
+// --> -1-2*a^2+4*a^3+5*a
+
+cout << (4*x*x*x-2*x*x+5*x-1).subs(pow(x,wild())==pow(a,wild()),
+ subs_options::algebraic) << endl;
+// --> -1+5*x+4*x^3-2*x^2
+// You should not really need this kind of patterns very often now.
+// But perhaps this it's-not-a-bug-it's-a-feature (c/sh)ould still change.
+
+cout << ex(sin(1+sin(x))).subs(sin(wild())==cos(wild()),
+ subs_options::algebraic) << endl;
+// --> cos(1+cos(x))
+
+cout << expand((a*sin(x+y)*sin(x+y)+a*cos(x+y)*cos(x+y)+b)
+ .subs((pow(cos(wild()),2)==1-pow(sin(wild()),2)),
+ subs_options::algebraic)) << endl;
+// --> b+a
+@end example
+
-@node Applying a Function on Subexpressions, Polynomial Arithmetic, Pattern Matching and Advanced Substitutions, Methods and Functions
+@node Applying a Function on Subexpressions, Visitors and Tree Traversal, Pattern Matching and Advanced Substitutions, Methods and Functions
@c node-name, next, previous, up
@section Applying a Function on Subexpressions
-@cindex Tree traversal
+@cindex tree traversal
@cindex @code{map()}
Sometimes you may want to perform an operation on specific parts of an
return ex_to<matrix>(e).trace();
else if (is_a<add>(e)) @{
ex sum = 0;
- for (unsigned i=0; i<e.nops(); i++)
+ for (size_t i=0; i<e.nops(); i++)
sum += calc_trace(e.op(i));
return sum;
@} else if (is_a<mul>)(e)) @{
operations:
@example
-static ex ex::map(map_function & f) const;
-static ex ex::map(ex (*f)(const ex & e)) const;
+ex ex::map(map_function & f) const;
+ex ex::map(ex (*f)(const ex & e)) const;
@end example
In the first (preferred) form, @code{map()} takes a function object that
@}
@end example
+Here is another example for you to meditate over. It removes quadratic
+terms in a variable from an expanded polynomial:
+
+@example
+struct map_rem_quad : public map_function @{
+ ex var;
+ map_rem_quad(const ex & var_) : var(var_) @{@}
+
+ ex operator()(const ex & e)
+ @{
+ if (is_a<add>(e) || is_a<mul>(e))
+ return e.map(*this);
+ else if (is_a<power>(e) &&
+ e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
+ return 0;
+ else
+ return e;
+ @}
+@};
+
+...
+
+@{
+ symbol x("x"), y("y");
+
+ ex e;
+ for (int i=0; i<8; i++)
+ e += pow(x, i) * pow(y, 8-i) * (i+1);
+ cout << e << endl;
+ // -> 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
+
+ map_rem_quad rem_quad(x);
+ cout << rem_quad(e) << endl;
+ // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
+@}
+@end example
+
@command{ginsh} offers a slightly different implementation of @code{map()}
that allows applying algebraic functions to operands. The second argument
to @code{map()} is an expression containing the wildcard @samp{$0} which
@end example
-@node Polynomial Arithmetic, Rational Expressions, Applying a Function on Subexpressions, Methods and Functions
+@node Visitors and Tree Traversal, Polynomial Arithmetic, Applying a Function on Subexpressions, Methods and Functions
@c node-name, next, previous, up
-@section Polynomial arithmetic
+@section Visitors and Tree Traversal
+@cindex tree traversal
+@cindex @code{visitor} (class)
+@cindex @code{accept()}
+@cindex @code{visit()}
+@cindex @code{traverse()}
+@cindex @code{traverse_preorder()}
+@cindex @code{traverse_postorder()}
-@subsection Expanding and collecting
-@cindex @code{expand()}
-@cindex @code{collect()}
+Suppose that you need a function that returns a list of all indices appearing
+in an arbitrary expression. The indices can have any dimension, and for
+indices with variance you always want the covariant version returned.
-A polynomial in one or more variables has many equivalent
-representations. Some useful ones serve a specific purpose. Consider
-for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
-21*y*z + 4*z^2} (written down here in output-style). It is equivalent
-to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
+You can't use @code{get_free_indices()} because you also want to include
+dummy indices in the list, and you can't use @code{find()} as it needs
+specific index dimensions (and it would require two passes: one for indices
+with variance, one for plain ones).
+
+The obvious solution to this problem is a tree traversal with a type switch,
+such as the following:
+
+@example
+void gather_indices_helper(const ex & e, lst & l)
+@{
+ if (is_a<varidx>(e)) @{
+ const varidx & vi = ex_to<varidx>(e);
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @} else if (is_a<idx>(e)) @{
+ l.append(e);
+ @} else @{
+ size_t n = e.nops();
+ for (size_t i = 0; i < n; ++i)
+ gather_indices_helper(e.op(i), l);
+ @}
+@}
+
+lst gather_indices(const ex & e)
+@{
+ lst l;
+ gather_indices_helper(e, l);
+ l.sort();
+ l.unique();
+ return l;
+@}
+@end example
+
+This works fine but fans of object-oriented programming will feel
+uncomfortable with the type switch. One reason is that there is a possibility
+for subtle bugs regarding derived classes. If we had, for example, written
+
+@example
+ if (is_a<idx>(e)) @{
+ ...
+ @} else if (is_a<varidx>(e)) @{
+ ...
+@end example
+
+in @code{gather_indices_helper}, the code wouldn't have worked because the
+first line "absorbs" all classes derived from @code{idx}, including
+@code{varidx}, so the special case for @code{varidx} would never have been
+executed.
+
+Also, for a large number of classes, a type switch like the above can get
+unwieldy and inefficient (it's a linear search, after all).
+@code{gather_indices_helper} only checks for two classes, but if you had to
+write a function that required a different implementation for nearly
+every GiNaC class, the result would be very hard to maintain and extend.
+
+The cleanest approach to the problem would be to add a new virtual function
+to GiNaC's class hierarchy. In our example, there would be specializations
+for @code{idx} and @code{varidx} while the default implementation in
+@code{basic} performed the tree traversal. Unfortunately, in C++ it's
+impossible to add virtual member functions to existing classes without
+changing their source and recompiling everything. GiNaC comes with source,
+so you could actually do this, but for a small algorithm like the one
+presented this would be impractical.
+
+One solution to this dilemma is the @dfn{Visitor} design pattern,
+which is implemented in GiNaC (actually, Robert Martin's Acyclic Visitor
+variation, described in detail in
+@uref{http://objectmentor.com/publications/acv.pdf}). Instead of adding
+virtual functions to the class hierarchy to implement operations, GiNaC
+provides a single "bouncing" method @code{accept()} that takes an instance
+of a special @code{visitor} class and redirects execution to the one
+@code{visit()} virtual function of the visitor that matches the type of
+object that @code{accept()} was being invoked on.
+
+Visitors in GiNaC must derive from the global @code{visitor} class as well
+as from the class @code{T::visitor} of each class @code{T} they want to
+visit, and implement the member functions @code{void visit(const T &)} for
+each class.
+
+A call of
+
+@example
+void ex::accept(visitor & v) const;
+@end example
+
+will then dispatch to the correct @code{visit()} member function of the
+specified visitor @code{v} for the type of GiNaC object at the root of the
+expression tree (e.g. a @code{symbol}, an @code{idx} or a @code{mul}).
+
+Here is an example of a visitor:
+
+@example
+class my_visitor
+ : public visitor, // this is required
+ public add::visitor, // visit add objects
+ public numeric::visitor, // visit numeric objects
+ public basic::visitor // visit basic objects
+@{
+ void visit(const add & x)
+ @{ cout << "called with an add object" << endl; @}
+
+ void visit(const numeric & x)
+ @{ cout << "called with a numeric object" << endl; @}
+
+ void visit(const basic & x)
+ @{ cout << "called with a basic object" << endl; @}
+@};
+@end example
+
+which can be used as follows:
+
+@example
+...
+ symbol x("x");
+ ex e1 = 42;
+ ex e2 = 4*x-3;
+ ex e3 = 8*x;
+
+ my_visitor v;
+ e1.accept(v);
+ // prints "called with a numeric object"
+ e2.accept(v);
+ // prints "called with an add object"
+ e3.accept(v);
+ // prints "called with a basic object"
+...
+@end example
+
+The @code{visit(const basic &)} method gets called for all objects that are
+not @code{numeric} or @code{add} and acts as an (optional) default.
+
+From a conceptual point of view, the @code{visit()} methods of the visitor
+behave like a newly added virtual function of the visited hierarchy.
+In addition, visitors can store state in member variables, and they can
+be extended by deriving a new visitor from an existing one, thus building
+hierarchies of visitors.
+
+We can now rewrite our index example from above with a visitor:
+
+@example
+class gather_indices_visitor
+ : public visitor, public idx::visitor, public varidx::visitor
+@{
+ lst l;
+
+ void visit(const idx & i)
+ @{
+ l.append(i);
+ @}
+
+ void visit(const varidx & vi)
+ @{
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @}
+
+public:
+ const lst & get_result() // utility function
+ @{
+ l.sort();
+ l.unique();
+ return l;
+ @}
+@};
+@end example
+
+What's missing is the tree traversal. We could implement it in
+@code{visit(const basic &)}, but GiNaC has predefined methods for this:
+
+@example
+void ex::traverse_preorder(visitor & v) const;
+void ex::traverse_postorder(visitor & v) const;
+void ex::traverse(visitor & v) const;
+@end example
+
+@code{traverse_preorder()} visits a node @emph{before} visiting its
+subexpressions, while @code{traverse_postorder()} visits a node @emph{after}
+visiting its subexpressions. @code{traverse()} is a synonym for
+@code{traverse_preorder()}.
+
+Here is a new implementation of @code{gather_indices()} that uses the visitor
+and @code{traverse()}:
+
+@example
+lst gather_indices(const ex & e)
+@{
+ gather_indices_visitor v;
+ e.traverse(v);
+ return v.get_result();
+@}
+@end example
+
+
+@node Polynomial Arithmetic, Rational Expressions, Visitors and Tree Traversal, Methods and Functions
+@c node-name, next, previous, up
+@section Polynomial arithmetic
+
+@subsection Expanding and collecting
+@cindex @code{expand()}
+@cindex @code{collect()}
+@cindex @code{collect_common_factors()}
+
+A polynomial in one or more variables has many equivalent
+representations. Some useful ones serve a specific purpose. Consider
+for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
+21*y*z + 4*z^2} (written down here in output-style). It is equivalent
+to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
representations are the recursive ones where one collects for exponents
in one of the three variable. Since the factors are themselves
polynomials in the remaining two variables the procedure can be
-repeated. In our expample, two possibilities would be @math{(4*y + z)*x
+repeated. In our example, two possibilities would be @math{(4*y + z)*x
+ 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
x*z}.
To bring an expression into expanded form, its method
@example
-ex ex::expand();
+ex ex::expand(unsigned options = 0);
@end example
may be called. In our example above, this corresponds to @math{4*x*y +
(1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
@end example
+Polynomials can often be brought into a more compact form by collecting
+common factors from the terms of sums. This is accomplished by the function
+
+@example
+ex collect_common_factors(const ex & e);
+@end example
+
+This function doesn't perform a full factorization but only looks for
+factors which are already explicitly present:
+
+@example
+> collect_common_factors(a*x+a*y);
+(x+y)*a
+> collect_common_factors(a*x^2+2*a*x*y+a*y^2);
+a*(2*x*y+y^2+x^2)
+> collect_common_factors(a*(b*(a+c)*x+b*((a+c)*x+(a+c)*y)*y));
+(c+a)*a*(x*y+y^2+x)*b
+@end example
+
@subsection Degree and coefficients
@cindex @code{degree()}
@cindex @code{ldegree()}
polynomial is analyzed:
@example
-#include <ginac/ginac.h>
-using namespace std;
-using namespace GiNaC;
-
-int main()
@{
symbol x("x"), y("y");
ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
The two functions
@example
-ex quo(const ex & a, const ex & b, const symbol & x);
-ex rem(const ex & a, const ex & b, const symbol & x);
+ex quo(const ex & a, const ex & b, const ex & x);
+ex rem(const ex & a, const ex & b, const ex & x);
@end example
compute the quotient and remainder of univariate polynomials in the variable
The additional function
@example
-ex prem(const ex & a, const ex & b, const symbol & x);
+ex prem(const ex & a, const ex & b, const ex & x);
@end example
computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
The methods
@example
-ex ex::unit(const symbol & x);
-ex ex::content(const symbol & x);
-ex ex::primpart(const symbol & x);
+ex ex::unit(const ex & x);
+ex ex::content(const ex & x);
+ex ex::primpart(const ex & x);
@end example
return the unit part, content part, and primitive polynomial of a multivariate
factorization is, however, easily implemented by noting that factors
appearing in a polynomial with power two or more also appear in the
derivative and hence can easily be found by computing the GCD of the
-original polynomial and its derivatives. Any system has an interface
-for this so called square-free factorization. So we provide one, too:
+original polynomial and its derivatives. Any decent system has an
+interface for this so called square-free factorization. So we provide
+one, too:
@example
ex sqrfree(const ex & a, const lst & l = lst());
@end example
-Here is an example that by the way illustrates how the result may depend
-on the order of differentiation:
+Here is an example that by the way illustrates how the exact form of the
+result may slightly depend on the order of differentiation, calling for
+some care with subsequent processing of the result:
@example
...
symbol x("x"), y("y");
- ex BiVarPol = expand(pow(x-2*y*x,3) * pow(x+y,2) * (x-y));
+ ex BiVarPol = expand(pow(2-2*y,3) * pow(1+x*y,2) * pow(x-2*y,2) * (x+y));
cout << sqrfree(BiVarPol, lst(x,y)) << endl;
- // -> (y+x)^2*(-1+6*y+8*y^3-12*y^2)*(y-x)*x^3
+ // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)
cout << sqrfree(BiVarPol, lst(y,x)) << endl;
- // -> (1-2*y)^3*(y+x)^2*(-y+x)*x^3
+ // -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x)
cout << sqrfree(BiVarPol) << endl;
// -> depending on luck, any of the above
...
@end example
+Note also, how factors with the same exponents are not fully factorized
+with this method.
@node Rational Expressions, Symbolic Differentiation, Polynomial Arithmetic, Methods and Functions
@code{.to_rational()}, described below.
This means that both expressions @code{t1} and @code{t2} are indeed
-simplified in this little program:
+simplified in this little code snippet:
@example
-#include <ginac/ginac.h>
-using namespace GiNaC;
-
-int main()
@{
symbol x("x");
ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
faster than using @code{numer()} and @code{denom()} separately.
-@subsection Converting to a rational expression
+@subsection Converting to a polynomial or rational expression
+@cindex @code{to_polynomial()}
@cindex @code{to_rational()}
Some of the methods described so far only work on polynomials or rational
general expressions by using the temporary replacement algorithm described
above. You do this by calling
+@example
+ex ex::to_polynomial(lst &l);
+@end example
+or
@example
ex ex::to_rational(lst &l);
@end example
with the generated temporary symbols and their replacement expressions in
a format that can be used directly for the @code{subs()} method. It can also
already contain a list of replacements from an earlier application of
-@code{.to_rational()}, so it's possible to use it on multiple expressions
-and get consistent results.
+@code{.to_polynomial()} or @code{.to_rational()}, so it's possible to use
+it on multiple expressions and get consistent results.
+
+The difference betwerrn @code{.to_polynomial()} and @code{.to_rational()}
+is probably best illustrated with an example:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex a = 2*x/sin(x) - y/(3*sin(x));
+ cout << a << endl;
+
+ lst lp;
+ ex p = a.to_polynomial(lp);
+ cout << " = " << p << "\n with " << lp << endl;
+ // = symbol3*symbol2*y+2*symbol2*x
+ // with @{symbol2==sin(x)^(-1),symbol3==-1/3@}
+
+ lst lr;
+ ex r = a.to_rational(lr);
+ cout << " = " << r << "\n with " << lr << endl;
+ // = -1/3*symbol4^(-1)*y+2*symbol4^(-1)*x
+ // with @{symbol4==sin(x)@}
+@}
+@end example
-For example,
+The following more useful example will print @samp{sin(x)-cos(x)}:
@example
@{
ex b = sin(x) + cos(x);
ex q;
lst l;
- divide(a.to_rational(l), b.to_rational(l), q);
+ divide(a.to_polynomial(l), b.to_polynomial(l), q);
cout << q.subs(l) << endl;
@}
@end example
-will print @samp{sin(x)-cos(x)}.
-
@node Symbolic Differentiation, Series Expansion, Rational Expressions, Methods and Functions
@c node-name, next, previous, up
the derivatives of all the monomials:
@example
-#include <ginac/ginac.h>
-using namespace GiNaC;
-
-int main()
@{
symbol x("x"), y("y"), z("z");
ex P = pow(x, 5) + pow(x, 2) + y;
- cout << P.diff(x,2) << endl; // 20*x^3 + 2
+ cout << P.diff(x,2) << endl;
+ // -> 20*x^3 + 2
cout << P.diff(y) << endl; // 1
+ // -> 1
cout << P.diff(z) << endl; // 0
+ // -> 0
@}
@end example
@cindex Taylor expansion
@cindex Laurent expansion
@cindex @code{pseries} (class)
+@cindex @code{Order()}
Expressions know how to expand themselves as a Taylor series or (more
generally) a Laurent series. As in most conventional Computer Algebra
@math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
series raised to the power @math{-2}.
-@cindex M@'echain's formula
+@cindex Machin's formula
As another instructive application, let us calculate the numerical
value of Archimedes' constant
@tex
$\pi$
@end tex
(for which there already exists the built-in constant @code{Pi})
-using M@'echain's amazing formula
+using John Machin's amazing formula
@tex
$\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
@end tex
@ifnottex
@math{Pi==16*atan(1/5)-4*atan(1/239)}.
@end ifnottex
-We may expand the arcus tangent around @code{0} and insert the fractions
-@code{1/5} and @code{1/239}. But, as we have seen, a series in GiNaC
-carries an order term with it and the question arises what the system is
-supposed to do when the fractions are plugged into that order term. The
-solution is to use the function @code{series_to_poly()} to simply strip
-the order term off:
+This equation (and similar ones) were used for over 200 years for
+computing digits of pi (see @cite{Pi Unleashed}). We may expand the
+arcus tangent around @code{0} and insert the fractions @code{1/5} and
+@code{1/239}. However, as we have seen, a series in GiNaC carries an
+order term with it and the question arises what the system is supposed
+to do when the fractions are plugged into that order term. The solution
+is to use the function @code{series_to_poly()} to simply strip the order
+term off:
@example
#include <ginac/ginac.h>
using namespace GiNaC;
-ex mechain_pi(int degr)
+ex machin_pi(int degr)
@{
symbol x;
ex pi_expansion = series_to_poly(atan(x).series(x,degr));
using std::endl; // ...dealing with this namespace std.
ex pi_frac;
for (int i=2; i<12; i+=2) @{
- pi_frac = mechain_pi(i);
+ pi_frac = machin_pi(i);
cout << i << ":\t" << pi_frac << endl
<< "\t" << pi_frac.evalf() << endl;
@}
@end example
-@node Built-in Functions, Input/Output, Symmetrization, Methods and Functions
+@node Built-in Functions, Solving Linear Systems of Equations, Symmetrization, Methods and Functions
@c node-name, next, previous, up
@section Predefined mathematical functions
@item @strong{Name} @tab @strong{Function}
@item @code{abs(x)}
@tab absolute value
+@cindex @code{abs()}
@item @code{csgn(x)}
@tab complex sign
+@cindex @code{csgn()}
@item @code{sqrt(x)}
-@tab square root (not a GiNaC function proper but equivalent to @code{pow(x, numeric(1, 2)})
+@tab square root (not a GiNaC function, rather an alias for @code{pow(x, numeric(1, 2))})
+@cindex @code{sqrt()}
@item @code{sin(x)}
@tab sine
+@cindex @code{sin()}
@item @code{cos(x)}
@tab cosine
+@cindex @code{cos()}
@item @code{tan(x)}
@tab tangent
+@cindex @code{tan()}
@item @code{asin(x)}
@tab inverse sine
+@cindex @code{asin()}
@item @code{acos(x)}
@tab inverse cosine
+@cindex @code{acos()}
@item @code{atan(x)}
@tab inverse tangent
+@cindex @code{atan()}
@item @code{atan2(y, x)}
@tab inverse tangent with two arguments
@item @code{sinh(x)}
@tab hyperbolic sine
+@cindex @code{sinh()}
@item @code{cosh(x)}
@tab hyperbolic cosine
+@cindex @code{cosh()}
@item @code{tanh(x)}
@tab hyperbolic tangent
+@cindex @code{tanh()}
@item @code{asinh(x)}
@tab inverse hyperbolic sine
+@cindex @code{asinh()}
@item @code{acosh(x)}
@tab inverse hyperbolic cosine
+@cindex @code{acosh()}
@item @code{atanh(x)}
@tab inverse hyperbolic tangent
+@cindex @code{atanh()}
@item @code{exp(x)}
@tab exponential function
+@cindex @code{exp()}
@item @code{log(x)}
@tab natural logarithm
+@cindex @code{log()}
@item @code{Li2(x)}
@tab Dilogarithm
+@cindex @code{Li2()}
@item @code{zeta(x)}
@tab Riemann's zeta function
+@cindex @code{zeta()}
@item @code{zeta(n, x)}
@tab derivatives of Riemann's zeta function
@item @code{tgamma(x)}
@tab Gamma function
+@cindex @code{tgamma()}
+@cindex Gamma function
@item @code{lgamma(x)}
@tab logarithm of Gamma function
+@cindex @code{lgamma()}
@item @code{beta(x, y)}
@tab Beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
+@cindex @code{beta()}
@item @code{psi(x)}
@tab psi (digamma) function
+@cindex @code{psi()}
@item @code{psi(n, x)}
@tab derivatives of psi function (polygamma functions)
@item @code{factorial(n)}
@tab factorial function
+@cindex @code{factorial()}
@item @code{binomial(n, m)}
@tab binomial coefficients
+@cindex @code{binomial()}
@item @code{Order(x)}
@tab order term function in truncated power series
-@item @code{Derivative(x, l)}
-@tab inert partial differentiation operator (used internally)
+@cindex @code{Order()}
+@item @code{Li(n,x)}
+@tab polylogarithm
+@cindex @code{Li()}
+@item @code{S(n,p,x)}
+@tab Nielsen's generalized polylogarithm
+@cindex @code{S()}
+@item @code{H(m_lst,x)}
+@tab harmonic polylogarithm
+@cindex @code{H()}
+@item @code{Li(m_lst,x_lst)}
+@tab multiple polylogarithm
+@cindex @code{Li()}
+@item @code{mZeta(m_lst)}
+@tab multiple zeta value
+@cindex @code{mZeta()}
@end multitable
@end cartouche
compatible with C99.
-@node Input/Output, Extending GiNaC, Built-in Functions, Methods and Functions
+@node Solving Linear Systems of Equations, Input/Output, Built-in Functions, Methods and Functions
+@c node-name, next, previous, up
+@section Solving Linear Systems of Equations
+@cindex @code{lsolve()}
+
+The function @code{lsolve()} provides a convenient wrapper around some
+matrix operations that comes in handy when a system of linear equations
+needs to be solved:
+
+@example
+ex lsolve(const ex &eqns, const ex &symbols, unsigned options=solve_algo::automatic);
+@end example
+
+Here, @code{eqns} is a @code{lst} of equalities (i.e. class
+@code{relational}) while @code{symbols} is a @code{lst} of
+indeterminates. (@xref{The Class Hierarchy}, for an exposition of class
+@code{lst}).
+
+It returns the @code{lst} of solutions as an expression. As an example,
+let us solve the two equations @code{a*x+b*y==3} and @code{x-y==b}:
+
+@example
+@{
+ symbol a("a"), b("b"), x("x"), y("y");
+ lst eqns;
+ eqns.append(a*x+b*y==3).append(x-y==b);
+ lst vars;
+ vars.append(x).append(y);
+ cout << lsolve(eqns, vars) << endl;
+ // -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
+@end example
+
+When the linear equations @code{eqns} are underdetermined, the solution
+will contain one or more tautological entries like @code{x==x},
+depending on the rank of the system. When they are overdetermined, the
+solution will be an empty @code{lst}. Note the third optional parameter
+to @code{lsolve()}: it accepts the same parameters as
+@code{matrix::solve()}. This is because @code{lsolve} is just a wrapper
+around that method.
+
+
+@node Input/Output, Extending GiNaC, Solving Linear Systems of Equations, Methods and Functions
@c node-name, next, previous, up
@section Input and output of expressions
@cindex I/O
@cindex printing
@cindex output of expressions
-The easiest way to print an expression is to write it to a stream:
+Expressions can simply be written to any stream:
@example
@{
symbol x("x");
- ex e = 4.5+pow(x,2)*3/2;
- cout << e << endl; // prints '(4.5)+3/2*x^2'
+ ex e = 4.5*I+pow(x,2)*3/2;
+ cout << e << endl; // prints '4.5*I+3/2*x^2'
// ...
@end example
-The output format is identical to the @command{ginsh} input syntax and
+The default output format is identical to the @command{ginsh} input syntax and
to that used by most computer algebra systems, but not directly pastable
into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
is printed as @samp{x^2}).
It is possible to print expressions in a number of different formats with
-the method
+a set of stream manipulators;
@example
-void ex::print(const print_context & c, unsigned level = 0);
+std::ostream & dflt(std::ostream & os);
+std::ostream & latex(std::ostream & os);
+std::ostream & tree(std::ostream & os);
+std::ostream & csrc(std::ostream & os);
+std::ostream & csrc_float(std::ostream & os);
+std::ostream & csrc_double(std::ostream & os);
+std::ostream & csrc_cl_N(std::ostream & os);
+std::ostream & index_dimensions(std::ostream & os);
+std::ostream & no_index_dimensions(std::ostream & os);
@end example
-@cindex @code{print_context} (class)
-The type of @code{print_context} object passed in determines the format
-of the output. The possible types are defined in @file{ginac/print.h}.
-All constructors of @code{print_context} and derived classes take an
-@code{ostream &} as their first argument.
+The @code{tree}, @code{latex} and @code{csrc} formats are also available in
+@command{ginsh} via the @code{print()}, @code{print_latex()} and
+@code{print_csrc()} functions, respectively.
-To print an expression in a way that can be directly used in a C or C++
-program, you pass a @code{print_csrc} object like this:
+@cindex @code{dflt}
+All manipulators affect the stream state permanently. To reset the output
+format to the default, use the @code{dflt} manipulator:
@example
// ...
- cout << "float f = ";
- e.print(print_csrc_float(cout));
- cout << ";\n";
+ cout << latex; // all output to cout will be in LaTeX format from now on
+ cout << e << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
+ cout << sin(x/2) << endl; // prints '\sin(\frac@{1@}@{2@} x)'
+ cout << dflt; // revert to default output format
+ cout << e << endl; // prints '4.5*I+3/2*x^2'
+ // ...
+@end example
- cout << "double d = ";
- e.print(print_csrc_double(cout));
- cout << ";\n";
+If you don't want to affect the format of the stream you're working with,
+you can output to a temporary @code{ostringstream} like this:
- cout << "cl_N n = ";
- e.print(print_csrc_cl_N(cout));
- cout << ";\n";
+@example
+ // ...
+ ostringstream s;
+ s << latex << e; // format of cout remains unchanged
+ cout << s.str() << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
// ...
@end example
-The three possible types mostly affect the way in which floating point
-numbers are written.
+@cindex @code{csrc}
+@cindex @code{csrc_float}
+@cindex @code{csrc_double}
+@cindex @code{csrc_cl_N}
+The @code{csrc} (an alias for @code{csrc_double}), @code{csrc_float},
+@code{csrc_double} and @code{csrc_cl_N} manipulators set the output to a
+format that can be directly used in a C or C++ program. The three possible
+formats select the data types used for numbers (@code{csrc_cl_N} uses the
+classes provided by the CLN library):
+
+@example
+ // ...
+ cout << "f = " << csrc_float << e << ";\n";
+ cout << "d = " << csrc_double << e << ";\n";
+ cout << "n = " << csrc_cl_N << e << ";\n";
+ // ...
+@end example
-The above example will produce (note the @code{x^2} being converted to @code{x*x}):
+The above example will produce (note the @code{x^2} being converted to
+@code{x*x}):
@example
-float f = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
-double d = (3.000000e+00/2.000000e+00)*(x*x)+4.500000e+00;
-cl_N n = (cln::cl_F("3.0")/cln::cl_F("2.0"))*(x*x)+cln::cl_F("4.5");
+f = (3.0/2.0)*(x*x)+std::complex<float>(0.0,4.5000000e+00);
+d = (3.0/2.0)*(x*x)+std::complex<double>(0.0,4.5000000000000000e+00);
+n = cln::cl_RA("3/2")*(x*x)+cln::complex(cln::cl_I("0"),cln::cl_F("4.5_17"));
@end example
-The @code{print_context} type @code{print_tree} provides a dump of the
-internal structure of an expression for debugging purposes:
+@cindex @code{tree}
+The @code{tree} manipulator allows dumping the internal structure of an
+expression for debugging purposes:
@example
// ...
- e.print(print_tree(cout));
+ cout << tree << e;
@}
@end example
@example
add, hash=0x0, flags=0x3, nops=2
- power, hash=0x9, flags=0x3, nops=2
- x (symbol), serial=3, hash=0x44a113a6, flags=0xf
- 2 (numeric), hash=0x80000042, flags=0xf
- 3/2 (numeric), hash=0x80000061, flags=0xf
+ power, hash=0x0, flags=0x3, nops=2
+ x (symbol), serial=0, hash=0xc8d5bcdd, flags=0xf
+ 2 (numeric), hash=0x6526b0fa, flags=0xf
+ 3/2 (numeric), hash=0xf9828fbd, flags=0xf
-----
overall_coeff
- 4.5L0 (numeric), hash=0x8000004b, flags=0xf
+ 4.5L0i (numeric), hash=0xa40a97e0, flags=0xf
=====
@end example
-This kind of output is also available in @command{ginsh} as the @code{print()}
-function.
-
-Another useful output format is for LaTeX parsing in mathematical mode.
-It is rather similar to the default @code{print_context} but provides
-some braces needed by LaTeX for delimiting boxes and also converts some
-common objects to conventional LaTeX names. It is possible to give symbols
-a special name for LaTeX output by supplying it as a second argument to
-the @code{symbol} constructor.
+@cindex @code{latex}
+The @code{latex} output format is for LaTeX parsing in mathematical mode.
+It is rather similar to the default format but provides some braces needed
+by LaTeX for delimiting boxes and also converts some common objects to
+conventional LaTeX names. It is possible to give symbols a special name for
+LaTeX output by supplying it as a second argument to the @code{symbol}
+constructor.
For example, the code snippet
@example
- // ...
- symbol x("x");
- ex foo = lgamma(x).series(x==0,3);
- foo.print(print_latex(std::cout));
+@{
+ symbol x("x", "\\circ");
+ ex e = lgamma(x).series(x==0,3);
+ cout << latex << e << endl;
+@}
+@end example
+
+will print
+
+@example
+ @{(-\ln(\circ))@}+@{(-\gamma_E)@} \circ+@{(\frac@{1@}@{12@} \pi^@{2@})@} \circ^@{2@}+\mathcal@{O@}(\circ^@{3@})
@end example
-will print out:
+@cindex @code{index_dimensions}
+@cindex @code{no_index_dimensions}
+Index dimensions are normally hidden in the output. To make them visible, use
+the @code{index_dimensions} manipulator. The dimensions will be written in
+square brackets behind each index value in the default and LaTeX output
+formats:
@example
- @{(-\ln(x))@}+@{(-\gamma_E)@} x+@{(1/12 \pi^2)@} x^@{2@}+\mathcal@{O@}(x^3)
+@{
+ symbol x("x"), y("y");
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
+ ex e = indexed(x, mu) * indexed(y, nu);
+
+ cout << e << endl;
+ // prints 'x~mu*y~nu'
+ cout << index_dimensions << e << endl;
+ // prints 'x~mu[4]*y~nu[4]'
+ cout << no_index_dimensions << e << endl;
+ // prints 'x~mu*y~nu'
+@}
@end example
+
@cindex Tree traversal
If you need any fancy special output format, e.g. for interfacing GiNaC
with other algebra systems or for producing code for different
else
cout << e.bp->class_name();
cout << "(";
- unsigned n = e.nops();
+ size_t n = e.nops();
if (n)
- for (unsigned i=0; i<n; i++) @{
+ for (size_t i=0; i<n; i++) @{
my_print(e.op(i));
if (i != n-1)
cout << ",";
cout << ")";
@}
-int main(void)
+int main()
@{
my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
return 0;
int main()
@{
- symbol x("x");
- string s;
-
- cout << "Enter an expression containing 'x': ";
- getline(cin, s);
-
- try @{
- ex e(s, lst(x));
- cout << "The derivative of " << e << " with respect to x is ";
- cout << e.diff(x) << ".\n";
- @} catch (exception &p) @{
- cerr << p.what() << endl;
- @}
+ symbol x("x");
+ string s;
+
+ cout << "Enter an expression containing 'x': ";
+ getline(cin, s);
+
+ try @{
+ ex e(s, lst(x));
+ cout << "The derivative of " << e << " with respect to x is ";
+ cout << e.diff(x) << ".\n";
+ @} catch (exception &p) @{
+ cerr << p.what() << endl;
+ @}
@}
@end example
the @code{syms} list above, the @code{ex1.subs(x == 2)} statement would
have had no effect because the @code{x} in @code{ex1} would have been a
different symbol than the @code{x} which was defined at the beginning of
-the program, altough both would appear as @samp{x} when printed.
+the program, although both would appear as @samp{x} when printed.
You can also use the information stored in an @code{archive} object to
output expressions in a format suitable for exact reconstruction. The
archive_node::propinfovector p;
n.get_properties(p);
- unsigned num = p.size();
- for (unsigned i=0; i<num; i++) @{
+ size_t num = p.size();
+ for (size_t i=0; i<num; i++) @{
const string &name = p[i].name;
if (name == "class")
continue;
switch (p[i].type) @{
case archive_node::PTYPE_BOOL: @{
bool x;
- n.find_bool(name, x);
+ n.find_bool(name, x, j);
cout << (x ? "true" : "false");
break;
@}
case archive_node::PTYPE_UNSIGNED: @{
unsigned x;
- n.find_unsigned(name, x);
+ n.find_unsigned(name, x, j);
cout << x;
break;
@}
case archive_node::PTYPE_STRING: @{
string x;
- n.find_string(name, x);
+ n.find_string(name, x, j);
cout << '\"' << x << '\"';
break;
@}
cout << ")";
@}
-int main(void)
+int main()
@{
ex e = pow(2, x) - y;
archive ar(e, "e");
@menu
* What does not belong into GiNaC:: What to avoid.
* Symbolic functions:: Implementing symbolic functions.
-* Adding classes:: Defining new algebraic classes.
+* Structures:: Defining new algebraic classes (the easy way).
+* Adding classes:: Defining new algebraic classes (the hard way).
@end menu
language. There are no loops or conditional expressions in
@command{ginsh}, it is merely a window into the library for the
programmer to test stuff (or to show off). Still, the design of a
-complete CAS with a language of its own, graphical capabilites and all
+complete CAS with a language of its own, graphical capabilities and all
this on top of GiNaC is possible and is without doubt a nice project for
the future.
generally. This ought to be fixed. However, doing numerical
computations with GiNaC's quite abstract classes is doomed to be
inefficient. For this purpose, the underlying foundation classes
-provided by @acronym{CLN} are much better suited.
+provided by CLN are much better suited.
-@node Symbolic functions, Adding classes, What does not belong into GiNaC, Extending GiNaC
+@node Symbolic functions, Structures, What does not belong into GiNaC, Extending GiNaC
@c node-name, next, previous, up
@section Symbolic functions
-The easiest and most instructive way to start with is probably to
-implement your own function. GiNaC's functions are objects of class
-@code{function}. The preprocessor is then used to convert the function
-names to objects with a corresponding serial number that is used
-internally to identify them. You usually need not worry about this
-number. New functions may be inserted into the system via a kind of
-`registry'. It is your responsibility to care for some functions that
-are called when the user invokes certain methods. These are usual
-C++-functions accepting a number of @code{ex} as arguments and returning
-one @code{ex}. As an example, if we have a look at a simplified
-implementation of the cosine trigonometric function, we first need a
-function that is called when one wishes to @code{eval} it. It could
-look something like this:
-
-@example
-static ex cos_eval_method(const ex & x)
+The easiest and most instructive way to start extending GiNaC is probably to
+create your own symbolic functions. These are implemented with the help of
+two preprocessor macros:
+
+@cindex @code{DECLARE_FUNCTION}
+@cindex @code{REGISTER_FUNCTION}
+@example
+DECLARE_FUNCTION_<n>P(<name>)
+REGISTER_FUNCTION(<name>, <options>)
+@end example
+
+The @code{DECLARE_FUNCTION} macro will usually appear in a header file. It
+declares a C++ function with the given @samp{name} that takes exactly @samp{n}
+parameters of type @code{ex} and returns a newly constructed GiNaC
+@code{function} object that represents your function.
+
+The @code{REGISTER_FUNCTION} macro implements the function. It must be passed
+the same @samp{name} as the respective @code{DECLARE_FUNCTION} macro, and a
+set of options that associate the symbolic function with C++ functions you
+provide to implement the various methods such as evaluation, derivative,
+series expansion etc. They also describe additional attributes the function
+might have, such as symmetry and commutation properties, and a name for
+LaTeX output. Multiple options are separated by the member access operator
+@samp{.} and can be given in an arbitrary order.
+
+(By the way: in case you are worrying about all the macros above we can
+assure you that functions are GiNaC's most macro-intense classes. We have
+done our best to avoid macros where we can.)
+
+@subsection A minimal example
+
+Here is an example for the implementation of a function with two arguments
+that is not further evaluated:
+
+@example
+DECLARE_FUNCTION_2P(myfcn)
+
+static ex myfcn_eval(const ex & x, const ex & y)
@{
- // if (!x%(2*Pi)) return 1
- // if (!x%Pi) return -1
- // if (!x%Pi/2) return 0
- // care for other cases...
- return cos(x).hold();
+ return myfcn(x, y).hold();
+@}
+
+REGISTER_FUNCTION(myfcn, eval_func(myfcn_eval))
+@end example
+
+Any code that has seen the @code{DECLARE_FUNCTION} line can use @code{myfcn()}
+in algebraic expressions:
+
+@example
+@{
+ ...
+ symbol x("x");
+ ex e = 2*myfcn(42, 3*x+1) - x;
+ // this calls myfcn_eval(42, 3*x+1), and inserts its return value into
+ // the actual expression
+ cout << e << endl;
+ // prints '2*myfcn(42,1+3*x)-x'
+ ...
@}
@end example
@cindex @code{hold()}
@cindex evaluation
-The last line returns @code{cos(x)} if we don't know what else to do and
-stops a potential recursive evaluation by saying @code{.hold()}, which
-sets a flag to the expression signaling that it has been evaluated. We
-should also implement a method for numerical evaluation and since we are
-lazy we sweep the problem under the rug by calling someone else's
-function that does so, in this case the one in class @code{numeric}:
+The @code{eval_func()} option specifies the C++ function that implements
+the @code{eval()} method, GiNaC's anonymous evaluator. This function takes
+the same number of arguments as the associated symbolic function (two in this
+case) and returns the (possibly transformed or in some way simplified)
+symbolically evaluated function (@xref{Automatic evaluation}, for a description
+of the automatic evaluation process). If no (further) evaluation is to take
+place, the @code{eval_func()} function must return the original function
+with @code{.hold()}, to avoid a potential infinite recursion. If your
+symbolic functions produce a segmentation fault or stack overflow when
+using them in expressions, you are probably missing a @code{.hold()}
+somewhere.
+
+There is not much you can do with the @code{myfcn} function. It merely acts
+as a kind of container for its arguments (which is, however, sometimes
+perfectly sufficient). Let's have a look at the implementation of GiNaC's
+cosine function.
+
+@subsection The cosine function
+
+The GiNaC header file @file{inifcns.h} contains the line
+
+@example
+DECLARE_FUNCTION_1P(cos)
+@end example
+
+which declares to all programs using GiNaC that there is a function @samp{cos}
+that takes one @code{ex} as an argument. This is all they need to know to use
+this function in expressions.
+
+The implementation of the cosine function is in @file{inifcns_trans.cpp}. The
+@code{eval_func()} function looks something like this (actually, it doesn't
+look like this at all, but it should give you an idea what is going on):
+
+@example
+static ex cos_eval(const ex & x)
+@{
+ if (<x is a multiple of 2*Pi>)
+ return 1;
+ else if (<x is a multiple of Pi>)
+ return -1;
+ else if (<x is a multiple of Pi/2>)
+ return 0;
+ // more rules...
+
+ else if (<x has the form 'acos(y)'>)
+ return y;
+ else if (<x has the form 'asin(y)'>)
+ return sqrt(1-y^2);
+ // more rules...
+
+ else
+ return cos(x).hold();
+@}
+@end example
+
+In this way, @code{cos(4*Pi)} automatically becomes @math{1},
+@code{cos(asin(a+b))} becomes @code{sqrt(1-(a+b)^2)}, etc. If no reasonable
+symbolic transformation can be done, the unmodified function is returned
+with @code{.hold()}.
+
+GiNaC doesn't automatically transform @code{cos(2)} to @samp{-0.416146...}.
+The user has to call @code{evalf()} for that. This is implemented in a
+different function:
@example
static ex cos_evalf(const ex & x)
@{
- return cos(ex_to<numeric>(x));
+ if (is_a<numeric>(x))
+ return cos(ex_to<numeric>(x));
+ else
+ return cos(x).hold();
@}
@end example
+Since we are lazy we defer the problem of numeric evaluation to somebody else,
+in this case the @code{cos()} function for @code{numeric} objects, which in
+turn hands it over to the @code{cos()} function in CLN. The @code{.hold()}
+isn't really needed here, but reminds us that the corresponding @code{eval()}
+function would require it in this place.
+
Differentiation will surely turn up and so we need to tell @code{cos}
-what the first derivative is (higher derivatives (@code{.diff(x,3)} for
-instance are then handled automatically by @code{basic::diff} and
+what its first derivative is (higher derivatives, @code{.diff(x,3)} for
+instance, are then handled automatically by @code{basic::diff} and
@code{ex::diff}):
@example
@cindex product rule
The second parameter is obligatory but uninteresting at this point. It
specifies which parameter to differentiate in a partial derivative in
-case the function has more than one parameter and its main application
-is for correct handling of the chain rule. For Taylor expansion, it is
-enough to know how to differentiate. But if the function you want to
-implement does have a pole somewhere in the complex plane, you need to
-write another method for Laurent expansion around that point.
+case the function has more than one parameter, and its main application
+is for correct handling of the chain rule.
+
+An implementation of the series expansion is not needed for @code{cos()} as
+it doesn't have any poles and GiNaC can do Taylor expansion by itself (as
+long as it knows what the derivative of @code{cos()} is). @code{tan()}, on
+the other hand, does have poles and may need to do Laurent expansion:
-Now that all the ingredients for @code{cos} have been set up, we need
-to tell the system about it. This is done by a macro and we are not
-going to descibe how it expands, please consult your preprocessor if you
-are curious:
+@example
+static ex tan_series(const ex & x, const relational & rel,
+ int order, unsigned options)
+@{
+ // Find the actual expansion point
+ const ex x_pt = x.subs(rel);
+
+ if (<x_pt is not an odd multiple of Pi/2>)
+ throw do_taylor(); // tell function::series() to do Taylor expansion
+
+ // On a pole, expand sin()/cos()
+ return (sin(x)/cos(x)).series(rel, order+2, options);
+@}
+@end example
+
+The @code{series()} implementation of a function @emph{must} return a
+@code{pseries} object, otherwise your code will crash.
+
+Now that all the ingredients have been set up, the @code{REGISTER_FUNCTION}
+macro is used to tell the system how the @code{cos()} function behaves:
@example
REGISTER_FUNCTION(cos, eval_func(cos_eval).
evalf_func(cos_evalf).
- derivative_func(cos_deriv));
-@end example
-
-The first argument is the function's name used for calling it and for
-output. The second binds the corresponding methods as options to this
-object. Options are separated by a dot and can be given in an arbitrary
-order. GiNaC functions understand several more options which are always
-specified as @code{.option(params)}, for example a method for series
-expansion @code{.series_func(cos_series)}. Again, if no series
-expansion method is given, GiNaC defaults to simple Taylor expansion,
-which is correct if there are no poles involved as is the case for the
-@code{cos} function. The way GiNaC handles poles in case there are any
-is best understood by studying one of the examples, like the Gamma
-(@code{tgamma}) function for instance. (In essence the function first
-checks if there is a pole at the evaluation point and falls back to
-Taylor expansion if there isn't. Then, the pole is regularized by some
-suitable transformation.) Also, the new function needs to be declared
-somewhere. This may also be done by a convenient preprocessor macro:
+ derivative_func(cos_deriv).
+ latex_name("\\cos"));
+@end example
+
+This registers the @code{cos_eval()}, @code{cos_evalf()} and
+@code{cos_deriv()} C++ functions with the @code{cos()} function, and also
+gives it a proper LaTeX name.
+
+@subsection Function options
+
+GiNaC functions understand several more options which are always
+specified as @code{.option(params)}. None of them are required, but you
+need to specify at least one option to @code{REGISTER_FUNCTION()} (usually
+the @code{eval()} method).
@example
-DECLARE_FUNCTION_1P(cos)
+eval_func(<C++ function>)
+evalf_func(<C++ function>)
+derivative_func(<C++ function>)
+series_func(<C++ function>)
@end example
-The suffix @code{_1P} stands for @emph{one parameter}. Of course, this
-implementation of @code{cos} is very incomplete and lacks several safety
-mechanisms. Please, have a look at the real implementation in GiNaC.
-(By the way: in case you are worrying about all the macros above we can
-assure you that functions are GiNaC's most macro-intense classes. We
-have done our best to avoid macros where we can.)
+These specify the C++ functions that implement symbolic evaluation,
+numeric evaluation, partial derivatives, and series expansion, respectively.
+They correspond to the GiNaC methods @code{eval()}, @code{evalf()},
+@code{diff()} and @code{series()}.
+The @code{eval_func()} function needs to use @code{.hold()} if no further
+automatic evaluation is desired or possible.
-@node Adding classes, A Comparison With Other CAS, Symbolic functions, Extending GiNaC
+If no @code{series_func()} is given, GiNaC defaults to simple Taylor
+expansion, which is correct if there are no poles involved. If the function
+has poles in the complex plane, the @code{series_func()} needs to check
+whether the expansion point is on a pole and fall back to Taylor expansion
+if it isn't. Otherwise, the pole usually needs to be regularized by some
+suitable transformation.
+
+@example
+latex_name(const string & n)
+@end example
+
+specifies the LaTeX code that represents the name of the function in LaTeX
+output. The default is to put the function name in an @code{\mbox@{@}}.
+
+@example
+do_not_evalf_params()
+@end example
+
+This tells @code{evalf()} to not recursively evaluate the parameters of the
+function before calling the @code{evalf_func()}.
+
+@example
+set_return_type(unsigned return_type, unsigned return_type_tinfo)
+@end example
+
+This allows you to explicitly specify the commutation properties of the
+function (@xref{Non-commutative objects}, for an explanation of
+(non)commutativity in GiNaC). For example, you can use
+@code{set_return_type(return_types::noncommutative, TINFO_matrix)} to make
+GiNaC treat your function like a matrix. By default, functions inherit the
+commutation properties of their first argument.
+
+@example
+set_symmetry(const symmetry & s)
+@end example
+
+specifies the symmetry properties of the function with respect to its
+arguments. @xref{Indexed objects}, for an explanation of symmetry
+specifications. GiNaC will automatically rearrange the arguments of
+symmetric functions into a canonical order.
+
+
+@node Structures, Adding classes, Symbolic functions, Extending GiNaC
+@c node-name, next, previous, up
+@section Structures
+
+If you are doing some very specialized things with GiNaC, or if you just
+need some more organized way to store data in your expressions instead of
+anonymous lists, you may want to implement your own algebraic classes.
+('algebraic class' means any class directly or indirectly derived from
+@code{basic} that can be used in GiNaC expressions).
+
+GiNaC offers two ways of accomplishing this: either by using the
+@code{structure<T>} template class, or by rolling your own class from
+scratch. This section will discuss the @code{structure<T>} template which
+is easier to use but more limited, while the implementation of custom
+GiNaC classes is the topic of the next section. However, you may want to
+read both sections because many common concepts and member functions are
+shared by both concepts, and it will also allow you to decide which approach
+is most suited to your needs.
+
+The @code{structure<T>} template, defined in the GiNaC header file
+@file{structure.h}, wraps a type that you supply (usually a C++ @code{struct}
+or @code{class}) into a GiNaC object that can be used in expressions.
+
+@subsection Example: scalar products
+
+Let's suppose that we need a way to handle some kind of abstract scalar
+product of the form @samp{<x|y>} in expressions. Objects of the scalar
+product class have to store their left and right operands, which can in turn
+be arbitrary expressions. Here is a possible way to represent such a
+product in a C++ @code{struct}:
+
+@example
+#include <iostream>
+using namespace std;
+
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+struct sprod_s @{
+ ex left, right;
+
+ sprod_s() @{@}
+ sprod_s(ex l, ex r) : left(l), right(r) @{@}
+@};
+@end example
+
+The default constructor is required. Now, to make a GiNaC class out of this
+data structure, we need only one line:
+
+@example
+typedef structure<sprod_s> sprod;
+@end example
+
+That's it. This line constructs an algebraic class @code{sprod} which
+contains objects of type @code{sprod_s}. We can now use @code{sprod} in
+expressions like any other GiNaC class:
+
+@example
+...
+ symbol a("a"), b("b");
+ ex e = sprod(sprod_s(a, b));
+...
+@end example
+
+Note the difference between @code{sprod} which is the algebraic class, and
+@code{sprod_s} which is the unadorned C++ structure containing the @code{left}
+and @code{right} data members. As shown above, an @code{sprod} can be
+constructed from an @code{sprod_s} object.
+
+If you find the nested @code{sprod(sprod_s())} constructor too unwieldy,
+you could define a little wrapper function like this:
+
+@example
+inline ex make_sprod(ex left, ex right)
+@{
+ return sprod(sprod_s(left, right));
+@}
+@end example
+
+The @code{sprod_s} object contained in @code{sprod} can be accessed with
+the GiNaC @code{ex_to<>()} function followed by the @code{->} operator or
+@code{get_struct()}:
+
+@example
+...
+ cout << ex_to<sprod>(e)->left << endl;
+ // -> a
+ cout << ex_to<sprod>(e).get_struct().right << endl;
+ // -> b
+...
+@end example
+
+You only have read access to the members of @code{sprod_s}.
+
+The type definition of @code{sprod} is enough to write your own algorithms
+that deal with scalar products, for example:
+
+@example
+ex swap_sprod(ex p)
+@{
+ if (is_a<sprod>(p)) @{
+ const sprod_s & sp = ex_to<sprod>(p).get_struct();
+ return make_sprod(sp.right, sp.left);
+ @} else
+ return p;
+@}
+
+...
+ f = swap_sprod(e);
+ // f is now <b|a>
+...
+@end example
+
+@subsection Structure output
+
+While the @code{sprod} type is useable it still leaves something to be
+desired, most notably proper output:
+
+@example
+...
+ cout << e << endl;
+ // -> [structure object]
+...
+@end example
+
+By default, any structure types you define will be printed as
+@samp{[structure object]}. To override this, you can specialize the
+template's @code{print()} member function. The member functions of
+GiNaC classes are described in more detail in the next section, but
+it shouldn't be hard to figure out what's going on here:
+
+@example
+void sprod::print(const print_context & c, unsigned level) const
+@{
+ // tree debug output handled by superclass
+ if (is_a<print_tree>(c))
+ inherited::print(c, level);
+
+ // get the contained sprod_s object
+ const sprod_s & sp = get_struct();
+
+ // print_context::s is a reference to an ostream
+ c.s << "<" << sp.left << "|" << sp.right << ">";
+@}
+@end example
+
+Now we can print expressions containing scalar products:
+
+@example
+...
+ cout << e << endl;
+ // -> <a|b>
+ cout << swap_sprod(e) << endl;
+ // -> <b|a>
+...
+@end example
+
+@subsection Comparing structures
+
+The @code{sprod} class defined so far still has one important drawback: all
+scalar products are treated as being equal because GiNaC doesn't know how to
+compare objects of type @code{sprod_s}. This can lead to some confusing
+and undesired behavior:
+
+@example
+...
+ cout << make_sprod(a, b) - make_sprod(a*a, b*b) << endl;
+ // -> 0
+ cout << make_sprod(a, b) + make_sprod(a*a, b*b) << endl;
+ // -> 2*<a|b> or 2*<a^2|b^2> (which one is undefined)
+...
+@end example
+
+To remedy this, we first need to define the operators @code{==} and @code{<}
+for objects of type @code{sprod_s}:
+
+@example
+inline bool operator==(const sprod_s & lhs, const sprod_s & rhs)
+@{
+ return lhs.left.is_equal(rhs.left) && lhs.right.is_equal(rhs.right);
+@}
+
+inline bool operator<(const sprod_s & lhs, const sprod_s & rhs)
+@{
+ return lhs.left.compare(rhs.left) < 0 ? true : lhs.right.compare(rhs.right) < 0;
+@}
+@end example
+
+The ordering established by the @code{<} operator doesn't have to make any
+algebraic sense, but it needs to be well defined. Note that we can't use
+expressions like @code{lhs.left == rhs.left} or @code{lhs.left < rhs.left}
+in the implementation of these operators because they would construct
+GiNaC @code{relational} objects which in the case of @code{<} do not
+establish a well defined ordering (for arbitrary expressions, GiNaC can't
+decide which one is algebraically 'less').
+
+Next, we need to change our definition of the @code{sprod} type to let
+GiNaC know that an ordering relation exists for the embedded objects:
+
+@example
+typedef structure<sprod_s, compare_std_less> sprod;
+@end example
+
+@code{sprod} objects then behave as expected:
+
+@example
+...
+ cout << make_sprod(a, b) - make_sprod(a*a, b*b) << endl;
+ // -> <a|b>-<a^2|b^2>
+ cout << make_sprod(a, b) + make_sprod(a*a, b*b) << endl;
+ // -> <a|b>+<a^2|b^2>
+ cout << make_sprod(a, b) - make_sprod(a, b) << endl;
+ // -> 0
+ cout << make_sprod(a, b) + make_sprod(a, b) << endl;
+ // -> 2*<a|b>
+...
+@end example
+
+The @code{compare_std_less} policy parameter tells GiNaC to use the
+@code{std::less} and @code{std::equal_to} functors to compare objects of
+type @code{sprod_s}. By default, these functors forward their work to the
+standard @code{<} and @code{==} operators, which we have overloaded.
+Alternatively, we could have specialized @code{std::less} and
+@code{std::equal_to} for class @code{sprod_s}.
+
+GiNaC provides two other comparison policies for @code{structure<T>}
+objects: the default @code{compare_all_equal}, and @code{compare_bitwise}
+which does a bit-wise comparison of the contained @code{T} objects.
+This should be used with extreme care because it only works reliably with
+built-in integral types, and it also compares any padding (filler bytes of
+undefined value) that the @code{T} class might have.
+
+@subsection Subexpressions
+
+Our scalar product class has two subexpressions: the left and right
+operands. It might be a good idea to make them accessible via the standard
+@code{nops()} and @code{op()} methods:
+
+@example
+size_t sprod::nops() const
+@{
+ return 2;
+@}
+
+ex sprod::op(size_t i) const
+@{
+ switch (i) @{
+ case 0:
+ return get_struct().left;
+ case 1:
+ return get_struct().right;
+ default:
+ throw std::range_error("sprod::op(): no such operand");
+ @}
+@}
+@end example
+
+Implementing @code{nops()} and @code{op()} for container types such as
+@code{sprod} has two other nice side effects:
+
+@itemize @bullet
+@item
+@code{has()} works as expected
+@item
+GiNaC generates better hash keys for the objects (the default implementation
+of @code{calchash()} takes subexpressions into account)
+@end itemize
+
+@cindex @code{let_op()}
+There is a non-const variant of @code{op()} called @code{let_op()} that
+allows replacing subexpressions:
+
+@example
+ex & sprod::let_op(size_t i)
+@{
+ // every non-const member function must call this
+ ensure_if_modifiable();
+
+ switch (i) @{
+ case 0:
+ return get_struct().left;
+ case 1:
+ return get_struct().right;
+ default:
+ throw std::range_error("sprod::let_op(): no such operand");
+ @}
+@}
+@end example
+
+Once we have provided @code{let_op()} we also get @code{subs()} and
+@code{map()} for free. In fact, every container class that returns a non-null
+@code{nops()} value must either implement @code{let_op()} or provide custom
+implementations of @code{subs()} and @code{map()}.
+
+In turn, the availability of @code{map()} enables the recursive behavior of a
+couple of other default method implementations, in particular @code{evalf()},
+@code{evalm()}, @code{normal()}, @code{diff()} and @code{expand()}. Although
+we probably want to provide our own version of @code{expand()} for scalar
+products that turns expressions like @samp{<a+b|c>} into @samp{<a|c>+<b|c>}.
+This is left as an exercise for the reader.
+
+The @code{structure<T>} template defines many more member functions that
+you can override by specialization to customize the behavior of your
+structures. You are referred to the next section for a description of
+some of these (especially @code{eval()}). There is, however, one topic
+that shall be addressed here, as it demonstrates one peculiarity of the
+@code{structure<T>} template: archiving.
+
+@subsection Archiving structures
+
+If you don't know how the archiving of GiNaC objects is implemented, you
+should first read the next section and then come back here. You're back?
+Good.
+
+To implement archiving for structures it is not enough to provide
+specializations for the @code{archive()} member function and the
+unarchiving constructor (the @code{unarchive()} function has a default
+implementation). You also need to provide a unique name (as a string literal)
+for each structure type you define. This is because in GiNaC archives,
+the class of an object is stored as a string, the class name.
+
+By default, this class name (as returned by the @code{class_name()} member
+function) is @samp{structure} for all structure classes. This works as long
+as you have only defined one structure type, but if you use two or more you
+need to provide a different name for each by specializing the
+@code{get_class_name()} member function. Here is a sample implementation
+for enabling archiving of the scalar product type defined above:
+
+@example
+const char *sprod::get_class_name() @{ return "sprod"; @}
+
+void sprod::archive(archive_node & n) const
+@{
+ inherited::archive(n);
+ n.add_ex("left", get_struct().left);
+ n.add_ex("right", get_struct().right);
+@}
+
+sprod::structure(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst)
+@{
+ n.find_ex("left", get_struct().left, sym_lst);
+ n.find_ex("right", get_struct().right, sym_lst);
+@}
+@end example
+
+Note that the unarchiving constructor is @code{sprod::structure} and not
+@code{sprod::sprod}, and that we don't need to supply an
+@code{sprod::unarchive()} function.
+
+
+@node Adding classes, A Comparison With Other CAS, Structures, Extending GiNaC
@c node-name, next, previous, up
@section Adding classes
-If you are doing some very specialized things with GiNaC you may find that
-you have to implement your own algebraic classes to fit your needs. This
-section will explain how to do this by giving the example of a simple
-'string' class. After reading this section you will know how to properly
-declare a GiNaC class and what the minimum required member functions are
-that you have to implement. We only cover the implementation of a 'leaf'
-class here (i.e. one that doesn't contain subexpressions). Creating a
-container class like, for example, a class representing tensor products is
-more involved but this section should give you enough information so you can
-consult the source to GiNaC's predefined classes if you want to implement
-something more complicated.
+The @code{structure<T>} template provides an way to extend GiNaC with custom
+algebraic classes that is easy to use but has its limitations, the most
+severe of which being that you can't add any new member functions to
+structures. To be able to do this, you need to write a new class definition
+from scratch.
+
+This section will explain how to implement new algebraic classes in GiNaC by
+giving the example of a simple 'string' class. After reading this section
+you will know how to properly declare a GiNaC class and what the minimum
+required member functions are that you have to implement. We only cover the
+implementation of a 'leaf' class here (i.e. one that doesn't contain
+subexpressions). Creating a container class like, for example, a class
+representing tensor products is more involved but this section should give
+you enough information so you can consult the source to GiNaC's predefined
+classes if you want to implement something more complicated.
@subsection GiNaC's run-time type information system
@end example
The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
-macros are defined in @file{registrar.h}. They take the name of the class
+macros are defined in @file{registrar.h}. They take the name of the class
and its direct superclass as arguments and insert all required declarations
for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
the first line after the opening brace of the class definition. The
source (at global scope, of course, not inside a function).
@code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
-declarations of the default and copy constructor, the destructor, the
-assignment operator and a couple of other functions that are required. It
-also defines a type @code{inherited} which refers to the superclass so you
-don't have to modify your code every time you shuffle around the class
-hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} implements the copy
-constructor, the destructor and the assignment operator.
-
-Now there are nine member functions we have to implement to get a working
+declarations of the default constructor and a couple of other functions that
+are required. It also defines a type @code{inherited} which refers to the
+superclass so you don't have to modify your code every time you shuffle around
+the class hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} registers the
+class with the GiNaC RTTI.
+
+Now there are seven member functions we have to implement to get a working
class:
@itemize
@item
@code{mystring()}, the default constructor.
-@item
-@code{void destroy(bool call_parent)}, which is used in the destructor and the
-assignment operator to free dynamically allocated members. The @code{call_parent}
-specifies whether the @code{destroy()} function of the superclass is to be
-called also.
-
-@item
-@code{void copy(const mystring &other)}, which is used in the copy constructor
-and assignment operator to copy the member variables over from another
-object of the same class.
-
@item
@code{void archive(archive_node &n)}, the archiving function. This stores all
information needed to reconstruct an object of this class inside an
@code{archive_node}.
@item
-@code{mystring(const archive_node &n, const lst &sym_lst)}, the unarchiving
+@code{mystring(const archive_node &n, lst &sym_lst)}, the unarchiving
constructor. This constructs an instance of the class from the information
found in an @code{archive_node}.
@item
-@code{ex unarchive(const archive_node &n, const lst &sym_lst)}, the static
+@code{ex unarchive(const archive_node &n, lst &sym_lst)}, the static
unarchiving function. It constructs a new instance by calling the unarchiving
constructor.
@item
+@cindex @code{compare_same_type()}
@code{int compare_same_type(const basic &other)}, which is used internally
by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
-1, depending on the relative order of this object and the @code{other}
Let's proceed step-by-step. The default constructor looks like this:
@example
-mystring::mystring() : inherited(TINFO_mystring)
-@{
- // dynamically allocate resources here if required
-@}
+mystring::mystring() : inherited(TINFO_mystring) @{@}
@end example
The golden rule is that in all constructors you have to set the
it will be set by the constructor of the superclass and all hell will break
loose in the RTTI. For your convenience, the @code{basic} class provides
a constructor that takes a @code{tinfo_key} value, which we are using here
-(remember that in our case @code{inherited = basic}). If the superclass
+(remember that in our case @code{inherited == basic}). If the superclass
didn't have such a constructor, we would have to set the @code{tinfo_key}
to the right value manually.
In the default constructor you should set all other member variables to
reasonable default values (we don't need that here since our @code{str}
-member gets set to an empty string automatically). The constructor(s) are of
-course also the right place to allocate any dynamic resources you require.
-
-Next, the @code{destroy()} function:
-
-@example
-void mystring::destroy(bool call_parent)
-@{
- // free dynamically allocated resources here if required
- if (call_parent)
- inherited::destroy(call_parent);
-@}
-@end example
-
-This function is where we free all dynamically allocated resources. We don't
-have any so we're not doing anything here, but if we had, for example, used
-a C-style @code{char *} to store our string, this would be the place to
-@code{delete[]} the string storage. If @code{call_parent} is true, we have
-to call the @code{destroy()} function of the superclass after we're done
-(to mimic C++'s automatic invocation of superclass destructors where
-@code{destroy()} is called from outside a destructor).
-
-The @code{copy()} function just copies over the member variables from
-another object:
-
-@example
-void mystring::copy(const mystring &other)
-@{
- inherited::copy(other);
- str = other.str;
-@}
-@end example
-
-We can simply overwrite the member variables here. There's no need to worry
-about dynamically allocated storage. The assignment operator (which is
-automatically defined by @code{GINAC_IMPLEMENT_REGISTERED_CLASS}, as you
-recall) calls @code{destroy()} before it calls @code{copy()}. You have to
-explicitly call the @code{copy()} function of the superclass here so
-all the member variables will get copied.
+member gets set to an empty string automatically).
Next are the three functions for archiving. You have to implement them even
if you don't plan to use archives, but the minimum required implementation
-is really simple. First, the archiving function:
+is really simple. First, the archiving function:
@example
void mystring::archive(archive_node &n) const
The only thing that is really required is calling the @code{archive()}
function of the superclass. Optionally, you can store all information you
deem necessary for representing the object into the passed
-@code{archive_node}. We are just storing our string here. For more
+@code{archive_node}. We are just storing our string here. For more
information on how the archiving works, consult the @file{archive.h} header
file.
function:
@example
-mystring::mystring(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
+mystring::mystring(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
@{
n.find_string("string", str);
@}
Finally, the unarchiving function:
@example
-ex mystring::unarchive(const archive_node &n, const lst &sym_lst)
+ex mystring::unarchive(const archive_node &n, lst &sym_lst)
@{
return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
@}
@end example
-You don't have to understand how exactly this works. Just copy these four
-lines into your code literally (replacing the class name, of course). It
-calls the unarchiving constructor of the class and unless you are doing
-something very special (like matching @code{archive_node}s to global
-objects) you don't need a different implementation. For those who are
-interested: setting the @code{dynallocated} flag puts the object under
-the control of GiNaC's garbage collection. It will get deleted automatically
-once it is no longer referenced.
+You don't have to understand how exactly this works. Just copy these
+four lines into your code literally (replacing the class name, of
+course). It calls the unarchiving constructor of the class and unless
+you are doing something very special (like matching @code{archive_node}s
+to global objects) you don't need a different implementation. For those
+who are interested: setting the @code{dynallocated} flag puts the object
+under the control of GiNaC's garbage collection. It will get deleted
+automatically once it is no longer referenced.
Our @code{compare_same_type()} function uses a provided function to compare
the string members:
Now the only thing missing is our two new constructors:
@example
-mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s)
-@{
- // dynamically allocate resources here if required
-@}
-
-mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s)
-@{
- // dynamically allocate resources here if required
-@}
+mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s) @{@}
+mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s) @{@}
@end example
No surprises here. We set the @code{str} member from the argument and
@subsection Automatic evaluation
-@cindex @code{hold()}
@cindex evaluation
+@cindex @code{eval()}
+@cindex @code{hold()}
When dealing with objects that are just a little more complicated than the
simple string objects we have implemented, chances are that you will want to
have some automatic simplifications or canonicalizations performed on them.
@end example
The @code{level} argument is used to limit the recursion depth of the
-evaluation. We don't have any subexpressions in the @code{mystring} class
-so we are not concerned with this. If we had, we would call the @code{eval()}
-functions of the subexpressions with @code{level - 1} as the argument if
-@code{level != 1}. The @code{hold()} member function sets a flag in the
-object that prevents further evaluation. Otherwise we might end up in an
-endless loop. When you want to return the object unmodified, use
-@code{return this->hold();}.
+evaluation. We don't have any subexpressions in the @code{mystring}
+class so we are not concerned with this. If we had, we would call the
+@code{eval()} functions of the subexpressions with @code{level - 1} as
+the argument if @code{level != 1}. The @code{hold()} member function
+sets a flag in the object that prevents further evaluation. Otherwise
+we might end up in an endless loop. When you want to return the object
+unmodified, use @code{return this->hold();}.
Let's confirm that it works:
// -> 3*"wow"
@end example
-@subsection Other member functions
+@subsection Optional member functions
We have implemented only a small set of member functions to make the class
-work in the GiNaC framework. For a real algebraic class, there are probably
-some more functions that you will want to re-implement, such as
-@code{evalf()}, @code{series()} or @code{op()}. Have a look at @file{basic.h}
-or the header file of the class you want to make a subclass of to see
-what's there. One member function that you will most likely want to
-implement for terminal classes like the described string class is
-@code{calcchash()} that returns an @code{unsigned} hash value for the object
-which will allow GiNaC to compare and canonicalize expressions much more
-efficiently.
-
-You can, of course, also add your own new member functions. Remember,
+work in the GiNaC framework. There are two functions that are not strictly
+required but will make operations with objects of the class more efficient:
+
+@cindex @code{calchash()}
+@cindex @code{is_equal_same_type()}
+@example
+unsigned calchash() const;
+bool is_equal_same_type(const basic &other) const;
+@end example
+
+The @code{calchash()} method returns an @code{unsigned} hash value for the
+object which will allow GiNaC to compare and canonicalize expressions much
+more efficiently. You should consult the implementation of some of the built-in
+GiNaC classes for examples of hash functions. The default implementation of
+@code{calchash()} calculates a hash value out of the @code{tinfo_key} of the
+class and all subexpressions that are accessible via @code{op()}.
+
+@code{is_equal_same_type()} works like @code{compare_same_type()} but only
+tests for equality without establishing an ordering relation, which is often
+faster. The default implementation of @code{is_equal_same_type()} just calls
+@code{compare_same_type()} and tests its result for zero.
+
+@subsection Other member functions
+
+For a real algebraic class, there are probably some more functions that you
+might want to provide:
+
+@example
+bool info(unsigned inf) const;
+ex evalf(int level = 0) const;
+ex series(const relational & r, int order, unsigned options = 0) const;
+ex derivative(const symbol & s) const;
+@end example
+
+If your class stores sub-expressions (see the scalar product example in the
+previous section) you will probably want to override
+
+@cindex @code{let_op()}
+@example
+size_t nops() cont;
+ex op(size_t i) const;
+ex & let_op(size_t i);
+ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ex map(map_function & f) const;
+@end example
+
+@code{let_op()} is a variant of @code{op()} that allows write access. The
+default implementations of @code{subs()} and @code{map()} use it, so you have
+to implement either @code{let_op()}, or @code{subs()} and @code{map()}.
+
+You can, of course, also add your own new member functions. Remember
that the RTTI may be used to get information about what kinds of objects
you are dealing with (the position in the class hierarchy) and that you
can always extract the bare object from an @code{ex} by stripping the
disadvantages over these systems.
@menu
-* Advantages:: Stengths of the GiNaC approach.
+* Advantages:: Strengths of the GiNaC approach.
* Disadvantages:: Weaknesses of the GiNaC approach.
* Why C++?:: Attractiveness of C++.
@end menu
windows into GiNaC have been implemented and many more are possible: the
tiny @command{ginsh} that is part of the distribution exposes GiNaC's
types to a command line and second, as a more consistent approach, an
-interactive interface to the @acronym{Cint} C++ interpreter has been put
-together (called @acronym{GiNaC-cint}) that allows an interactive
-scripting interface consistent with the C++ language.
+interactive interface to the Cint C++ interpreter has been put together
+(called GiNaC-cint) that allows an interactive scripting interface
+consistent with the C++ language. It is available from the usual GiNaC
+FTP-site.
@item
-seemless integration: it is somewhere between difficult and impossible
+seamless integration: it is somewhere between difficult and impossible
to call CAS functions from within a program written in C++ or any other
programming language and vice versa. With GiNaC, your symbolic routines
are part of your program. You can easily call third party libraries,
portability: While the GiNaC library itself is designed to avoid any
platform dependent features (it should compile on any ANSI compliant C++
compiler), the currently used version of the CLN library (fast large
-integer and arbitrary precision arithmetics) can be compiled only on
-systems with a recently new C++ compiler from the GNU Compiler
-Collection (@acronym{GCC}).@footnote{This is because CLN uses
-PROVIDE/REQUIRE like macros to let the compiler gather all static
-initializations, which works for GNU C++ only.} GiNaC uses recent
-language features like explicit constructors, mutable members, RTTI,
-@code{dynamic_cast}s and STL, so ANSI compliance is meant literally.
-Recent @acronym{GCC} versions starting at 2.95, although itself not yet
-ANSI compliant, support all needed features.
+integer and arbitrary precision arithmetics) can only by compiled
+without hassle on systems with the C++ compiler from the GNU Compiler
+Collection (GCC).@footnote{This is because CLN uses PROVIDE/REQUIRE like
+macros to let the compiler gather all static initializations, which
+works for GNU C++ only. Feel free to contact the authors in case you
+really believe that you need to use a different compiler. We have
+occasionally used other compilers and may be able to give you advice.}
+GiNaC uses recent language features like explicit constructors, mutable
+members, RTTI, @code{dynamic_cast}s and STL, so ANSI compliance is meant
+literally. Recent GCC versions starting at 2.95.3, although itself not
+yet ANSI compliant, support all needed features.
@end itemize
@cindex reference counting
@cindex copy-on-write
@cindex garbage collection
-An expression is extremely light-weight since internally it works like a
-handle to the actual representation and really holds nothing more than a
-pointer to some other object. What this means in practice is that
-whenever you create two @code{ex} and set the second equal to the first
-no copying process is involved. Instead, the copying takes place as soon
-as you try to change the second. Consider the simple sequence of code:
+In GiNaC, there is an @emph{intrusive reference-counting} mechanism at work
+where the counter belongs to the algebraic objects derived from class
+@code{basic} but is maintained by the smart pointer class @code{ptr}, of
+which @code{ex} contains an instance. If you understood that, you can safely
+skip the rest of this passage.
+
+Expressions are extremely light-weight since internally they work like
+handles to the actual representation. They really hold nothing more
+than a pointer to some other object. What this means in practice is
+that whenever you create two @code{ex} and set the second equal to the
+first no copying process is involved. Instead, the copying takes place
+as soon as you try to change the second. Consider the simple sequence
+of code:
@example
+#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
can be:
@example
-#include <ginac/ginac.h>
-using namespace std;
-using namespace GiNaC;
-
-int main()
@{
symbol x("x"), y("y");
@example
#include <ginac/ginac.h>
-int main(void)
+int main()
@{
GiNaC::symbol x("x");
GiNaC::ex a = GiNaC::sin(x);
AC_PROG_INSTALL
AC_LANG_CPLUSPLUS
-AM_PATH_GINAC(0.7.0, [
+AM_PATH_GINAC(0.9.0, [
LIBS="$LIBS $GINACLIB_LIBS"
CPPFLAGS="$CPPFLAGS $GINACLIB_CPPFLAGS"
], AC_MSG_ERROR([need to have GiNaC installed]))
@end example
This @file{Makefile.am}, says that we are building a single executable,
-from a single sourcefile @file{simple.cpp}. Since every program
+from a single source file @file{simple.cpp}. Since every program
we are building uses GiNaC we simply added the GiNaC options
to @env{$LIBS} and @env{$CPPFLAGS}, but in other circumstances, we might
want to specify them on a per-program basis: for instance by
@item
@cite{Computer Algebra: Systems and Algorithms for Algebraic Computation},
-J.H. Davenport, Y. Siret, and E. Tournier, ISBN 0-12-204230-1, 1988,
+James H. Davenport, Yvon Siret and Evelyne Tournier, ISBN 0-12-204230-1, 1988,
Academic Press, London
@item
-@cite{The Role of gamma5 in Dimensional Regularization}, D. Kreimer, hep-ph/9401354
+@cite{Computer Algebra Systems - A Practical Guide},
+Michael J. Wester (editor), ISBN 0-471-98353-5, 1999, Wiley, Chichester
+
+@item
+@cite{The Art of Computer Programming, Vol 2: Seminumerical Algorithms},
+Donald E. Knuth, ISBN 0-201-89684-2, 1998, Addison Wesley
+
+@item
+@cite{Pi Unleashed}, J@"org Arndt and Christoph Haenel,
+ISBN 3-540-66572-2, 2001, Springer, Heidelberg
+
+@item
+@cite{The Role of gamma5 in Dimensional Regularization}, Dirk Kreimer, hep-ph/9401354
@end itemize
git://www.ginac.de / ginac.git / blobdiff