This is a tutorial that documents GiNaC @value{VERSION}, an open
framework for symbolic computation within the C++ programming language.
-Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+Copyright (C) 1999-2020 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
@title GiNaC @value{VERSION}
@subtitle An open framework for symbolic computation within the C++ programming language
@subtitle @value{UPDATED}
-@author @uref{http://www.ginac.de}
+@author @uref{https://www.ginac.de}
@page
@vskip 0pt plus 1filll
-Copyright @copyright{} 1999-2015 Johannes Gutenberg University Mainz, Germany
+Copyright @copyright{} 1999-2020 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
the development, the actual documentation is inside the sources in the
form of comments. That documentation may be parsed by one of the many
Javadoc-like documentation systems. If you fail at generating it you
-may access it from @uref{http://www.ginac.de/reference/, the GiNaC home
+may access it from @uref{https://www.ginac.de/reference/, the GiNaC home
page}. It is an invaluable resource not only for the advanced user who
wishes to extend the system (or chase bugs) but for everybody who wants
to comprehend the inner workings of GiNaC. This little tutorial on the
@section License
The GiNaC framework for symbolic computation within the C++ programming
-language is Copyright @copyright{} 1999-2015 Johannes Gutenberg
+language is Copyright @copyright{} 1999-2020 Johannes Gutenberg
University Mainz, Germany.
This program is free software; you can redistribute it and/or
and run it like this:
@example
-$ c++ hello.cc -o hello -lcln -lginac
+$ c++ hello.cc -o hello -lginac -lcln
$ ./hello
355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
@end example
@end example
Multivariate polynomials and rational functions may be expanded,
-collected and normalized (i.e. converted to a ratio of two coprime
-polynomials):
+collected, factorized, and normalized (i.e. converted to a ratio of
+two coprime polynomials):
@example
> a = x^4 + 2*x^2*y^2 + 4*x^3*y + 12*x*y^3 - 3*y^4;
4*x*y-y^2+x^2
> expand(a*b);
8*x^5*y+17*x^4*y^2+43*x^2*y^4-24*x*y^5+16*x^3*y^3+3*y^6+x^6
+> factor(%);
+(4*x*y+x^2-y^2)^2*(x^2+3*y^2)
> collect(a+b,x);
4*x^3*y-y^2-3*y^4+(12*y^3+4*y)*x+x^4+x^2*(1+2*y^2)
> collect(a+b,y);
3*y^2+x^2
@end example
+Here we have made use of the @command{ginsh}-command @code{%} to pop the
+previously evaluated element from @command{ginsh}'s internal stack.
+
You can differentiate functions and expand them as Taylor or Laurent
series in a very natural syntax (the second argument of @code{series} is
a relation defining the evaluation point, the third specifies the
-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
-previously evaluated element from @command{ginsh}'s internal stack.
-
Often, functions don't have roots in closed form. Nevertheless, it's
quite easy to compute a solution numerically, to arbitrary precision:
@uref{http://pkg-config.freedesktop.org}.
Last but not least, the CLN library
is used extensively and needs to be installed on your system.
-Please get it from @uref{http://www.ginac.de/CLN/} (it is licensed under
+Please get it from @uref{https://www.ginac.de/CLN/} (it is licensed under
the GPL) and install it prior to trying to install GiNaC. The configure
script checks if it can find it and if it cannot, it will refuse to
continue.
* Matrices:: Matrices.
* Indexed objects:: Handling indexed quantities.
* Non-commutative objects:: Algebras with non-commutative products.
-* Hash maps:: A faster alternative to std::map<>.
@end menu
as "@code{\Box}" in LaTeX code (@xref{Input/output}, for more
information about the different output formats of expressions in GiNaC).
GiNaC automatically creates proper LaTeX code for symbols having names of
-greek letters (@samp{alpha}, @samp{mu}, etc.).
+greek letters (@samp{alpha}, @samp{mu}, etc.). You can retrieve the name
+and the LaTeX name of a symbol using the respective methods:
+@cindex @code{get_name()}
+@cindex @code{get_TeX_name()}
+@example
+symbol::get_name() const;
+symbol::get_TeX_name() const;
+@end example
@cindex @code{subs()}
Symbols in GiNaC can't be assigned values. If you need to store results of
@code{to_int()/to_long()} and @code{to_double()} discard the imaginary
part of complex numbers.
+Note the signature of the above methods, you may need to apply a type
+conversion and call @code{evalf()} as shown in the following example:
+@example
+ ...
+ ex e1 = 1, e2 = sin(Pi/5);
+ cout << ex_to<numeric>(e1).to_int() << endl
+ << ex_to<numeric>(e2.evalf()).to_double() << endl;
+ ...
+@end example
@node Constants, Fundamental containers, Numbers, Basic concepts
@c node-name, next, previous, up
@example
@{
symbol x("x"), y("y");
- lst l;
- l = @{x, 2, y, x+y@};
+ lst l = @{x, 2, y, x+y@};
// now, l is a list holding the expressions 'x', '2', 'y', and 'x+y',
// in that order
...
@example
...
- lst l1, l2;
- l1 = x, 2, y, x+y;
- l2 = 2, x+y, x, y;
+ lst l1 = @{x, 2, y, x+y@};
+ lst l2 = @{2, x+y, x, y@};
l1.sort();
l2.sort();
// l1 and l2 are now equal
@example
...
- lst l3;
- l3 = x, 2, 2, 2, y, x+y, y+x;
+ lst l3 = @{x, 2, 2, 2, y, x+y, y+x@};
l3.unique(); // l3 is now @{x, 2, y, x+y@}
@}
@end example
ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
ex matrix::trace() const;
ex matrix::charpoly(const ex & lambda) const;
-unsigned matrix::rank() const;
+unsigned matrix::rank(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
-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.
+The optional @samp{algo} argument of @code{determinant()} and @code{rank()}
+functions allows to select between different algorithms for calculating the
+determinant and rank respectively. 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()} (matrix)
@cindex @code{solve()}
-Matrices may also be inverted using the @code{ex matrix::inverse()}
-method and linear systems may be solved with:
+Linear systems can be solved with:
@example
matrix matrix::solve(const matrix & vars, const matrix & rhs,
contain some of the indeterminates from @code{vars}. If the system is
overdetermined, an exception is thrown.
+@cindex @code{inverse()} (matrix)
+To invert a matrix, use the method:
+
+@example
+matrix matrix::inverse(unsigned algo=solve_algo::automatic) const;
+@end example
+
+The @samp{algo} argument is optional. If given, it must be one of
+@code{solve_algo} defined in @file{flags.h}.
@node Indexed objects, Non-commutative objects, Matrices, Basic concepts
@c node-name, next, previous, up
of the metric tensor.
-@node Non-commutative objects, Hash maps, Indexed objects, Basic concepts
+@node Non-commutative objects, Methods and functions, Indexed objects, Basic concepts
@c node-name, next, previous, up
@section Non-commutative objects
something very close to @code{dirac_gamma(mu)}, although
@code{dirac_gamma} have more efficient simplification mechanism.
@cindex @code{get_metric()}
+Also, the object created by @code{clifford_unit(mu, minkmetric())} is
+not aware about the symmetry of its metric, see the start of the previous
+paragraph. A more accurate analog of 'dirac_gamma(mu)' should be
+specifies as follows:
+
+@example
+ clifford_unit(mu, indexed(minkmetric(),sy_symm(),varidx(symbol("i"),4),varidx(symbol("j"),4)));
+@end example
+
The method @code{clifford::get_metric()} returns a metric defining this
Clifford number.
@example
ex clifford_prime(const ex & e)
- inline ex clifford_star(const ex & e) @{ return e.conjugate(); @}
- inline ex clifford_bar(const ex & e) @{ return clifford_prime(e.conjugate()); @}
+ inline ex clifford_star(const ex & e)
+ inline ex clifford_bar(const ex & e)
@end example
The automorphism of a Clifford algebra @code{clifford_prime()} simply
changes signs of all Clifford units in the expression. The reversion
-of a Clifford algebra @code{clifford_star()} coincides with the
-@code{conjugate()} method and effectively reverses the order of Clifford
+of a Clifford algebra @code{clifford_star()} reverses the order of Clifford
units in any product. Finally the main anti-automorphism
of a Clifford algebra @code{clifford_bar()} is the composition of the
previous two, i.e. it makes the reversion and changes signs of all Clifford units
@end example
-@node Hash maps, Methods and functions, Non-commutative objects, Basic concepts
-@c node-name, next, previous, up
-@section Hash Maps
-@cindex hash maps
-@cindex @code{exhashmap} (class)
-
-For your convenience, GiNaC offers the container template @code{exhashmap<T>}
-that can be used as a drop-in replacement for the STL
-@code{std::map<ex, T, ex_is_less>}, using hash tables to provide faster,
-typically constant-time, element look-up than @code{map<>}.
-
-@code{exhashmap<>} supports all @code{map<>} members and operations, with the
-following differences:
-
-@itemize @bullet
-@item
-no @code{lower_bound()} and @code{upper_bound()} methods
-@item
-no reverse iterators, no @code{rbegin()}/@code{rend()}
-@item
-no @code{operator<(exhashmap, exhashmap)}
-@item
-the comparison function object @code{key_compare} is hardcoded to
-@code{ex_is_less}
-@item
-the constructor @code{exhashmap(size_t n)} allows specifying the minimum
-initial hash table size (the actual table size after construction may be
-larger than the specified value)
-@item
-the method @code{size_t bucket_count()} returns the current size of the hash
-table
-@item
-@code{insert()} and @code{erase()} operations invalidate all iterators
-@end itemize
-
-
-@node Methods and functions, Information about expressions, Hash maps, Top
+@node Methods and functions, Information about expressions, Non-commutative objects, Top
@c node-name, next, previous, up
@chapter Methods and functions
@cindex polynomial
* Symmetrization::
* Built-in functions:: List of predefined mathematical functions.
* Multiple polylogarithms::
+* Iterated integrals::
* Complex expressions::
* Solving linear systems of equations::
* Input/output:: Input and output of expressions.
predicates to the STL:
@example
-class ex_is_less : public std::binary_function<ex, ex, bool> @{
+class ex_is_less @{
public:
bool operator()(const ex &lh, const ex &rh) const;
@};
-class ex_is_equal : public std::binary_function<ex, ex, bool> @{
+class ex_is_equal @{
public:
bool operator()(const ex &lh, const ex &rh) const;
@};
// count the number of expressions equal to '1'
unsigned num_ones = std::count_if(v.begin(), v.end(),
- std::bind2nd(ex_is_equal(), 1));
+ [](const ex& e) @{ return ex_is_equal()(e, 1); @});
@end example
The implementation of @code{ex_is_less} uses the member function
To evaluate them using floating-point arithmetic you need to call
@example
-ex ex::evalf(int level = 0) const;
+ex ex::evalf() const;
@end example
@cindex @code{Digits}
@cindex @code{ldegree()}
@cindex @code{coeff()}
-The degree and low degree of a polynomial can be obtained using the two
-methods
+The degree and low degree of a polynomial in expanded form can be obtained
+using the two methods
@example
int ex::degree(const ex & s);
int ex::ldegree(const ex & s);
@end example
-which also work reliably on non-expanded input polynomials (they even work
-on rational functions, returning the asymptotic degree). By definition, the
-degree of zero is zero. To extract a coefficient with a certain power from
-an expanded polynomial you use
+These functions even work on rational functions, returning the asymptotic
+degree. By definition, the degree of zero is zero. To extract a coefficient
+with a certain power from an expanded polynomial you use
@example
ex ex::coeff(const ex & s, int n);
These functions will first normalize the expression as described above and
then return the numerator, denominator, or both as a list, respectively.
-If you need both numerator and denominator, calling @code{numer_denom()} is
-faster than using @code{numer()} and @code{denom()} separately.
+If you need both numerator and denominator, call @code{numer_denom()}: it
+is faster than using @code{numer()} and @code{denom()} separately. And even
+more important: a separate evaluation of @code{numer()} and @code{denom()}
+may result in a spurious sign, e.g. for $x/(x^2-1)$ @code{numer()} may
+return $x$ and @code{denom()} $1-x^2$.
@subsection Converting to a polynomial or rational expression
@example
ex ex::to_polynomial(exmap & m);
-ex ex::to_polynomial(lst & l);
@end example
or
@example
ex ex::to_rational(exmap & m);
-ex ex::to_rational(lst & l);
@end example
-on the expression to be converted. The supplied @code{exmap} or @code{lst}
-will be filled 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_polynomial()} or @code{.to_rational()}, so it's
-possible to use it on multiple expressions and get consistent results.
+on the expression to be converted. The supplied @code{exmap} will be filled
+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_polynomial()} or @code{.to_rational()}, so it's possible to use
+it on multiple expressions and get consistent results.
The difference between @code{.to_polynomial()} and @code{.to_rational()}
is probably best illustrated with an example:
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;
+ exmap mp;
+ ex p = a.to_polynomial(mp);
+ cout << " = " << p << "\n with " << mp << 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;
+ exmap mr;
+ ex r = a.to_rational(mr);
+ cout << " = " << r << "\n with " << mr << endl;
// = -1/3*symbol4^(-1)*y+2*symbol4^(-1)*x
// with @{symbol4==sin(x)@}
@}
@cindex @code{zeta()}
@item @code{zetaderiv(n, x)}
@tab derivatives of Riemann's zeta function
+@item @code{iterated_integral(a, y)}
+@tab iterated integral
+@cindex @code{iterated_integral()}
+@item @code{iterated_integral(a, y, N)}
+@tab iterated integral with explicit truncation parameter
+@cindex @code{iterated_integral()}
@item @code{tgamma(x)}
@tab gamma function
@cindex @code{tgamma()}
@cindex @code{psi()}
@item @code{psi(n, x)}
@tab derivatives of psi function (polygamma functions)
+@item @code{EllipticK(x)}
+@tab complete elliptic integral of the first kind
+@cindex @code{EllipticK()}
+@item @code{EllipticE(x)}
+@tab complete elliptic integral of the second kind
+@cindex @code{EllipticE()}
@item @code{factorial(n)}
@tab factorial function @math{n!}
@cindex @code{factorial()}
the maximal expansion. For example, for the exponent GiNaC firstly expands
the argument and then the function. For the logarithm and absolute value,
GiNaC uses the opposite order: firstly expands the function and then its
-argument. Of course, a user can fine-tune this behaviour by sequential
+argument. Of course, a user can fine-tune this behavior by sequential
calls of several @code{expand()} methods with desired flags.
-@node Multiple polylogarithms, Complex expressions, Built-in functions, Methods and functions
+@node Multiple polylogarithms, Iterated integrals, Built-in functions, Methods and functions
@c node-name, next, previous, up
@subsection Multiple polylogarithms
The multiple polylogarithm is the most generic member of a family of functions,
to which others like the harmonic polylogarithm, Nielsen's generalized
polylogarithm and the multiple zeta value belong.
-Everyone of these functions can also be written as a multiple polylogarithm with specific
+Each of these functions can also be written as a multiple polylogarithm with specific
parameters. This whole family of functions is therefore often referred to simply as
multiple polylogarithms, containing @code{Li}, @code{G}, @code{H}, @code{S} and @code{zeta}.
The multiple polylogarithm itself comes in two variants: @code{Li} and @code{G}. While
@cite{Numerical Evaluation of Multiple Polylogarithms},
J.Vollinga, S.Weinzierl, hep-ph/0410259
-@node Complex expressions, Solving linear systems of equations, Multiple polylogarithms, Methods and functions
+@node Iterated integrals, Complex expressions, Multiple polylogarithms, Methods and functions
+@c node-name, next, previous, up
+@subsection Iterated integrals
+
+Multiple polylogarithms are a particular example of iterated integrals.
+An iterated integral is defined by the function @code{iterated_integral(a,y)}.
+The variable @code{y} gives the upper integration limit for the outermost integration, by convention the lower integration limit is always set to zero.
+The variable @code{a} must be a GiNaC @code{lst} containing sub-classes of @code{integration_kernel} as elements.
+The depth of the iterated integral corresponds to the number of elements of @code{a}.
+The available integrands for iterated integrals are
+(for a more detailed description the user is referred to the publications listed at the end of this section)
+@cartouche
+@multitable @columnfractions .40 .60
+@item @strong{Class} @tab @strong{Description}
+@item @code{integration_kernel()}
+@tab Base class, represents the one-form @math{dy}
+@cindex @code{integration_kernel()}
+@item @code{basic_log_kernel()}
+@tab Logarithmic one-form @math{dy/y}
+@cindex @code{basic_log_kernel()}
+@item @code{multiple_polylog_kernel(z_j)}
+@tab The one-form @math{dy/(y-z_j)}
+@cindex @code{multiple_polylog_kernel()}
+@item @code{ELi_kernel(n, m, x, y)}
+@tab The one form @math{ELi_{n;m}(x;y;q) dq/q}
+@cindex @code{ELi_kernel()}
+@item @code{Ebar_kernel(n, m, x, y)}
+@tab The one form @math{\overline{E}_{n;m}(x;y;q) dq/q}
+@cindex @code{Ebar_kernel()}
+@item @code{Kronecker_dtau_kernel(k, z_j, K, C_k)}
+@tab The one form @math{C_k K (k-1)/(2 \pi i)^k g^{(k)}(z_j,K \tau) dq/q}
+@cindex @code{Kronecker_dtau_kernel()}
+@item @code{Kronecker_dz_kernel(k, z_j, tau, K, C_k)}
+@tab The one form @math{C_k (2 \pi i)^{2-k} g^{(k-1)}(z-z_j,K \tau) dz}
+@cindex @code{Kronecker_dz_kernel()}
+@item @code{Eisenstein_kernel(k, N, a, b, K, C_k)}
+@tab The one form @math{C_k E_{k,N,a,b,K}(\tau) dq/q}
+@cindex @code{Eisenstein_kernel()}
+@item @code{Eisenstein_h_kernel(k, N, r, s, C_k)}
+@tab The one form @math{C_k h_{k,N,r,s}(\tau) dq/q}
+@cindex @code{Eisenstein_h_kernel()}
+@item @code{modular_form_kernel(k, P, C_k)}
+@tab The one form @math{C_k P dq/q}
+@cindex @code{modular_form_kernel()}
+@item @code{user_defined_kernel(f, y)}
+@tab The one form @math{f(y) dy}
+@cindex @code{user_defined_kernel()}
+@end multitable
+@end cartouche
+All parameters are assumed to be such that all integration kernels have a convergent Laurent expansion
+around zero with at most a simple pole at zero.
+The iterated integral may also be called with an optional third parameter
+@code{iterated_integral(a,y,N_trunc)}, in which case the numerical evaluation will truncate the series
+expansion at order @code{N_trunc}.
+
+The classes @code{Eisenstein_kernel()}, @code{Eisenstein_h_kernel()} and @code{modular_form_kernel()}
+provide a method @code{q_expansion_modular_form(q, order)}, which can used to obtain the q-expansion
+of @math{E_{k,N,a,b,K}(\tau)}, @math{h_{k,N,r,s}(\tau)} or @math{P} to the specified order.
+
+Useful publications:
+
+@cite{Numerical evaluation of iterated integrals related to elliptic Feynman integrals},
+M.Walden, S.Weinzierl, arXiv:2010.05271
+
+@node Complex expressions, Solving linear systems of equations, Iterated integrals, Methods and functions
@c node-name, next, previous, up
@section Complex expressions
@c
@example
@{
symbol a("a"), b("b"), x("x"), y("y");
- lst eqns, vars;
- eqns = a*x+b*y==3, x-y==b;
- vars = x, y;
+ lst eqns = @{a*x+b*y==3, x-y==b@};
+ lst vars = @{x, y@};
cout << lsolve(eqns, vars) << endl;
// -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
@end example
@cindex ginac-excompiler
@code{compile_ex} uses the shell script @code{ginac-excompiler} to start the C
compiler and produce the object files. This shell script comes with GiNaC and
-will be installed together with GiNaC in the configured @code{$PREFIX/bin}
-directory. You can also export additional compiler flags via the $CXXFLAGS
-variable:
+will be installed together with GiNaC in the configured @code{$LIBEXECDIR}
+(typically @code{$PREFIX/libexec} or @code{$PREFIX/lib/ginac}). You can also
+export additional compiler flags via the @env{$CXXFLAGS} variable:
@example
setenv("CXXFLAGS", "-O3 -fomit-frame-pointer -ffast-math", 1);
@example
#include <fstream>
-using namespace std;
#include <ginac/ginac.h>
+using namespace std;
using namespace GiNaC;
int main()
@example
// ...
- ofstream out("foobar.gar");
+ ofstream out("foobar.gar", ios::binary);
out << a;
out.close();
// ...
@end example
The file @file{foobar.gar} contains all information that is needed to
-reconstruct the expressions @code{foo} and @code{bar}.
+reconstruct the expressions @code{foo} and @code{bar}. The flag
+@code{ios::binary} prevents locales setting of your OS tampers the
+archive file structure.
@cindex @command{viewgar}
The tool @command{viewgar} that comes with GiNaC can be used to view
@example
// ...
archive a2;
- ifstream in("foobar.gar");
+ ifstream in("foobar.gar", ios::binary);
in >> a2;
// ...
@end example
@example
// ...
- lst syms;
- syms = x, y;
+ lst syms = @{x, y@};
ex ex1 = a2.unarchive_ex(syms, "foo");
ex ex2 = a2.unarchive_ex(syms, "the second one");
@example
#include <iostream>
-using namespace std;
-
#include <ginac/ginac.h>
+using namespace std;
using namespace GiNaC;
struct sprod_s @{
#include <iostream>
#include <string>
#include <stdexcept>
-using namespace std;
-
#include <ginac/ginac.h>
+using namespace std;
using namespace GiNaC;
@end example
@example
bool info(unsigned inf) const override;
-ex evalf(int level = 0) const override;
+ex evalf() const override;
ex series(const relational & r, int order, unsigned options = 0) const override;
ex derivative(const symbol & s) const override;
@end example
git://www.ginac.de / ginac.git / blobdiff