X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=doc%2Ftutorial%2Fginac.texi;h=4bc6410259ac11e1e0404d318fbd67af566bebb3;hb=c97593c0903a847e84c7ae4aad182bafb6430ffd;hp=360acb05a0e98b1a38f8923f7b43920aa266cabc;hpb=1994570795185f89cc848843cc53664740d20fa0;p=ginac.git diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 360acb05..4bc64102 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -24,7 +24,7 @@ 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-2024 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 @@ -48,11 +48,11 @@ notice identical to this one. @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-2024 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 @@ -126,7 +126,7 @@ hand-made documentation like this one is difficult to keep in sync with 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 @@ -135,7 +135,7 @@ the near future. @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-2024 Johannes Gutenberg University Mainz, Germany. This program is free software; you can redistribute it and/or @@ -201,7 +201,7 @@ Assuming the file is called @file{hello.cc}, on our system we can compile 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 @@ -372,8 +372,8 @@ lambda^2-3*lambda+11 @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; @@ -382,6 +382,8 @@ polynomials): 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); @@ -390,6 +392,9 @@ polynomials): 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 @@ -414,9 +419,6 @@ x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x -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: @@ -476,7 +478,7 @@ installation. 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 GCC for development +ISO standard @cite{ISO/IEC 14882:2011(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. The pkg-config utility is @@ -484,10 +486,10 @@ required for the configuration, it can be downloaded from @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{ftp://ftpthep.physik.uni-mainz.de/pub/gnu/} -(it is covered by 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. +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. @node Configuration, Building GiNaC, Prerequisites, Installation @@ -717,7 +719,6 @@ meta-class for storing all mathematical objects. * 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 @@ -838,8 +839,8 @@ some immediate simplifications. 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; +ex ex::eval() const; +ex basic::eval() const; @end example but unless you are extending GiNaC with your own classes or functions, there @@ -1146,7 +1147,14 @@ This creates a symbol that is printed as "@code{x}" in normal output, but 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 @@ -1186,7 +1194,7 @@ 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 +the documentation of that library. @xref{Top,,, cln, The CLN Manual}, 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 @@ -1526,6 +1534,15 @@ 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. +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(e1).to_int() << endl + << ex_to(e2.evalf()).to_double() << endl; + ... +@end example @node Constants, Fundamental containers, Numbers, Basic concepts @c node-name, next, previous, up @@ -1647,30 +1664,17 @@ packages, but are sometimes used to supply a variable number of arguments of the same type to GiNaC methods such as @code{subs()} and some @code{matrix} constructors, so you should have a basic understanding of them. -Lists can be constructed by assigning a comma-separated sequence of -expressions: +Lists can be constructed from an initializer list of expressions: @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 ... @end example -There are also constructors that allow direct creation of lists of up to -16 expressions, which is often more convenient but slightly less efficient: - -@example - ... - // This produces the same list 'l' as above: - // lst l(x, 2, y, x+y); - // lst l = lst(x, 2, y, x+y); - ... -@end example - Use the @code{nops()} method to determine the size (number of expressions) of a list and the @code{op()} method or the @code{[]} operator to access individual elements: @@ -1772,9 +1776,8 @@ You can bring the elements of a list into a canonical order with @code{sort()}: @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 @@ -1786,8 +1789,7 @@ elements with @code{unique()}: @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 @@ -1864,13 +1866,31 @@ substitutions. They are also used as arguments to the @code{ex::series} method, where the left hand side of the relation specifies the variable to expand in and the right hand side the expansion point. They can also be used for creating systems of equations that are to be solved for -unknown variables. But the most common usage of objects of this class +unknown variables. + +But the most common usage of objects of this class is rather inconspicuous in statements of the form @code{if (expand(pow(a+b,2))==a*a+2*a*b+b*b) @{...@}}. Here, an implicit conversion from @code{relational} to @code{bool} takes place. Note, however, that @code{==} here does not perform any simplifications, hence @code{expand()} must be called explicitly. +Simplifications of +relationals may be more efficient if preceded by a call to +@example +ex relational::canonical() const +@end example +which returns an equivalent relation with the zero +right-hand side. For example: +@example +possymbol p("p"); +relational rel = (p >= (p*p-1)/p); +if (ex_to(rel.canonical().normal())) + cout << "correct inequality" << endl; +@end example +However, a user shall not expect that any inequality can be fully +resolved by GiNaC. + @node Integrals, Matrices, Relations, Basic concepts @c node-name, next, previous, up @section Integrals @@ -1952,9 +1972,17 @@ matrix::matrix(unsigned r, unsigned c); creates a matrix with @samp{r} rows and @samp{c} columns with all elements set to zero. -The fastest way to create a matrix with preinitialized elements is to assign -a list of comma-separated expressions to an empty matrix (see below for an -example). But you can also specify the elements as a (flat) list with +The easiest way to create a matrix is using an initializer list of +initializer lists, all of the same size: + +@example +@{ + matrix m = @{@{1, -a@}, + @{a, 1@}@}; +@} +@end example + +You can also specify the elements as a (flat) list with @example matrix::matrix(unsigned r, unsigned c, const lst & l); @@ -1977,6 +2005,7 @@ matrices: @cindex @code{symbolic_matrix()} @example ex diag_matrix(const lst & l); +ex diag_matrix(initializer_list l); ex unit_matrix(unsigned x); ex unit_matrix(unsigned r, unsigned c); ex symbolic_matrix(unsigned r, unsigned c, const string & base_name); @@ -1984,7 +2013,7 @@ ex symbolic_matrix(unsigned r, unsigned c, const string & base_name, const string & tex_base_name); @end example -@code{diag_matrix()} constructs a diagonal matrix given the list of diagonal +@code{diag_matrix()} constructs a square diagonal matrix given the 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 @@ -2011,10 +2040,9 @@ that specify which row and column to remove: @example @{ - matrix m(3,3); - m = 11, 12, 13, - 21, 22, 23, - 31, 32, 33; + matrix m = @{@{11, 12, 13@}, + @{21, 22, 23@}, + @{31, 32, 33@}@}; cout << reduced_matrix(m, 1, 1) << endl; // -> [[11,13],[31,33]] cout << sub_matrix(m, 1, 2, 1, 2) << endl; @@ -2040,9 +2068,8 @@ Here are a couple of examples for constructing matrices: @{ symbol a("a"), b("b"); - matrix M(2, 2); - M = a, 0, - 0, b; + matrix M = @{@{a, 0@}, + @{0, b@}@}; cout << M << endl; // -> [[a,0],[0,b]] @@ -2052,13 +2079,13 @@ Here are a couple of examples for constructing matrices: cout << M2 << endl; // -> [[a,0],[0,b]] - cout << matrix(2, 2, lst(a, 0, 0, b)) << endl; + cout << matrix(2, 2, lst@{a, 0, 0, b@}) << endl; // -> [[a,0],[0,b]] - cout << lst_to_matrix(lst(lst(a, 0), lst(0, b))) << endl; + cout << lst_to_matrix(lst@{lst@{a, 0@}, lst@{0, b@}@}) << endl; // -> [[a,0],[0,b]] - cout << diag_matrix(lst(a, b)) << endl; + cout << diag_matrix(lst@{a, b@}) << endl; // -> [[a,0],[0,b]] cout << unit_matrix(3) << endl; @@ -2094,13 +2121,12 @@ and @math{C}: @example @{ - matrix A(2, 2), B(2, 2), C(2, 2); - A = 1, 2, - 3, 4; - B = -1, 0, - 2, 1; - C = 8, 4, - 2, 1; + matrix A = @{@{ 1, 2@}, + @{ 3, 4@}@}; + matrix B = @{@{-1, 0@}, + @{ 2, 1@}@}; + matrix C = @{@{ 8, 4@}, + @{ 2, 1@}@}; matrix result = A.mul(B).sub(C.mul_scalar(2)); cout << result << endl; @@ -2171,22 +2197,20 @@ computing determinants, traces, characteristic polynomials and ranks: 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, @@ -2201,6 +2225,15 @@ 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. +@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 @@ -2940,10 +2973,9 @@ and scalar products): symbol x("x"), y("y"); // A is a 2x2 matrix, X is a 2x1 vector - matrix A(2, 2), X(2, 1); - A = 1, 2, - 3, 4; - X = x, y; + matrix A = @{@{1, 2@}, + @{3, 4@}@}; + matrix X = @{@{x, y@}@}; cout << indexed(A, i, i) << endl; // -> 5 @@ -2974,7 +3006,7 @@ one form for @samp{F} and explicitly multiply it with a matrix representation 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 @@ -3033,6 +3065,8 @@ canonicalize themselves according to rules specified in the implementation of the non-commutative classes. The drawback is that to work with other than the built-in algebras you have to implement new classes yourself. Both symbols and user-defined functions can be specified as being non-commutative. +For symbols, this is done by subclassing class symbol; for functions, +by explicitly setting the return type (@pxref{Symbolic functions}). @cindex @code{return_type()} @cindex @code{return_type_tinfo()} @@ -3265,7 +3299,7 @@ QED: dirac_gamma(mu.toggle_variance()) * (dirac_slash(l, D) + m * dirac_ONE()); e = dirac_trace(e).simplify_indexed(sp); - e = e.collect(lst(l, ldotq, m)); + e = e.collect(lst@{l, ldotq, m@}); cout << e << endl; // -> (8-4*D)*l^2+(8-4*D)*ldotq+4*D*m^2 @} @@ -3335,6 +3369,15 @@ Note that the call @code{clifford_unit(mu, minkmetric())} creates 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. @@ -3357,7 +3400,7 @@ ways. For example ... idx i(symbol("i"), 4); realsymbol s("s"); - ex M = diag_matrix(lst(1, -1, 0, s)); + ex M = diag_matrix(lst@{1, -1, 0, s@}); ex e = clifford_unit(i, M); ex e0 = e.subs(i == 0); ex e1 = e.subs(i == 1); @@ -3423,11 +3466,11 @@ The previous code may be rewritten with the help of @code{lst_to_clifford()} as ... idx i(symbol("i"), 4); realsymbol s("s"); - ex M = diag_matrix(lst(1, -1, 0, s)); - ex e0 = lst_to_clifford(lst(1, 0, 0, 0), i, M); - ex e1 = lst_to_clifford(lst(0, 1, 0, 0), i, M); - ex e2 = lst_to_clifford(lst(0, 0, 1, 0), i, M); - ex e3 = lst_to_clifford(lst(0, 0, 0, 1), i, M); + ex M = diag_matrix(@{1, -1, 0, s@}); + ex e0 = lst_to_clifford(lst@{1, 0, 0, 0@}, i, M); + ex e1 = lst_to_clifford(lst@{0, 1, 0, 0@}, i, M); + ex e2 = lst_to_clifford(lst@{0, 0, 1, 0@}, i, M); + ex e3 = lst_to_clifford(lst@{0, 0, 0, 1@}, i, M); ... @} @end example @@ -3482,14 +3525,13 @@ There are several functions for (anti-)automorphisms of Clifford algebras: @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 @@ -3748,43 +3790,7 @@ example: @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} -that can be used as a drop-in replacement for the STL -@code{std::map}, 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 @@ -3838,6 +3844,7 @@ avoided. * 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. @@ -3974,8 +3981,6 @@ table: @tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I}) @item @code{rational_function} @tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)}) -@item @code{algebraic} -@tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1}) @end multitable @end cartouche @@ -4068,7 +4073,7 @@ The following example illustrates the differences between @example @{ symbol A("A"), B("B"), C("C"); - ex e = lst(lst(A, B), C); + ex e = lst@{lst@{A, B@}, C@}; std::copy(e.begin(), e.end(), std::ostream_iterator(cout, "\n")); @@ -4155,12 +4160,12 @@ provides two functors that can be supplied as proper binary comparison predicates to the STL: @example -class ex_is_less : public std::binary_function @{ +class ex_is_less @{ public: bool operator()(const ex &lh, const ex &rh) const; @}; -class ex_is_equal : public std::binary_function @{ +class ex_is_equal @{ public: bool operator()(const ex &lh, const ex &rh) const; @}; @@ -4188,7 +4193,7 @@ 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)); + [](const ex& e) @{ return ex_is_equal()(e, 1); @}); @end example The implementation of @code{ex_is_less} uses the member function @@ -4211,7 +4216,7 @@ 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; +ex ex::evalf() const; @end example @cindex @code{Digits} @@ -4268,13 +4273,13 @@ In the first form, @code{subs()} accepts a relational of the form // -> 73 ex e2 = x*y + x; - cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl; + cout << "e2(-2, 4) = " << e2.subs(lst@{x == -2, y == 4@}) << endl; // -> -10 @} @end example If you specify multiple substitutions, they are performed in parallel, so e.g. -@code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}. +@code{subs(lst@{x == y, y == x@})} exchanges @samp{x} and @samp{y}. The second form of @code{subs()} takes an @code{exmap} object which is a pair associative container that maps expressions to expressions (currently @@ -4305,7 +4310,7 @@ contain the same number of elements). Using this form, you would write symbol x("x"), y("y"); ex e2 = x*y + x; - cout << "e2(-2, 4) = " << e2.subs(lst(x, y), lst(-2, 4)) << endl; + cout << "e2(-2, 4) = " << e2.subs(lst@{x, y@}, lst@{-2, 4@}) << endl; @} @end example @@ -4868,7 +4873,7 @@ 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 +@uref{https://condor.depaul.edu/dmumaugh/OOT/Design-Principles/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 @@ -5025,7 +5030,7 @@ one variable, the variables are given as a list. @example (x*y*sin(y)).is_polynomial(x) // Returns true. -(x*y*sin(y)).is_polynomial(lst(x,y)) // Returns false. +(x*y*sin(y)).is_polynomial(lst@{x,y@}) // Returns false. @end example @subsection Expanding and collecting @@ -5116,18 +5121,17 @@ a*(2*x*y+y^2+x^2) @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); @@ -5352,7 +5356,7 @@ int main() Square-free decomposition is available in GiNaC: @example -ex sqrfree(const ex & a, const lst & l = lst()); +ex sqrfree(const ex & a, const lst & l = lst@{@}); @end example 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 @@ -5362,10 +5366,10 @@ some care with subsequent processing of the result: symbol x("x"), y("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; + cout << sqrfree(BiVarPol, lst@{x,y@}) << endl; // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x) - cout << sqrfree(BiVarPol, lst(y,x)) << endl; + cout << sqrfree(BiVarPol, lst@{y,x@}) << endl; // -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x) cout << sqrfree(BiVarPol) << endl; @@ -5375,6 +5379,27 @@ some care with subsequent processing of the result: Note also, how factors with the same exponents are not fully factorized with this method. +@subsection Square-free partial fraction decomposition +@cindex square-free partial fraction decomposition +@cindex partial fraction decomposition +@cindex @code{sqrfree_parfrac()} + +GiNaC also supports square-free partial fraction decomposition of +rational functions: +@example +ex sqrfree_parfrac(const ex & a, const symbol & x); +@end example +It is called square-free because it assumes a square-free +factorization of the input's denominator: +@example + ... + symbol x("x"); + + ex rat = (x-4)/(pow(x,2)*(x+2)); + cout << sqrfree_parfrac(rat, x) << endl; + // -> -2*x^(-2)+3/2*x^(-1)-3/2*(2+x)^(-1) +@end example + @subsection Polynomial factorization @cindex factorization @cindex polynomial factorization @@ -5485,8 +5510,11 @@ ex ex::numer_denom(); 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 @@ -5500,20 +5528,18 @@ above. You do this by calling @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: @@ -5524,15 +5550,15 @@ 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)@} @} @@ -5695,7 +5721,7 @@ using namespace GiNaC; ex machin_pi(int degr) @{ symbol x; - ex pi_expansion = series_to_poly(atan(x).series(x,degr)); + ex pi_expansion = series_to_poly(atan(x).series(x==0,degr)); ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5)) -4*pi_expansion.subs(x==numeric(1,239)); return pi_approx; @@ -5774,9 +5800,9 @@ almost any kind of object (anything that is @code{subs()}able): cout << ex(indexed(A, i, j)).symmetrize() << endl; // -> 1/2*A.j.i+1/2*A.i.j - cout << ex(indexed(A, i, j, k)).antisymmetrize(lst(i, j)) << endl; + cout << ex(indexed(A, i, j, k)).antisymmetrize(lst@{i, j@}) << endl; // -> -1/2*A.j.i.k+1/2*A.i.j.k - cout << ex(lst(a, b, c)).symmetrize_cyclic(lst(a, b, c)) << endl; + cout << ex(lst@{a, b, c@}).symmetrize_cyclic(lst@{a, b, c@}) << endl; // -> 1/3*@{a,b,c@}+1/3*@{b,c,a@}+1/3*@{c,a,b@} @} @end example @@ -5887,6 +5913,12 @@ GiNaC contains the following predefined mathematical functions: @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()} @@ -5902,6 +5934,12 @@ GiNaC contains the following predefined mathematical functions: @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()} @@ -5979,10 +6017,10 @@ If both flags are set (as in the last call), then GiNaC tries to get 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 @@ -5996,7 +6034,7 @@ calls of several @code{expand()} methods with desired flags. 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 @@ -6079,7 +6117,7 @@ The functions only evaluate if the indices are integers greater than zero, excep will be interpreted as the sequence of signs for the corresponding indices @code{m} or the sign of the imaginary part for the corresponding arguments @code{a}, it must contain 1 or -1, e.g. -@code{zeta(lst(3,4), lst(-1,1))} means +@code{zeta(lst@{3,4@}, lst@{-1,1@})} means @tex $\zeta(\overline{3},4)$ @end tex @@ -6087,7 +6125,7 @@ $\zeta(\overline{3},4)$ @command{zeta(\overline@{3@},4)} @end ifnottex and -@code{G(lst(a,b), lst(-1,1), c)} means +@code{G(lst@{a,b@}, lst@{-1,1@}, c)} means @tex $G(a-0\epsilon,b+0\epsilon;c)$. @end tex @@ -6096,7 +6134,7 @@ $G(a-0\epsilon,b+0\epsilon;c)$. @end ifnottex The definition of @code{H} allows indices to be 0, 1 or -1 (in expanded notation) or equally to be any integer (in compact notation). With GiNaC expanded and compact notation can be mixed, -e.g. @code{lst(0,0,-1,0,1,0,0)}, @code{lst(0,0,-1,2,0,0)} and @code{lst(-3,2,0,0)} are equivalent as +e.g. @code{lst@{0,0,-1,0,1,0,0@}}, @code{lst@{0,0,-1,2,0,0@}} and @code{lst@{-3,2,0,0@}} are equivalent as indices. The anonymous evaluator @code{eval()} tries to reduce the functions, if possible, to the least-generic multiple polylogarithm. If all arguments are unit, it returns @code{zeta}. Arguments equal to zero get considered, too. Riemann's zeta function @code{zeta} (with depth one) @@ -6165,7 +6203,71 @@ J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001) @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 @@ -6237,9 +6339,8 @@ 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, 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 @@ -6505,7 +6606,7 @@ to map input (sub)strings to arbitrary expressions: table["x"] = x+log(y)+1; parser reader(table); ex e = reader("5*x^3 - x^2"); - // e = 5*(x+log(y)+1)^3 + (x+log(y)+1)^2 + // e = 5*(x+log(y)+1)^3 - (x+log(y)+1)^2 @} @end example @@ -6518,8 +6619,8 @@ with @code{get_syms()} method: parser reader; ex e = reader("2*x+sin(y)"); symtab table = reader.get_syms(); - symbol x = reader["x"]; - symbol y = reader["y"]; + symbol x = ex_to(table["x"]); + symbol y = ex_to(table["y"]); @} @end example @@ -6691,8 +6792,14 @@ ones supplied to @code{compile_ex} should appear in the expression. @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. +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); +compile_ex(...); +@end example @subsection Archiving @cindex @code{archive} (class) @@ -6705,8 +6812,8 @@ expression a unique name: @example #include -using namespace std; #include +using namespace std; using namespace GiNaC; int main() @@ -6726,14 +6833,16 @@ The archive can then be written to a file: @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 @@ -6751,7 +6860,7 @@ read in again: @example // ... archive a2; - ifstream in("foobar.gar"); + ifstream in("foobar.gar", ios::binary); in >> a2; // ... @end example @@ -6760,8 +6869,7 @@ And the stored expressions can be retrieved by their name: @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"); @@ -7621,9 +7729,8 @@ product in a C++ @code{struct}: @example #include -using namespace std; - #include +using namespace std; using namespace GiNaC; struct sprod_s @{ @@ -8015,9 +8122,8 @@ as follows: #include #include #include -using namespace std; - #include +using namespace std; using namespace GiNaC; @end example @@ -8240,11 +8346,11 @@ class mystring : public basic @{ ... public: - ex eval(int level = 0) const; + ex eval() const override; ... @}; -ex mystring::eval(int level) const +ex mystring::eval() const @{ string new_str; for (size_t i=0; ihold();}. +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();}. + +If our class had subobjects, we would have to evaluate them first (unless +they are all of type @code{ex}, which are automatically evaluated). We don't +have any subexpressions in the @code{mystring} class, so we are not concerned +with this. Let's confirm that it works: @@ -8292,8 +8398,8 @@ 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; +unsigned calchash() const override; +bool is_equal_same_type(const basic & other) const override; @end example The @code{calchash()} method returns an @code{unsigned} hash value for the @@ -8314,10 +8420,10 @@ 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; +bool info(unsigned inf) 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 If your class stores sub-expressions (see the scalar product example in the @@ -8325,11 +8431,11 @@ 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; +size_t nops() const override; +ex op(size_t i) const override; +ex & let_op(size_t i) override; +ex subs(const lst & ls, const lst & lr, unsigned options = 0) const override; +ex map(map_function & f) const override; @end example @code{let_op()} is a variant of @code{op()} that allows write access. The @@ -8428,14 +8534,9 @@ fix bugs in a traditional system. multiple interfaces: Though real GiNaC programs have to be written in some editor, then be compiled, linked and executed, there are more ways to work with the GiNaC engine. Many people want to play with -expressions interactively, as in traditional CASs. Currently, two such -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 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. +expressions interactively, as in traditional CASs: The tiny +@command{ginsh} that comes with the distribution exposes many, but not +all, of GiNaC's types to a command line. @item seamless integration: it is somewhere between difficult and impossible @@ -8486,8 +8587,7 @@ 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. +literally. @end itemize @@ -8602,9 +8702,9 @@ inserted. But it may be useful to remember that this is not what happens. Knowing this will enable you to write much more efficient code. If you still have an uncertain feeling with copy-on-write semantics, we recommend you have a look at the -@uref{http://www.parashift.com/c++-faq-lite/, C++-FAQ lite} by -Marshall Cline. Chapter 16 covers this issue and presents an -implementation which is pretty close to the one in GiNaC. +@uref{https://isocpp.org/faq, C++-FAQ's} chapter on memory management. +It covers this issue and presents an implementation which is pretty +close to the one in GiNaC. @node Internal representation of products and sums, Package tools, Expressions are reference counted, Internal structures @@ -8703,12 +8803,12 @@ program use @footnote{If GiNaC is installed into some non-standard directory @var{prefix} one should set the @var{PKG_CONFIG_PATH} environment variable to @var{prefix}/lib/pkgconfig for this to work.} @example -g++ -o simple `pkg-config --cflags --libs ginac` simple.cpp +g++ -o simple simple.cpp `pkg-config --cflags --libs ginac` @end example This command line might expand to (for example): @example -g++ -o simple -lginac -lcln simple.cpp +g++ -o simple simple.cpp -lginac -lcln @end example Not only is the form using @command{pkg-config} easier to type, it will @@ -8908,7 +9008,7 @@ $ make install @itemize @minus{} @item -@cite{ISO/IEC 14882:1998: Programming Languages: C++} +@cite{ISO/IEC 14882:2011: Programming Languages: C++} @item @cite{CLN: A Class Library for Numbers}, @email{haible@@ilog.fr, Bruno Haible}