4.1.1 Note: Expressions and STL containers

GiNaC expressions (ex objects) have value semantics (they can be assigned, reassigned and copied like integral types) but the operator < doesn’t provide a well-defined ordering on them. In STL-speak, expressions are ‘Assignable’ but not ‘LessThanComparable’.

This implies that in order to use expressions in sorted containers such as std::map<> and std::set<> you have to supply a suitable comparison predicate. GiNaC provides such a predicate, called ex_is_less. For example, a set of expressions should be defined as std::set<ex, ex_is_less>.

Unsorted containers such as std::vector<> and std::list<> don’t pose a problem. A std::vector<ex> works as expected.

See Getting information about expressions, for more about comparing and ordering expressions.

4.6.1 Tests on numbers

Once you have declared some numbers, assigned them to expressions and done some arithmetic with them it is frequently desired to retrieve some kind of information from them like asking whether that number is integer, rational, real or complex. For those cases GiNaC provides several useful methods. (Internally, they fall back to invocations of certain CLN functions.)

As an example, let’s construct some rational number, multiply it with some multiple of its denominator and test what comes out:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

// some very important constants:
const numeric twentyone(21);
const numeric ten(10);
const numeric five(5);

int main()
{
    numeric answer = twentyone;

    answer /= five;
    cout << answer.is_integer() << endl;  // false, it's 21/5
    answer *= ten;
    cout << answer.is_integer() << endl;  // true, it's 42 now!
}

Note that the variable answer is constructed here as an integer by numeric’s copy constructor, but in an intermediate step it 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 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.

4.6.2 Numeric functions

The following functions can be applied to numeric objects and will be evaluated immediately:

Most of these functions are also available as symbolic functions that can be used in expressions (see Mathematical functions) or, like gcd(), as polynomial algorithms.

4.6.3 Converting numbers

Sometimes it is desirable to convert a numeric object back to a built-in arithmetic type (int, double, etc.). The numeric class provides a couple of methods for this purpose:

int numeric::to_int() const;
long numeric::to_long() const;
double numeric::to_double() const;
cln::cl_N numeric::to_cl_N() const;

to_int() and to_long() only work when the number they are applied on is an exact integer. Otherwise the program will halt with a message like ‘Not a 32-bit integer’. to_double() applied on a rational number will return a floating-point approximation. Both to_int()/to_long() and 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 evalf() as shown in the following example:

    ...
    ex e1 = 1, e2 = sin(Pi/5);
    cout << ex_to<numeric>(e1).to_int() << endl
         << ex_to<numeric>(e2.evalf()).to_double() << endl;
    ...

4.14.1 Indexed quantities and their indices

Indexed expressions in GiNaC are constructed of two special types of objects, index objects and indexed objects.

• Index objects are of class idx or a subclass. Every index has a value and a dimension (which is the dimension of the space the index lives in) which can both be arbitrary expressions but are usually a number or a simple symbol. In addition, indices of class varidx have a variance (they can be co- or contravariant), and indices of class spinidx have a variance and can be dotted or undotted.
• Indexed objects are of class indexed or a subclass. They contain a base expression (which is the expression being indexed), and one or more indices.

Please notice: when printing expressions, covariant indices and indices without variance are denoted ‘.i’ while contravariant indices are denoted ‘~i’. Dotted indices have a ‘*’ in front of the index value. In the following, we are going to use that notation in the text so instead of A^i_jk we will write ‘A~i.j.k’. Index dimensions are not visible in the output.

A simple example shall illustrate the concepts:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main()
{
    symbol i_sym("i"), j_sym("j");
    idx i(i_sym, 3), j(j_sym, 3);

    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)
    ...

The idx constructor takes two arguments, the index value and the index dimension. First we define two index objects, i and j, both with the numeric dimension 3. The value of the index i is the symbol i_sym (which prints as ‘i’) and the value of the index j is the symbol j_sym (which prints as ‘j’). Next we construct an expression containing one indexed object, ‘A.i.j’. It has the symbol A as its base expression and the two indices i and j.

The dimensions of indices are normally not visible in the output, but one can request them to be printed with the index_dimensions manipulator, as shown above.

Note the difference between the indices i and j which are of class idx, and the index values which are the symbols i_sym and 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:

symbol i("i"), j("j");
e = indexed(A, i, j); // ERROR: indices must be of type idx

You can have multiple indexed objects in an expression, index values can be numeric, and index dimensions symbolic:

    ...
    symbol B("B"), dim("dim");
    cout << 4 * indexed(A, i)
          + indexed(B, idx(j_sym, 4), idx(2, 3), idx(i_sym, dim)) << endl;
     // -> B.j.2.i+4*A.i
    ...

B has a 4-dimensional symbolic index ‘k’, a 3-dimensional numeric index of value 2, and a symbolic index ‘i’ with the symbolic dimension ‘dim’. Note that GiNaC doesn’t automatically notify you that the free indices of ‘A’ and ‘B’ in the sum don’t match (you have to call simplify_indexed() for that, see below).

In fact, base expressions, index values and index dimensions can be arbitrary expressions:

    ...
    cout << indexed(A+B, idx(2*i_sym+1, dim/2)) << endl;
     // -> (B+A).(1+2*i)
    ...

It’s also possible to construct nonsense like ‘Pi.sin(x)’. You will not get an error message from this but you will probably not be able to do anything useful with it.

The methods

ex idx::get_value();
ex idx::get_dim();

return the value and dimension of an idx object. If you have an index in an expression, such as returned by calling .op() on an indexed object, you can get a reference to the idx object with the function ex_to<idx>() on the expression.

There are also the methods

bool idx::is_numeric();
bool idx::is_symbolic();
bool idx::is_dim_numeric();
bool idx::is_dim_symbolic();

for checking whether the value and dimension are numeric or symbolic (non-numeric). Using the info() method of an index (see Getting information about expressions) returns information about the index value.

If you need co- and contravariant indices, use the varidx class:

    ...
    symbol mu_sym("mu"), nu_sym("nu");
    varidx mu(mu_sym, 4), nu(nu_sym, 4); // default is contravariant ~mu, ~nu
    varidx mu_co(mu_sym, 4, true);       // covariant index .mu

    cout << indexed(A, mu, nu) << endl;
     // -> A~mu~nu
    cout << indexed(A, mu_co, nu) << endl;
     // -> A.mu~nu
    cout << indexed(A, mu.toggle_variance(), nu) << endl;
     // -> A.mu~nu
    ...

A varidx is an idx with an additional flag that marks it as co- or contravariant. The default is a contravariant (upper) index, but this can be overridden by supplying a third argument to the varidx constructor. The two methods

bool varidx::is_covariant();
bool varidx::is_contravariant();

allow you to check the variance of a varidx object (use ex_to<varidx>() to get the object reference from an expression). There’s also the very useful method

ex varidx::toggle_variance();

which makes a new index with the same value and dimension but the opposite variance. By using it you only have to define the index once.

The spinidx class provides dotted and undotted variant indices, as used in the Weyl-van-der-Waerden spinor formalism:

    ...
    symbol K("K"), C_sym("C"), D_sym("D");
    spinidx C(C_sym, 2), D(D_sym);          // default is 2-dimensional,
                                            // contravariant, undotted
    spinidx C_co(C_sym, 2, true);           // covariant index
    spinidx D_dot(D_sym, 2, false, true);   // contravariant, dotted
    spinidx D_co_dot(D_sym, 2, true, true); // covariant, dotted

    cout << indexed(K, C, D) << endl;
     // -> K~C~D
    cout << indexed(K, C_co, D_dot) << endl;
     // -> K.C~*D
    cout << indexed(K, D_co_dot, D) << endl;
     // -> K.*D~D
    ...

A spinidx is a varidx with an additional flag that marks it as dotted or undotted. The default is undotted but this can be overridden by supplying a fourth argument to the spinidx constructor. The two methods

bool spinidx::is_dotted();
bool spinidx::is_undotted();

allow you to check whether or not a spinidx object is dotted (use ex_to<spinidx>() to get the object reference from an expression). Finally, the two methods

ex spinidx::toggle_dot();
ex spinidx::toggle_variance_dot();

create a new index with the same value and dimension but opposite dottedness and the same or opposite variance.

4.14.2 Substituting indices

Sometimes you will want to substitute one symbolic index with another symbolic or numeric index, for example when calculating one specific element of a tensor expression. This is done with the .subs() method, as it is done for symbols (see Substituting expressions).

You have two possibilities here. You can either substitute the whole index by another index or expression:

    ...
    ex e = indexed(A, mu_co);
    cout << e << " becomes " << e.subs(mu_co == nu) << endl;
     // -> A.mu becomes A~nu
    cout << e << " becomes " << e.subs(mu_co == varidx(0, 4)) << endl;
     // -> A.mu becomes A~0
    cout << e << " becomes " << e.subs(mu_co == 0) << endl;
     // -> A.mu becomes A.0
    ...

The third example shows that trying to replace an index with something that is not an index will substitute the index value instead.

Alternatively, you can substitute the symbol of a symbolic index by another expression:

    ...
    ex e = indexed(A, mu_co);
    cout << e << " becomes " << e.subs(mu_sym == nu_sym) << endl;
     // -> A.mu becomes A.nu
    cout << e << " becomes " << e.subs(mu_sym == 0) << endl;
     // -> A.mu becomes A.0
    ...

As you see, with the second method only the value of the index will get substituted. Its other properties, including its dimension, remain unchanged. If you want to change the dimension of an index you have to substitute the whole index by another one with the new dimension.

Finally, substituting the base expression of an indexed object works as expected:

    ...
    ex e = indexed(A, mu_co);
    cout << e << " becomes " << e.subs(A == A+B) << endl;
     // -> A.mu becomes (B+A).mu
    ...

4.14.3 Symmetries

Indexed objects can have certain symmetry properties with respect to their indices. Symmetries are specified as a tree of objects of class symmetry that is constructed with the helper functions

symmetry sy_none(...);
symmetry sy_symm(...);
symmetry sy_anti(...);
symmetry sy_cycl(...);

sy_none() stands for no symmetry, sy_symm() and sy_anti() specify fully symmetric or antisymmetric, respectively, and sy_cycl() represents a cyclic symmetry. Each of these functions accepts up to four arguments which can be either symmetry objects themselves or unsigned integer numbers that represent an index position (counting from 0). A symmetry specification that consists of only a single sy_symm(), sy_anti() or sy_cycl() with no arguments specifies the respective symmetry for all indices.

Here are some examples of symmetry definitions:

    ...
    // No symmetry:
    e = indexed(A, i, j);
    e = indexed(A, sy_none(), i, j);     // equivalent
    e = indexed(A, sy_none(0, 1), i, j); // equivalent

    // Symmetric in all three indices:
    e = indexed(A, sy_symm(), i, j, k);
    e = indexed(A, sy_symm(0, 1, 2), i, j, k); // equivalent
    e = indexed(A, sy_symm(2, 0, 1), i, j, k); // same symmetry, but yields a
                                               // different canonical order

    // Symmetric in the first two indices only:
    e = indexed(A, sy_symm(0, 1), i, j, k);
    e = indexed(A, sy_none(sy_symm(0, 1), 2), i, j, k); // equivalent

    // Antisymmetric in the first and last index only (index ranges need not
    // be contiguous):
    e = indexed(A, sy_anti(0, 2), i, j, k);
    e = indexed(A, sy_none(sy_anti(0, 2), 1), i, j, k); // equivalent

    // An example of a mixed symmetry: antisymmetric in the first two and
    // last two indices, symmetric when swapping the first and last index
    // pairs (like the Riemann curvature tensor):
    e = indexed(A, sy_symm(sy_anti(0, 1), sy_anti(2, 3)), i, j, k, l);

    // Cyclic symmetry in all three indices:
    e = indexed(A, sy_cycl(), i, j, k);
    e = indexed(A, sy_cycl(0, 1, 2), i, j, k); // equivalent

    // The following examples are invalid constructions that will throw
    // an exception at run time.

    // An index may not appear multiple times:
    e = indexed(A, sy_symm(0, 0, 1), i, j, k); // ERROR
    e = indexed(A, sy_none(sy_symm(0, 1), sy_anti(0, 2)), i, j, k); // ERROR

    // Every child of sy_symm(), sy_anti() and sy_cycl() must refer to the
    // same number of indices:
    e = indexed(A, sy_symm(sy_anti(0, 1), 2), i, j, k); // ERROR

    // And of course, you cannot specify indices which are not there:
    e = indexed(A, sy_symm(0, 1, 2, 3), i, j, k); // ERROR
    ...

If you need to specify more than four indices, you have to use the .add() method of the symmetry class. For example, to specify full symmetry in the first six indices you would write sy_symm(0, 1, 2, 3).add(4).add(5).

If an indexed object has a symmetry, GiNaC will automatically bring the indices into a canonical order which allows for some immediate simplifications:

    ...
    cout << indexed(A, sy_symm(), i, j)
          + indexed(A, sy_symm(), j, i) << endl;
     // -> 2*A.j.i
    cout << indexed(B, sy_anti(), i, j)
          + indexed(B, sy_anti(), j, i) << endl;
     // -> 0
    cout << indexed(B, sy_anti(), i, j, k)
          - indexed(B, sy_anti(), j, k, i) << endl;
     // -> 0
    ...

4.14.4 Dummy indices

GiNaC treats certain symbolic index pairs as dummy indices meaning that a summation over the index range is implied. Symbolic indices which are not dummy indices are called free indices. Numeric indices are neither dummy nor free indices.

To be recognized as a dummy index pair, the two indices must be of the same class and their value must be the same single symbol (an index like ‘2*n+1’ is never a dummy index). If the indices are of class varidx they must also be of opposite variance; if they are of class spinidx they must be both dotted or both undotted.

The method .get_free_indices() returns a vector containing the free indices of an expression. It also checks that the free indices of the terms of a sum are consistent:

{
    symbol A("A"), B("B"), C("C");

    symbol i_sym("i"), j_sym("j"), k_sym("k"), l_sym("l");
    idx i(i_sym, 3), j(j_sym, 3), k(k_sym, 3), l(l_sym, 3);

    ex e = indexed(A, i, j) * indexed(B, j, k) + indexed(C, k, l, i, l);
    cout << exprseq(e.get_free_indices()) << endl;
     // -> (.i,.k)
     // 'j' and 'l' are dummy indices

    symbol mu_sym("mu"), nu_sym("nu"), rho_sym("rho"), sigma_sym("sigma");
    varidx mu(mu_sym, 4), nu(nu_sym, 4), rho(rho_sym, 4), sigma(sigma_sym, 4);

    e = indexed(A, mu, nu) * indexed(B, nu.toggle_variance(), rho)
      + indexed(C, mu, sigma, rho, sigma.toggle_variance());
    cout << exprseq(e.get_free_indices()) << endl;
     // -> (~mu,~rho)
     // 'nu' is a dummy index, but 'sigma' is not

    e = indexed(A, mu, mu);
    cout << exprseq(e.get_free_indices()) << endl;
     // -> (~mu)
     // 'mu' is not a dummy index because it appears twice with the same
     // variance

    e = indexed(A, mu, nu) + 42;
    cout << exprseq(e.get_free_indices()) << endl; // ERROR
     // this will throw an exception:
     // "add::get_free_indices: inconsistent indices in sum"
}

A dummy index summation like a.i b~i can be expanded for indices with numeric dimensions (e.g. 3) into the explicit sum like a.1 b~1 + a.2 b~2 + a.3 b~3. This is performed by the function

    ex expand_dummy_sum(const ex & e, bool subs_idx = false);

which takes an expression e and returns the expanded sum for all dummy indices with numeric dimensions. If the parameter subs_idx is set to true then all substitutions are made by idx class indices, i.e. without variance. In this case the above sum a.i b~i will be expanded to a.1 b.1 + a.2 b.2 + a.3 b.3.

4.14.5 Simplifying indexed expressions

In addition to the few automatic simplifications that GiNaC performs on indexed expressions (such as re-ordering the indices of symmetric tensors and calculating traces and convolutions of matrices and predefined tensors) there is the method

ex ex::simplify_indexed();
ex ex::simplify_indexed(const scalar_products & sp);

that performs some more expensive operations:

• it checks the consistency of free indices in sums in the same way get_free_indices() does
• it tries to give dummy indices that appear in different terms of a sum the same name to allow simplifications like a_i*b_i-a_j*b_j=0
• it (symbolically) calculates all possible dummy index summations/contractions with the predefined tensors (this will be explained in more detail in the next section)
• it detects contractions that vanish for symmetry reasons, for example the contraction of a symmetric and a totally antisymmetric tensor
• as a special case of dummy index summation, it can replace scalar products of two tensors with a user-defined value

The last point is done with the help of the scalar_products class which is used to store scalar products with known values (this is not an arithmetic class, you just pass it to simplify_indexed()):

{
    symbol A("A"), B("B"), C("C"), i_sym("i");
    idx i(i_sym, 3);

    scalar_products sp;
    sp.add(A, B, 0); // A and B are orthogonal
    sp.add(A, C, 0); // A and C are orthogonal
    sp.add(A, A, 4); // A^2 = 4 (A has length 2)

    e = indexed(A + B, i) * indexed(A + C, i);
    cout << e << endl;
     // -> (B+A).i*(A+C).i

    cout << e.expand(expand_options::expand_indexed).simplify_indexed(sp)
         << endl;
     // -> 4+C.i*B.i
}

The scalar_products object sp acts as a storage for the scalar products added to it with the .add() method. This method takes three arguments: the two expressions of which the scalar product is taken, and the expression to replace it with.

The example above also illustrates a feature of the expand() method: if passed the expand_indexed option it will distribute indices over sums, so ‘(A+B).i’ becomes ‘A.i+B.i’.

4.14.6 Predefined tensors

Some frequently used special tensors such as the delta, epsilon and metric tensors are predefined in GiNaC. They have special properties when contracted with other tensor expressions and some of them have constant matrix representations (they will evaluate to a number when numeric indices are specified).

• Delta tensor
• General metric tensor
• Minkowski metric tensor
• Spinor metric tensor
• Epsilon tensor

{
    symbol A("A"), B("B");

    idx i(symbol("i"), 3), j(symbol("j"), 3),
        k(symbol("k"), 3), l(symbol("l"), 3);

    ex e = indexed(A, i, j) * indexed(B, k, l)
         * delta_tensor(i, k) * delta_tensor(j, l);
    cout << e.simplify_indexed() << endl;
     // -> B.i.j*A.i.j

    cout << delta_tensor(i, i) << endl;
     // -> 3
}

{
    symbol A("A");

    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);

    ex e = metric_tensor(mu, nu) * indexed(A, nu.toggle_variance(), rho);
    cout << e.simplify_indexed() << endl;
     // -> A~mu~rho

    e = delta_tensor(mu, nu.toggle_variance()) * metric_tensor(nu, rho);
    cout << e.simplify_indexed() << endl;
     // -> g~mu~rho

    e = metric_tensor(mu.toggle_variance(), nu.toggle_variance())
      * metric_tensor(nu, rho);
    cout << e.simplify_indexed() << endl;
     // -> delta.mu~rho

    e = metric_tensor(nu.toggle_variance(), rho.toggle_variance())
      * metric_tensor(mu, nu) * (delta_tensor(mu.toggle_variance(), rho)
        + indexed(A, mu.toggle_variance(), rho));
    cout << e.simplify_indexed() << endl;
     // -> 4+A.rho~rho
}

{
    varidx mu(symbol("mu"), 4);

    e = delta_tensor(varidx(0, 4), mu.toggle_variance())
      * lorentz_g(mu, varidx(0, 4));       // negative signature
    cout << e.simplify_indexed() << endl;
     // -> 1

    e = delta_tensor(varidx(0, 4), mu.toggle_variance())
      * lorentz_g(mu, varidx(0, 4), true); // positive signature
    cout << e.simplify_indexed() << endl;
     // -> -1
}

{
    symbol psi("psi");

    spinidx A(symbol("A")), B(symbol("B")), C(symbol("C"));
    ex A_co = A.toggle_variance(), B_co = B.toggle_variance();

    e = spinor_metric(A, B) * indexed(psi, B_co);
    cout << e.simplify_indexed() << endl;
     // -> psi~A

    e = spinor_metric(A, B) * indexed(psi, A_co);
    cout << e.simplify_indexed() << endl;
     // -> -psi~B

    e = spinor_metric(A_co, B_co) * indexed(psi, B);
    cout << e.simplify_indexed() << endl;
     // -> -psi.A

    e = spinor_metric(A_co, B_co) * indexed(psi, A);
    cout << e.simplify_indexed() << endl;
     // -> psi.B

    e = spinor_metric(A_co, B_co) * spinor_metric(A, B);
    cout << e.simplify_indexed() << endl;
     // -> 2

    e = spinor_metric(A_co, B_co) * spinor_metric(B, C);
    cout << e.simplify_indexed() << endl;
     // -> -delta.A~C
}

ex epsilon_tensor(const ex & i1, const ex & i2);
ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4,
               bool pos_sig = false);

{
    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
}

4.14.7 Linear algebra

The matrix class can be used with indices to do some simple linear algebra (linear combinations and products of vectors and matrices, traces and scalar products):

{
    idx i(symbol("i"), 2), j(symbol("j"), 2);
    symbol x("x"), y("y");

    // A is a 2x2 matrix, X is a 2x1 vector
    matrix A = {{1, 2},
                {3, 4}};
    matrix X = {{x, y}};

    cout << indexed(A, i, i) << endl;
     // -> 5

    ex e = indexed(A, i, j) * indexed(X, j);
    cout << e.simplify_indexed() << endl;
     // -> [[2*y+x],[4*y+3*x]].i

    e = indexed(A, i, j) * indexed(X, i) + indexed(X, j) * 2;
    cout << e.simplify_indexed() << endl;
     // -> [[3*y+3*x,6*y+2*x]].j
}

You can of course obtain the same results with the matrix::add(), matrix::mul() and matrix::trace() methods (see Matrices) but with indices you don’t have to worry about transposing matrices.

Matrix indices always start at 0 and their dimension must match the number of rows/columns of the matrix. Matrices with one row or one column are vectors and can have one or two indices (it doesn’t matter whether it’s a row or a column vector). Other matrices must have two indices.

You should be careful when using indices with variance on matrices. GiNaC doesn’t look at the variance and doesn’t know that ‘F~mu~nu’ and ‘F.mu.nu’ are different matrices. In this case you should use only one form for ‘F’ and explicitly multiply it with a matrix representation of the metric tensor.

4.15.1 Clifford algebra

Clifford algebras are supported in two flavours: Dirac gamma matrices (more physical) and generic Clifford algebras (more mathematical).

• Dirac gamma matrices
• A generic Clifford algebra

ex dirac_gamma(const ex & mu, unsigned char rl = 0);

ex dirac_ONE(unsigned char rl = 0);

ex dirac_gamma5(unsigned char rl = 0);

ex dirac_gammaL(unsigned char rl = 0);
ex dirac_gammaR(unsigned char rl = 0);

ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);

{
    ...
    symbol a("a"), b("b"), D("D");
    varidx mu(symbol("mu"), D);
    ex e = dirac_gamma(mu) * dirac_slash(a, D)
         * dirac_gamma(mu.toggle_variance());
    cout << e << endl;
     // -> gamma~mu*a\*gamma.mu
    e = e.simplify_indexed();
    cout << e << endl;
     // -> -D*a\+2*a\
    cout << e.subs(D == 4) << endl;
     // -> -2*a\
    ...
}

ex dirac_trace(const ex & e, const std::set<unsigned char> & rls,
               const ex & trONE = 4);
ex dirac_trace(const ex & e, const lst & rll, const ex & trONE = 4);
ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4);

{
    // 4 dimensions
    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
    ex e = dirac_gamma(mu) * dirac_gamma(nu) *
           dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
    cout << dirac_trace(e).simplify_indexed() << endl;
     // -> -8*eta~rho~nu
}
...
{
    // D dimensions
    symbol D("D");
    varidx mu(symbol("mu"), D), nu(symbol("nu"), D), rho(symbol("rho"), D);
    ex e = dirac_gamma(mu) * dirac_gamma(nu) *
           dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
    cout << dirac_trace(e).simplify_indexed() << endl;
     // -> 8*eta~rho~nu-4*eta~rho~nu*D
}

{
    symbol q("q"), l("l"), m("m"), ldotq("ldotq"), D("D");
    varidx mu(symbol("mu"), D), nu(symbol("nu"), D);

    scalar_products sp;
    sp.add(l, l, pow(l, 2));
    sp.add(l, q, ldotq);

    ex e = dirac_gamma(mu) *
           (dirac_slash(l, D) + dirac_slash(q, D) + m * dirac_ONE()) *    
           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});
    cout << e << endl;
     // -> (8-4*D)*l^2+(8-4*D)*ldotq+4*D*m^2
}

{
    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
    ex e = dirac_gamma(mu) * dirac_gamma(nu) + dirac_gamma(nu) * dirac_gamma(mu);
    cout << e << endl;
     // -> gamma~mu*gamma~nu+gamma~nu*gamma~mu

    e = canonicalize_clifford(e);
    cout << e << endl;
     // -> 2*ONE*eta~mu~nu
}

    ex clifford_unit(const ex & mu, const ex & metr, unsigned char rl = 0);    

    clifford_unit(mu, indexed(minkmetric(),sy_symm(),varidx(symbol("i"),4),varidx(symbol("j"),4)));

    ex e = clifford_unit(mu, indexed(M, sy_symm(), i, j));

{
    ... 
    idx i(symbol("i"), 4);
    realsymbol s("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);
    ex e2 = e.subs(i == 2);
    ex e3 = e.subs(i == 3);
    ...
}

    ex lst_to_clifford(const ex & v, const ex & mu,  const ex & metr,
                       unsigned char rl = 0);
    ex lst_to_clifford(const ex & v, const ex & e);

{
    ...
    idx i(symbol("i"), 4);
    realsymbol s("s");
    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);
  ...
}

    lst clifford_to_lst(const ex & e, const ex & c, bool algebraic = true);

    ex clifford_prime(const ex & e)
    inline ex clifford_star(const ex & e)
    inline ex clifford_bar(const ex & e)

    ex clifford_norm(const ex & e);

    ex clifford_inverse(const ex & e);

    ex remove_dirac_ONE(const ex & e);

    ex clifford_moebius_map(const ex & a, const ex & b, const ex & c,
                            const ex & d, const ex & v, const ex & G,
                            unsigned char rl = 0);
    ex clifford_moebius_map(const ex & M, const ex & v, const ex & G,
                            unsigned char rl = 0);

char clifford_max_label(const ex & e, bool ignore_ONE = false);

    \newcommand{\clifford}[1][]{}

    \newcommand{\clifford}[2][]{\ifcase #1 #2\or \tilde{#2} \or \breve{#2} \fi}

4.15.2 Color algebra

For computations in quantum chromodynamics, GiNaC implements the base elements and structure constants of the su(3) Lie algebra (color algebra). The base elements T_a are constructed by the function

ex color_T(const ex & a, unsigned char rl = 0);

which takes two arguments: the index and a representation label in the range 0 to 255 which is used to distinguish elements of different color algebras. Objects with different labels commutate with each other. The dimension of the index must be exactly 8 and it should be of class idx, not varidx.

The unity element of a color algebra is constructed by

ex color_ONE(unsigned char rl = 0);

Please notice: You must always use color_ONE() when referring to multiples of the unity element, even though it’s customary to omit it. E.g. instead of color_T(a)*(color_T(b)*indexed(X,b)+1) you have to write color_T(a)*(color_T(b)*indexed(X,b)+color_ONE()). Otherwise, GiNaC may produce incorrect results.

The functions

ex color_d(const ex & a, const ex & b, const ex & c);
ex color_f(const ex & a, const ex & b, const ex & c);

create the symmetric and antisymmetric structure constants d_abc and f_abc which satisfy {T_a, T_b} = 1/3 delta_ab + d_abc T_c and [T_a, T_b] = i f_abc T_c.

These functions evaluate to their numerical values, if you supply numeric indices to them. The index values should be in the range from 1 to 8, not from 0 to 7. This departure from usual conventions goes along better with the notations used in physical literature.

There’s an additional function

ex color_h(const ex & a, const ex & b, const ex & c);

which returns the linear combination ‘color_d(a, b, c)+I*color_f(a, b, c)’.

The function simplify_indexed() performs some simplifications on expressions containing color objects:

{
    ...
    idx a(symbol("a"), 8), b(symbol("b"), 8), c(symbol("c"), 8),
        k(symbol("k"), 8), l(symbol("l"), 8);

    e = color_d(a, b, l) * color_f(a, b, k);
    cout << e.simplify_indexed() << endl;
     // -> 0

    e = color_d(a, b, l) * color_d(a, b, k);
    cout << e.simplify_indexed() << endl;
     // -> 5/3*delta.k.l

    e = color_f(l, a, b) * color_f(a, b, k);
    cout << e.simplify_indexed() << endl;
     // -> 3*delta.k.l

    e = color_h(a, b, c) * color_h(a, b, c);
    cout << e.simplify_indexed() << endl;
     // -> -32/3

    e = color_h(a, b, c) * color_T(b) * color_T(c);
    cout << e.simplify_indexed() << endl;
     // -> -2/3*T.a

    e = color_h(a, b, c) * color_T(a) * color_T(b) * color_T(c);
    cout << e.simplify_indexed() << endl;
     // -> -8/9*ONE

    e = color_T(k) * color_T(a) * color_T(b) * color_T(k);
    cout << e.simplify_indexed() << endl;
     // -> 1/4*delta.b.a*ONE-1/6*T.a*T.b
    ...

To calculate the trace of an expression containing color objects you use one of the functions

ex color_trace(const ex & e, const std::set<unsigned char> & rls);
ex color_trace(const ex & e, const lst & rll);
ex color_trace(const ex & e, unsigned char rl = 0);

These functions take the trace over all color ‘T’ objects in the specified set rls or list rll of representation labels, or the single label rl; ‘T’s with other labels are left standing. For example:

    ...
    e = color_trace(4 * color_T(a) * color_T(b) * color_T(c));
    cout << e << endl;
     // -> -I*f.a.c.b+d.a.c.b
}

5.1.1 Checking expression types

Sometimes it’s useful to check whether a given expression is a plain number, a sum, a polynomial with integer coefficients, or of some other specific type. GiNaC provides a couple of functions for this:

bool is_a<T>(const ex & e);
bool is_exactly_a<T>(const ex & e);
bool ex::info(unsigned flag);
unsigned ex::return_type() const;
return_type_t ex::return_type_tinfo() const;

When the test made by is_a<T>() returns true, it is safe to call one of the functions ex_to<T>(), where T is one of the class names (See The class hierarchy, for a list of all classes). For example, assuming e is an ex:

{
    …
    if (is_a<numeric>(e))
        numeric n = ex_to<numeric>(e);
    …
}

is_a<T>(e) allows you to check whether the top-level object of an expression ‘e’ is an instance of the GiNaC class ‘T’ (See The class hierarchy, for a list of all classes). This is most useful, e.g., for checking whether an expression is a number, a sum, or a product:

{
    symbol x("x");
    ex e1 = 42;
    ex e2 = 4*x - 3;
    is_a<numeric>(e1);  // true
    is_a<numeric>(e2);  // false
    is_a<add>(e1);      // false
    is_a<add>(e2);      // true
    is_a<mul>(e1);      // false
    is_a<mul>(e2);      // false
}

In contrast, is_exactly_a<T>(e) allows you to check whether the top-level object of an expression ‘e’ is an instance of the GiNaC class ‘T’, not including parent classes.

The info() method is used for checking certain attributes of expressions. The possible values for the flag argument are defined in ginac/flags.h, the most important being explained in the following table:

To determine whether an expression is commutative or non-commutative and if so, with which other expressions it would commutate, you use the methods return_type() and return_type_tinfo(). See Non-commutative objects, for an explanation of these.

5.1.2 Accessing subexpressions

Many GiNaC classes, like add, mul, lst, and function, act as containers for subexpressions. For example, the subexpressions of a sum (an add object) are the individual terms, and the subexpressions of a function are the function’s arguments.

GiNaC provides several ways of accessing subexpressions. The first way is to use the two methods

size_t ex::nops();
ex ex::op(size_t i);

nops() determines the number of subexpressions (operands) contained in the expression, while op(i) returns the i-th (0..nops()-1) subexpression. In the case of a power object, op(0) will return the basis and op(1) the exponent. For indexed objects, op(0) is the base expression and op(i), i>0 are the indices.

The second way to access subexpressions is via the STL-style random-access iterator class const_iterator and the methods

const_iterator ex::begin();
const_iterator ex::end();

begin() returns an iterator referring to the first subexpression; end() returns an iterator which is one-past the last subexpression. If the expression has no subexpressions, then begin() == end(). These iterators can also be used in conjunction with non-modifying STL algorithms.

Here is an example that (non-recursively) prints the subexpressions of a given expression in three different ways:

{
    ex e = ...

    // with nops()/op()
    for (size_t i = 0; i != e.nops(); ++i)
        cout << e.op(i) << endl;

    // with iterators
    for (const_iterator i = e.begin(); i != e.end(); ++i)
        cout << *i << endl;

    // with iterators and STL copy()
    std::copy(e.begin(), e.end(), std::ostream_iterator<ex>(cout, "\n"));
}

op()/nops() and const_iterator only access an expression’s immediate children. GiNaC provides two additional iterator classes, const_preorder_iterator and const_postorder_iterator, that iterate over all objects in an expression tree, in preorder or postorder, respectively. They are STL-style forward iterators, and are created with the methods

const_preorder_iterator ex::preorder_begin();
const_preorder_iterator ex::preorder_end();
const_postorder_iterator ex::postorder_begin();
const_postorder_iterator ex::postorder_end();

The following example illustrates the differences between const_iterator, const_preorder_iterator, and const_postorder_iterator:

{
    symbol A("A"), B("B"), C("C");
    ex e = lst{lst{A, B}, C};

    std::copy(e.begin(), e.end(),
              std::ostream_iterator<ex>(cout, "\n"));
    // {A,B}
    // C

    std::copy(e.preorder_begin(), e.preorder_end(),
              std::ostream_iterator<ex>(cout, "\n"));
    // {{A,B},C}
    // {A,B}
    // A
    // B
    // C

    std::copy(e.postorder_begin(), e.postorder_end(),
              std::ostream_iterator<ex>(cout, "\n"));
    // A
    // B
    // {A,B}
    // C
    // {{A,B},C}
}

Finally, the left-hand side and right-hand side expressions of objects of class relational (and only of these) can also be accessed with the methods

ex ex::lhs();
ex ex::rhs();

5.1.3 Comparing expressions

Expressions can be compared with the usual C++ relational operators like ==, >, and < but if the expressions contain symbols, the result is usually not determinable and the result will be false, except in the case of the != operator. You should also be aware that GiNaC will only do the most trivial test for equality (subtracting both expressions), so something like (pow(x,2)+x)/x==x+1 will return false.

Actually, if you construct an expression like a == b, this will be represented by an object of the relational class (see Relations) which is not evaluated until (explicitly or implicitly) cast to a bool.

There are also two methods

bool ex::is_equal(const ex & other);
bool ex::is_zero();

for checking whether one expression is equal to another, or equal to zero, respectively. See also the method ex::is_zero_matrix(), see Matrices.

5.1.4 Ordering expressions

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 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 <, > 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 ‘a’ and ‘b’ are different expressions, and a < b yields false, then b < a doesn’t necessarily yield true.

By default, STL classes and algorithms use the < and == 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:

class ex_is_less {
public:
    bool operator()(const ex &lh, const ex &rh) const;
};

class ex_is_equal {
public:
    bool operator()(const ex &lh, const ex &rh) const;
};

For example, to define a map that maps expressions to strings you have to use

std::map<ex, std::string, ex_is_less> myMap;

Omitting the ex_is_less template parameter will introduce spurious bugs because the map operates improperly.

Other examples for the use of the functors:

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(),
                                  [](const ex& e) { return ex_is_equal()(e, 1); });

The implementation of ex_is_less uses the member function

int ex::compare(const ex & other) const;

which returns 0 if *this and other are equal, -1 if *this sorts before other, and 1 if *this sorts after other.

5.4.1 Matching expressions

The most basic application of patterns is to check whether an expression matches a given pattern. This is done by the function

bool ex::match(const ex & pattern);
bool ex::match(const ex & pattern, exmap& repls);

This function returns true when the expression matches the pattern and false if it doesn’t. If used in the second form, the actual subexpressions matched by the wildcards get returned in the associative array repls with ‘wildcard’ as a key. If match() returns false, repls remains unmodified.

The matching algorithm works as follows:

• A single wildcard matches any expression. If one wildcard appears multiple times in a pattern, it must match the same expression in all places (e.g. ‘$0’ matches anything, and ‘$0*($0+1)’ matches ‘x*(x+1)’ but not ‘x*(y+1)’).
• If the expression is not of the same class as the pattern, the match fails (i.e. a sum only matches a sum, a function only matches a function, etc.).
• If the pattern is a function, it only matches the same function (i.e. ‘sin($0)’ matches ‘sin(x)’ but doesn’t match ‘exp(x)’).
• Except for sums and products, the match fails if the number of subexpressions (nops()) is not equal to the number of subexpressions of the pattern.
• If there are no subexpressions, the expressions and the pattern must be equal (in the sense of is_equal()).
• Except for sums and products, each subexpression (op()) must match the corresponding subexpression of the pattern.

Sums (add) and products (mul) are treated in a special way to account for their commutativity and associativity:

• If the pattern contains a term or factor that is a single wildcard, this one is used as the global wildcard. If there is more than one such wildcard, one of them is chosen as the global wildcard in a random way.
• Every term/factor of the pattern, except the global wildcard, is matched against every term of the expression in sequence. If no match is found, the whole match fails. Terms that did match are not considered in further matches.
• If there are no unmatched terms left, the match succeeds. Otherwise the match fails unless there is a global wildcard in the pattern, in which case this wildcard matches the remaining terms.

In general, having more than one single wildcard as a term of a sum or a factor of a product (such as ‘a+$0+$1’) will lead to unpredictable or ambiguous results.

Here are some examples in ginsh to demonstrate how it works (the match() function in ginsh returns ‘FAIL’ if the match fails, and the list of wildcard replacements otherwise):

> match((x+y)^a,(x+y)^a);
{}
> match((x+y)^a,(x+y)^b);
FAIL
> match((x+y)^a,$1^$2);
{$1==x+y,$2==a}
> match((x+y)^a,$1^$1);
FAIL
> match((x+y)^(x+y),$1^$1);
{$1==x+y}
> match((x+y)^(x+y),$1^$2);
{$1==x+y,$2==x+y}
> match((a+b)*(a+c),($1+b)*($1+c));
{$1==a}
> match((a+b)*(a+c),(a+$1)*(a+$2));
{$1==b,$2==c}
  (Unpredictable. The result might also be [$1==c,$2==b].)
> match((a+b)*(a+c),($1+$2)*($1+$3));
  (The result is undefined. Due to the sequential nature of the algorithm
   and the re-ordering of terms in GiNaC, the match for the first factor
   may be {$1==a,$2==b} in which case the match for the second factor
   succeeds, or it may be {$1==b,$2==a} which causes the second match to
   fail.)
> match(a*(x+y)+a*z+b,a*$1+$2);
  (This is also ambiguous and may return either {$1==z,$2==a*(x+y)+b} or
   {$1=x+y,$2=a*z+b}.)
> match(a+b+c+d+e+f,c);
FAIL
> match(a+b+c+d+e+f,c+$0);
{$0==a+e+b+f+d}
> match(a+b+c+d+e+f,c+e+$0);
{$0==a+b+f+d}
> match(a+b,a+b+$0);
{$0==0}
> match(a*b^2,a^$1*b^$2);
FAIL
  (The matching is syntactic, not algebraic, and "a" doesn't match "a^$1"
   even though a==a^1.)
> match(x*atan2(x,x^2),$0*atan2($0,$0^2));
{$0==x}
> match(atan2(y,x^2),atan2(y,$0));
{$0==x^2}

5.4.2 Matching parts of expressions

A more general way to look for patterns in expressions is provided by the member function

bool ex::has(const ex & pattern);

This function checks whether a pattern is matched by an expression itself or by any of its subexpressions.

Again some examples in ginsh for illustration (in ginsh, has() returns ‘1’ for true and ‘0’ for false):

> has(x*sin(x+y+2*a),y);
1
> has(x*sin(x+y+2*a),x+y);
0
  (This is because in GiNaC, "x+y" is not a subexpression of "x+y+2*a" (which
   has the subexpressions "x", "y" and "2*a".)
> has(x*sin(x+y+2*a),x+y+$1);
1
  (But this is possible.)
> has(x*sin(2*(x+y)+2*a),x+y);
0
  (This fails because "2*(x+y)" automatically gets converted to "2*x+2*y" of
   which "x+y" is not a subexpression.)
> has(x+1,x^$1);
0
  (Although x^1==x and x^0==1, neither "x" nor "1" are actually of the form
   "x^something".)
> has(4*x^2-x+3,$1*x);
1
> has(4*x^2+x+3,$1*x);
0
  (Another possible pitfall. The first expression matches because the term
   "-x" has the form "(-1)*x" in GiNaC. To check whether a polynomial
   contains a linear term you should use the coeff() function instead.)

The method

bool ex::find(const ex & pattern, exset& found);

works a bit like has() but it doesn’t stop upon finding the first match. Instead, it appends all found matches to the specified list. If there are multiple occurrences of the same expression, it is entered only once to the list. find() returns false if no matches were found (in ginsh, it returns an empty list):

> find(1+x+x^2+x^3,x);
{x}
> find(1+x+x^2+x^3,y);
{}
> find(1+x+x^2+x^3,x^$1);
{x^3,x^2}
  (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));
{sin(y),sin(x)}

5.4.3 Substituting expressions

Probably the most useful application of patterns is to use them for substituting expressions with the subs() method. Wildcards can be used in the search patterns as well as in the replacement expressions, where they get replaced by the expressions matched by them. subs() doesn’t know anything about algebra; it performs purely syntactic substitutions.

Some examples:

> subs(a^2+b^2+(x+y)^2,$1^2==$1^3);
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);
(a+b+c)^2
> subs((a+b+c)^2,a+b+$1==x+$1);
(x+c)^2
> 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
> subs(4*x^3-2*x^2+5*x-1,x^$0==a^$0);
-1+5*x-2*a^2+4*a^3
> subs(sin(1+sin(x)),sin($1)==cos($1));
cos(1+cos(x))
> expand(subs(a*sin(x+y)^2+a*cos(x+y)^2+b,cos($1)^2==1-sin($1)^2));
a+b

The last example would be written in C++ in this way:

{
    symbol a("a"), b("b"), x("x"), y("y");
    e = a*pow(sin(x+y), 2) + a*pow(cos(x+y), 2) + b;
    e = e.subs(pow(cos(wild()), 2) == 1-pow(sin(wild()), 2));
    cout << e.expand() << endl;
     // -> a+b
}

5.4.4 The option algebraic

Both has() and subs() take an optional argument to pass them extra options. This section describes what happens if you give the former the option has_options::algebraic or the latter subs_options::algebraic. In that case the matching condition for powers and multiplications is changed in such a way that they become more intuitive. Intuition says that x*y is a part of x*y*z. If you use these options you will find that (x*y*z).has(x*y, has_options::algebraic) indeed returns true. Besides matching some of the factors of a product also powers match as often as is possible without getting negative exponents. For example (x^5*y^2*z).subs(x^2*y^2==c, subs_options::algebraic) will return x*c^2*z. This also works with negative powers: (x^(-3)*y^(-2)*z).subs(1/(x*y)==c, subs_options::algebraic) will return x^(-1)*c^2*z.

Please notice: this only works for multiplications and not for locating x+y within x+y+z.

5.7.1 Testing whether an expression is a polynomial

Testing whether an expression is a polynomial in one or more variables can be done with the method

bool ex::is_polynomial(const ex & vars) const;

In the case of more than one variable, the variables are given as a list.

(x*y*sin(y)).is_polynomial(x)         // Returns true.
(x*y*sin(y)).is_polynomial(lst{x,y})  // Returns false.

5.7.2 Expanding and collecting

A polynomial in one or more variables has many equivalent representations. Some useful ones serve a specific purpose. Consider for example the trivariate polynomial 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 (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 example, two possibilities would be (4*y + z)*x + 20*y^2 + 21*y*z + 4*z^2 and 20*y^2 + (21*z + 4*x)*y + 4*z^2 + x*z.

To bring an expression into expanded form, its method

ex ex::expand(unsigned options = 0);

may be called. In our example above, this corresponds to 4*x*y + x*z + 20*y^2 + 21*y*z + 4*z^2. Again, since the canonical form in GiNaC is not easy to guess you should be prepared to see different orderings of terms in such sums!

Another useful representation of multivariate polynomials is as a univariate polynomial in one of the variables with the coefficients being polynomials in the remaining variables. The method collect() accomplishes this task:

ex ex::collect(const ex & s, bool distributed = false);

The first argument to collect() can also be a list of objects in which case the result is either a recursively collected polynomial, or a polynomial in a distributed form with terms like c*x1^e1*...*xn^en, as specified by the distributed flag.

Note that the original polynomial needs to be in expanded form (for the variables concerned) in order for collect() to be able to find the coefficients properly.

The following ginsh transcript shows an application of collect() together with find():

> a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
d*p*sin(x)+p*sin(x)+q*d*sin(x)+q*sin(y)+d*sin(x)+q*d*sin(y)+sin(y)+d*sin(y)
+q*sin(x)+d*sin(y)*p+sin(x)+sin(y)*p
> collect(a,{p,q});
d*sin(x)+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*p
+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*q+sin(y)+d*sin(y)+sin(x)
> collect(a,find(a,sin($1)));
(1+q+d+q*d+d*p+p)*sin(y)+(1+q+d+q*d+d*p+p)*sin(x)
> collect(a,{find(a,sin($1)),p,q});
(1+(1+d)*p+d+q*(1+d))*sin(x)+(1+(1+d)*p+d+q*(1+d))*sin(y)
> collect(a,{find(a,sin($1)),d});
(1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)

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

ex collect_common_factors(const ex & e);

This function doesn’t perform a full factorization but only looks for factors which are already explicitly present:

> 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

5.7.3 Degree and coefficients

The degree and low degree of a polynomial in expanded form can be obtained using the two methods

int ex::degree(const ex & s);
int ex::ldegree(const ex & s);

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

ex ex::coeff(const ex & s, int n);

You can also obtain the leading and trailing coefficients with the methods

ex ex::lcoeff(const ex & s);
ex ex::tcoeff(const ex & s);

which are equivalent to coeff(s, degree(s)) and coeff(s, ldegree(s)), respectively.

An application is illustrated in the next example, where a multivariate polynomial is analyzed:

{
    symbol x("x"), y("y");
    ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
                 - pow(x+y,2) + 2*pow(y+2,2) - 8;
    ex Poly = PolyInp.expand();
    
    for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) {
        cout << "The x^" << i << "-coefficient is "
             << Poly.coeff(x,i) << endl;
    }
    cout << "As polynomial in y: " 
         << Poly.collect(y) << endl;
}

When run, it returns an output in the following fashion:

The x^0-coefficient is y^2+11*y
The x^1-coefficient is 5*y^2-2*y
The x^2-coefficient is -1
The x^3-coefficient is 4*y
As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y

As always, the exact output may vary between different versions of GiNaC or even from run to run since the internal canonical ordering is not within the user’s sphere of influence.

degree(), ldegree(), coeff(), lcoeff(), tcoeff() and collect() can also be used to a certain degree with non-polynomial expressions as they not only work with symbols but with constants, functions and indexed objects as well:

{
    symbol a("a"), b("b"), c("c"), x("x");
    idx i(symbol("i"), 3);

    ex e = pow(sin(x) - cos(x), 4);
    cout << e.degree(cos(x)) << endl;
     // -> 4
    cout << e.expand().coeff(sin(x), 3) << endl;
     // -> -4*cos(x)

    e = indexed(a+b, i) * indexed(b+c, i); 
    e = e.expand(expand_options::expand_indexed);
    cout << e.collect(indexed(b, i)) << endl;
     // -> a.i*c.i+(a.i+c.i)*b.i+b.i^2
}

5.7.4 Polynomial division

The two functions

ex quo(const ex & a, const ex & b, const ex & x);
ex rem(const ex & a, const ex & b, const ex & x);

compute the quotient and remainder of univariate polynomials in the variable ‘x’. The results satisfy a = b*quo(a, b, x) + rem(a, b, x).

The additional function

ex prem(const ex & a, const ex & b, const ex & x);

computes the pseudo-remainder of ‘a’ and ‘b’ which satisfies c*a = b*q + prem(a, b, x), where c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1).

Exact division of multivariate polynomials is performed by the function

bool divide(const ex & a, const ex & b, ex & q);

If ‘b’ divides ‘a’ over the rationals, this function returns true and returns the quotient in the variable q. Otherwise it returns false in which case the value of q is undefined.

5.7.5 Unit, content and primitive part

The methods

ex ex::unit(const ex & x);
ex ex::content(const ex & x);
ex ex::primpart(const ex & x);
ex ex::primpart(const ex & x, const ex & c);

return the unit part, content part, and primitive polynomial of a multivariate polynomial with respect to the variable ‘x’ (the unit part being the sign of the leading coefficient, the content part being the GCD of the coefficients, and the primitive polynomial being the input polynomial divided by the unit and content parts). The second variant of primpart() expects the previously calculated content part of the polynomial in c, which enables it to work faster in the case where the content part has already been computed. The product of unit, content, and primitive part is the original polynomial.

Additionally, the method

void ex::unitcontprim(const ex & x, ex & u, ex & c, ex & p);

computes the unit, content, and primitive parts in one go, returning them in u, c, and p, respectively.

5.7.6 GCD, LCM and resultant

The functions for polynomial greatest common divisor and least common multiple have the synopsis

ex gcd(const ex & a, const ex & b);
ex lcm(const ex & a, const ex & b);

The functions gcd() and lcm() accept two expressions a and b as arguments and return a new expression, their greatest common divisor or least common multiple, respectively. If the polynomials a and b are coprime gcd(a,b) returns 1 and lcm(a,b) returns the product of a and b. Note that all the coefficients must be rationals.

#include <ginac/ginac.h>
using namespace GiNaC;

int main()
{
    symbol x("x"), y("y"), z("z");
    ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
    ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);

    ex P_gcd = gcd(P_a, P_b);
    // x + 5*y + 4*z
    ex P_lcm = lcm(P_a, P_b);
    // 4*x*y^2 + 13*y*x*z + 20*y^3 + 81*y^2*z + 67*y*z^2 + 3*x*z^2 + 12*z^3
}

The resultant of two expressions only makes sense with polynomials. It is always computed with respect to a specific symbol within the expressions. The function has the interface

ex resultant(const ex & a, const ex & b, const ex & s);

Resultants are symmetric in a and b. The following example computes the resultant of two expressions with respect to x and y, respectively:

#include <ginac/ginac.h>
using namespace GiNaC;

int main()
{
    symbol x("x"), y("y");

    ex e1 = x+pow(y,2), e2 = 2*pow(x,3)-1; // x+y^2, 2*x^3-1
    ex r;
    
    r = resultant(e1, e2, x); 
    // -> 1+2*y^6
    r = resultant(e1, e2, y); 
    // -> 1-4*x^3+4*x^6
}

5.7.7 Square-free decomposition

Square-free decomposition is available in GiNaC:

ex sqrfree(const ex & a, const lst & l = lst{});

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:

    ...
    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;
     // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)

    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;
     // -> depending on luck, any of the above
    ...

Note also, how factors with the same exponents are not fully factorized with this method.

5.7.8 Square-free partial fraction decomposition

GiNaC also supports square-free partial fraction decomposition of rational functions:

ex sqrfree_parfrac(const ex & a, const symbol & x);

It is called square-free because it assumes a square-free factorization of the input’s denominator:

    ...
    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)

5.7.9 Polynomial factorization

Polynomials can also be fully factored with a call to the function

ex factor(const ex & a, unsigned int options = 0);

The factorization works for univariate and multivariate polynomials with rational coefficients. The following code snippet shows its capabilities:

    ...
    cout << factor(pow(x,2)-1) << endl;
     // -> (1+x)*(-1+x)
    cout << factor(expand((x-y*z)*(x-pow(y,2)-pow(z,3))*(x+y+z))) << endl;
     // -> (y+z+x)*(y*z-x)*(y^2-x+z^3)
    cout << factor(pow(x,2)-1+sin(pow(x,2)-1)) << endl;
     // -> -1+sin(-1+x^2)+x^2
    ...

The results are as expected except for the last one where no factorization seems to have been done. This is due to the default option factor_options::polynomial (equals zero) to factor(), which tells GiNaC to try a factorization only if the expression is a valid polynomial. In the shown example this is not the case, because one term is a function.

There exists a second option factor_options::all, which tells GiNaC to ignore non-polynomial parts of an expression and also to look inside function arguments. With this option the example gives:

    ...
    cout << factor(pow(x,2)-1+sin(pow(x,2)-1), factor_options::all)
         << endl;
     // -> (-1+x)*(1+x)+sin((-1+x)*(1+x))
    ...

GiNaC’s factorization functions cannot handle algebraic extensions. Therefore the following example does not factor:

    ...
    cout << factor(pow(x,2)-2) << endl;
     // -> -2+x^2  and not  (x-sqrt(2))*(x+sqrt(2))
    ...

Factorization is useful in many applications. A lot of algorithms in computer algebra depend on the ability to factor a polynomial. Of course, factorization can also be used to simplify expressions, but it is costly and applying it to complicated expressions (high degrees or many terms) may consume far too much time. So usually, looking for a GCD at strategic points in a calculation is the cheaper and more appropriate alternative.

5.8.1 The normal method

Some basic form of simplification of expressions is called for frequently. GiNaC provides the method .normal(), which converts a rational function into an equivalent rational function of the form ‘numerator/denominator’ where numerator and denominator are coprime. If the input expression is already a fraction, it just finds the GCD of numerator and denominator and cancels it, otherwise it performs fraction addition and multiplication.

.normal() can also be used on expressions which are not rational functions as it will replace all non-rational objects (like functions or non-integer powers) by temporary symbols to bring the expression to the domain of rational functions before performing the normalization, and re-substituting these symbols afterwards. This algorithm is also available as a separate method .to_rational(), described below.

This means that both expressions t1 and t2 are indeed simplified in this little code snippet:

{
    symbol x("x");
    ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
    ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
    std::cout << "t1 is " << t1.normal() << std::endl;
    std::cout << "t2 is " << t2.normal() << std::endl;
}

Of course this works for multivariate polynomials too, so the ratio of the sample-polynomials from the section about GCD and LCM above would be normalized to P_a/P_b = (4*y+z)/(y+3*z).

5.8.2 Numerator and denominator

The numerator and denominator of an expression can be obtained with

ex ex::numer();
ex ex::denom();
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, call numer_denom(): it is faster than using numer() and denom() separately. And even more important: a separate evaluation of numer() and denom() may result in a spurious sign, e.g. for $x/(x^2-1)$ numer() may return $x$ and denom() $1-x^2$.

5.8.3 Converting to a polynomial or rational expression

Some of the methods described so far only work on polynomials or rational functions. GiNaC provides a way to extend the domain of these functions to general expressions by using the temporary replacement algorithm described above. You do this by calling

ex ex::to_polynomial(exmap & m);

or

ex ex::to_rational(exmap & m);

on the expression to be converted. The supplied exmap will be filled with the generated temporary symbols and their replacement expressions in a format that can be used directly for the subs() method. It can also already contain a list of replacements from an earlier application of .to_polynomial() or .to_rational(), so it’s possible to use it on multiple expressions and get consistent results.

The difference between .to_polynomial() and .to_rational() is probably best illustrated with an example:

{
    symbol x("x"), y("y");
    ex a = 2*x/sin(x) - y/(3*sin(x));
    cout << a << 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}

    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)}
}

The following more useful example will print ‘sin(x)-cos(x)’:

{
    symbol x("x");
    ex a = pow(sin(x), 2) - pow(cos(x), 2);
    ex b = sin(x) + cos(x);
    ex q;
    exmap m;
    divide(a.to_polynomial(m), b.to_polynomial(m), q);
    cout << q.subs(m) << endl;
}

5.12.1 Overview

GiNaC contains the following predefined mathematical functions:

For functions that have a branch cut in the complex plane, GiNaC follows the conventions of C/C++ for systems that do not support a signed zero. In particular: the natural logarithm (log) and the square root (sqrt) both have their branch cuts running along the negative real axis. The asin, acos, and atanh functions all have two branch cuts starting at +/-1 and running away towards infinity along the real axis. The atan and asinh functions have two branch cuts starting at +/-i and running away towards infinity along the imaginary axis. The acosh function has one branch cut starting at +1 and running towards -infinity. These functions are continuous as the branch cut is approached coming around the finite endpoint of the cut in a counter clockwise direction.

5.12.2 Expanding functions

GiNaC knows several expansion laws for trancedent functions, e.g. exp(a+b)=exp(a) exp(b), |zw|=|z| |w| or log(cd)=log(c)+log(d) (for positive c, d ). In order to use these rules you need to call expand() method with the option expand_options::expand_transcendental. Another relevant option is expand_options::expand_function_args. Their usage and interaction can be seen from the following example:

{
	symbol x("x"),  y("y");
	ex e=exp(pow(x+y,2));
	cout << e.expand() << endl;
	// -> exp((x+y)^2)
	cout << e.expand(expand_options::expand_transcendental) << endl;
	// -> exp((x+y)^2)
	cout << e.expand(expand_options::expand_function_args) << endl;
	// -> exp(2*x*y+x^2+y^2)
	cout << e.expand(expand_options::expand_function_args
			| expand_options::expand_transcendental) << endl;
	// -> exp(y^2)*exp(2*x*y)*exp(x^2)
}

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 behavior by sequential calls of several expand() methods with desired flags.

5.12.3 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. 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 Li, G, H, S and zeta. The multiple polylogarithm itself comes in two variants: Li and G. While Li and G in principle represent the same function, the different notations are more natural to the series representation or the integral representation, respectively.

To facilitate the discussion of these functions we distinguish between indices and arguments as parameters. In the table above indices are printed as m, s, n or p, whereas arguments are printed as x, a and y.

To define a Li, H or zeta with a depth greater than one, you have to pass a GiNaC lst for the indices m and s, and in the case of Li for the argument x as well. The parameter a of G must always be a lst containing the arguments in expanded form. If G is used with a third parameter s, s must have the same length as a. It contains then the signs of the imaginary parts of the arguments. If s is not given, the signs default to +1. Note that Li and zeta are polymorphic in this respect. They can stand in for the classical polylogarithm and Riemann’s zeta function (if depth is one), as well as for the multiple polylogarithm and the multiple zeta value, respectively. Note also, that GiNaC doesn’t check whether the lsts for two parameters do have the same length. It is up to the user to ensure this, otherwise evaluating will result in undefined behavior.

The functions print in LaTeX format as \mbox{Li}_{m_1,m_2,...,m_k}(x_1,x_2,...,x_k), \mbox{S}_{n,p}(x), \mbox{H}_{m_1,m_2,...,m_k}(x) and \zeta(m_1,m_2,...,m_k) (with the dots replaced by actual parameters). If zeta is an alternating zeta sum, i.e. zeta(m,s), the indices with negative sign are printed with a line above, e.g. \zeta(5,\overline{2}). The order of indices and arguments in the GiNaC lsts and in the output is the same.

Definitions and analytical as well as numerical properties of multiple polylogarithms are too numerous to be covered here. Instead, the user is referred to the publications listed at the end of this section. The implementation in GiNaC adheres to the definitions and conventions therein, except for a few differences which will be explicitly stated in the following.

One difference is about the order of the indices and arguments. For GiNaC we adopt the convention that the indices and arguments are understood to be in the same order as in which they appear in the series representation. This means Li_{m_1,m_2,m_3}(x,1,1) = H_{m_1,m_2,m_3}(x) and Li_{2,1}(1,1) = zeta(2,1) = zeta(3), but zeta(1,2) evaluates to infinity. So in comparison to the older ones of the referenced publications the order of indices and arguments for Li is reversed.

The functions only evaluate if the indices are integers greater than zero, except for the indices s in zeta and G as well as m in H. Since s will be interpreted as the sequence of signs for the corresponding indices m or the sign of the imaginary part for the corresponding arguments a, it must contain 1 or -1, e.g. zeta(lst{3,4}, lst{-1,1}) means zeta(\overline{3},4) and G(lst{a,b}, lst{-1,1}, c) means G(a-0\epsilon,b+0\epsilon;c). The definition of 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. lst{0,0,-1,0,1,0,0}, lst{0,0,-1,2,0,0} and lst{-3,2,0,0} are equivalent as indices. The anonymous evaluator eval() tries to reduce the functions, if possible, to the least-generic multiple polylogarithm. If all arguments are unit, it returns zeta. Arguments equal to zero get considered, too. Riemann’s zeta function zeta (with depth one) evaluates also for negative integers and positive even integers. For example:

> Li({3,1},{x,1});
S(2,2,x)
> H({-3,2},1);
-zeta({3,2},{-1,-1})
> S(3,1,1);
1/90*Pi^4

It is easy to tell for a given function into which other function it can be rewritten, may it be a less-generic or a more-generic one, except for harmonic polylogarithms H with negative indices or trailing zeros (the example above gives a hint). Signs can quickly be messed up, for example. Therefore GiNaC offers a C++ function convert_H_to_Li() to deal with the upgrade of a H to a multiple polylogarithm Li (eval() already cares for the possible downgrade):

> convert_H_to_Li({0,-2,-1,3},x);
Li({3,1,3},{-x,1,-1})
> convert_H_to_Li({2,-1,0},x);
-Li({2,1},{x,-1})*log(x)+2*Li({3,1},{x,-1})+Li({2,2},{x,-1})

Every function can be numerically evaluated for arbitrary real or complex arguments. The precision is arbitrary and can be set through the global variable Digits:

> Digits=100;
100
> evalf(zeta({3,1,3,1}));
0.005229569563530960100930652283899231589890420784634635522547448972148869544...

Note that the convention for arguments on the branch cut in GiNaC as stated above is different from the one Remiddi and Vermaseren have chosen for the harmonic polylogarithm.

If a function evaluates to infinity, no exceptions are raised, but the function is returned unevaluated, e.g. zeta(1). In long expressions this helps a lot with debugging, because you can easily spot the divergencies. But on the other hand, you have to make sure for yourself, that no illegal cancellations of divergencies happen.

Useful publications:

Nested Sums, Expansion of Transcendental Functions and Multi-Scale Multi-Loop Integrals, S.Moch, P.Uwer, S.Weinzierl, hep-ph/0110083

Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754

Special Values of Multiple Polylogarithms, J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941

Numerical Evaluation of Multiple Polylogarithms, J.Vollinga, S.Weinzierl, hep-ph/0410259

5.12.4 Iterated integrals

Multiple polylogarithms are a particular example of iterated integrals. An iterated integral is defined by the function iterated_integral(a,y). The variable y gives the upper integration limit for the outermost integration, by convention the lower integration limit is always set to zero. The variable a must be a GiNaC lst containing sub-classes of integration_kernel as elements. The depth of the iterated integral corresponds to the number of elements of 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)

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 iterated_integral(a,y,N_trunc), in which case the numerical evaluation will truncate the series expansion at order N_trunc.

The classes Eisenstein_kernel(), Eisenstein_h_kernel() and modular_form_kernel() provide a method q_expansion_modular_form(q, order), which can used to obtain the q-expansion of E_{k,N,a,b,K}(\tau), h_{k,N,r,s}(\tau) or P to the specified order.

Useful publications:

Numerical evaluation of iterated integrals related to elliptic Feynman integrals, M.Walden, S.Weinzierl, arXiv:2010.05271

5.15.1 Expression output

Expressions can simply be written to any stream:

{
    symbol x("x");
    ex e = 4.5*I+pow(x,2)*3/2;
    cout << e << endl;    // prints '4.5*I+3/2*x^2'
    // ...

The default output format is identical to the 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, pow(x,2) is printed as ‘x^2’).

It is possible to print expressions in a number of different formats with a set of stream manipulators;

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);

The tree, latex and csrc formats are also available in ginsh via the print(), print_latex() and print_csrc() functions, respectively.

All manipulators affect the stream state permanently. To reset the output format to the default, use the dflt manipulator:

    // ...
    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'
    // ...

If you don’t want to affect the format of the stream you’re working with, you can output to a temporary ostringstream like this:

    // ...
    ostringstream s;
    s << latex << e;         // format of cout remains unchanged
    cout << s.str() << endl; // prints '4.5 i+\frac{3}{2} x^{2}'
    // ...

The csrc (an alias for csrc_double), csrc_float, csrc_double and 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 (csrc_cl_N uses the classes provided by the CLN library):

    // ...
    cout << "f = " << csrc_float << e << ";\n";
    cout << "d = " << csrc_double << e << ";\n";
    cout << "n = " << csrc_cl_N << e << ";\n";
    // ...

The above example will produce (note the x^2 being converted to x*x):

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"));

The tree manipulator allows dumping the internal structure of an expression for debugging purposes:

    // ...
    cout << tree << e;
}

produces

add, hash=0x0, flags=0x3, nops=2
    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.5L0i (numeric), hash=0xa40a97e0, flags=0xf
    =====

The 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 symbol constructor.

For example, the code snippet

{
    symbol x("x", "\\circ");
    ex e = lgamma(x).series(x==0,3);
    cout << latex << e << endl;
}

will print

    {(-\ln(\circ))}+{(-\gamma_E)} \circ+{(\frac{1}{12} \pi^{2})} \circ^{2}
    +\mathcal{O}(\circ^{3})

Index dimensions are normally hidden in the output. To make them visible, use the index_dimensions manipulator. The dimensions will be written in square brackets behind each index value in the default and LaTeX output formats:

{
    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'
}

If you need any fancy special output format, e.g. for interfacing GiNaC with other algebra systems or for producing code for different programming languages, you can always traverse the expression tree yourself:

static void my_print(const ex & e)
{
    if (is_a<function>(e))
        cout << ex_to<function>(e).get_name();
    else
        cout << ex_to<basic>(e).class_name();
    cout << "(";
    size_t n = e.nops();
    if (n)
        for (size_t i=0; i<n; i++) {
            my_print(e.op(i));
            if (i != n-1)
                cout << ",";
        }
    else
        cout << e;
    cout << ")";
}

int main()
{
    my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
    return 0;
}

This will produce

add(power(numeric(3),symbol(x)),mul(sin(mul(power(constant(Pi),numeric(-1)),
symbol(y))),numeric(-2)))

If you need an output format that makes it possible to accurately reconstruct an expression by feeding the output to a suitable parser or object factory, you should consider storing the expression in an archive object and reading the object properties from there. See the section on archiving for more information.

5.15.2 Expression input

GiNaC provides no way to directly read an expression from a stream because you will usually want the user to be able to enter something like ‘2*x+sin(y)’ and have the ‘x’ and ‘y’ correspond to the symbols x and y you defined in your program and there is no way to specify the desired symbols to the >> stream input operator.

Instead, GiNaC lets you read an expression from a stream or a string, specifying the mapping between the input strings and symbols to be used:

{
    symbol x, y;
    symtab table;
    table["x"] = x;
    table["y"] = y;
    parser reader(table);
    ex e = reader("2*x+sin(y)");
}

The input syntax is the same as that used by ginsh and the stream output operator <<. Matching between the input strings and expressions is given by ‘table’. The ‘table’ in this example instructs GiNaC to substitute any input substring “x” with symbol x. Likewise, the substring “y” will be replaced with symbol y. It’s also possible to map input (sub)strings to arbitrary expressions:

{
    symbol x, y;
    symtab table;
    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
}

If no mapping is specified for a particular string GiNaC will create a symbol with corresponding name. Later on you can obtain all parser generated symbols with get_syms() method:

{
    parser reader;
    ex e = reader("2*x+sin(y)");
    symtab table = reader.get_syms();
    symbol x = ex_to<symbol>(table["x"]);
    symbol y = ex_to<symbol>(table["y"]);
}

Sometimes you might want to prevent GiNaC from inserting these extra symbols (for example, you want treat an unexpected string in the input as an error).

{
	symtab table;
	table["x"] = symbol();
	parser reader(table);
	parser.strict = true;
	ex e;
	try {
		e = reader("2*x+sin(y)");
	} catch (parse_error& err) {
		cerr << err.what() << endl;
		// prints "unknown symbol "y" in the input"
	}
}

With this parser, it’s also easy to implement interactive GiNaC programs. When running the following program interactively, remember to send an EOF marker after the input, e.g. by pressing Ctrl-D on an empty line:

#include <iostream>
#include <string>
#include <stdexcept>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main()
{
	cout << "Enter an expression containing 'x': " << flush;
	parser reader;

	try {
		ex e = reader(cin);
		symtab table = reader.get_syms();
		symbol x = table.find("x") != table.end() ? 
			   ex_to<symbol>(table["x"]) : symbol("x");
		cout << "The derivative of " << e << " with respect to x is ";
		cout << e.diff(x) << "." << endl;
	} catch (exception &p) {
		cerr << p.what() << endl;
	}
}

5.15.3 Compiling expressions to C function pointers

Numerical evaluation of algebraic expressions is seamlessly integrated into GiNaC by help of the CLN library. While CLN allows for very fast arbitrary precision numerics, which is more than sufficient for most users, sometimes only the speed of built-in floating point numbers is fast enough, e.g. for Monte Carlo integration. The only viable option then is the following: print the expression in C syntax format, manually add necessary C code, compile that program and run is as a separate application. This is not only cumbersome and involves a lot of manual intervention, but it also separates the algebraic and the numerical evaluation into different execution stages.

GiNaC offers a couple of functions that help to avoid these inconveniences and problems. The functions automatically perform the printing of a GiNaC expression and the subsequent compiling of its associated C code. The created object code is then dynamically linked to the currently running program. A function pointer to the C function that performs the numerical evaluation is returned and can be used instantly. This all happens automatically, no user intervention is needed.

The following example demonstrates the use of compile_ex:

    // ...
    symbol x("x");
    ex myexpr = sin(x) / x;

    FUNCP_1P fp;
    compile_ex(myexpr, x, fp);

    cout << fp(3.2) << endl;
    // ...

The function compile_ex is called with the expression to be compiled and its only free variable x. Upon successful completion the third parameter contains a valid function pointer to the corresponding C code module. If called like in the last line only built-in double precision numerics is involved.

The function pointer has to be defined in advance. GiNaC offers three function pointer types at the moment:

    typedef double (*FUNCP_1P) (double);
    typedef double (*FUNCP_2P) (double, double);
    typedef void (*FUNCP_CUBA) (const int*, const double[], const int*, double[]);

FUNCP_2P allows for two variables in the expression. FUNCP_CUBA is the correct type to be used with the CUBA library (http://www.feynarts.de/cuba) for numerical integrations. The details for the parameters of FUNCP_CUBA are explained in the CUBA manual.

For every function pointer type there is a matching compile_ex available:

    void compile_ex(const ex& expr, const symbol& sym, FUNCP_1P& fp,
                    const std::string filename = "");
    void compile_ex(const ex& expr, const symbol& sym1, const symbol& sym2,
                    FUNCP_2P& fp, const std::string filename = "");
    void compile_ex(const lst& exprs, const lst& syms, FUNCP_CUBA& fp,
                    const std::string filename = "");

When the last parameter filename is not supplied, compile_ex will choose a unique random name for the intermediate source and object files it produces. On program termination these files will be deleted. If one wishes to keep the C code and the object files, one can supply the filename parameter. The intermediate files will use that filename and will not be deleted.

link_ex is a function that allows to dynamically link an existing object file and to make it available via a function pointer. This is useful if you have already used compile_ex on an expression and want to avoid the compilation step to be performed over and over again when you restart your program. The precondition for this is of course, that you have chosen a filename when you did call compile_ex. For every above mentioned function pointer type there exists a corresponding link_ex function:

    void link_ex(const std::string filename, FUNCP_1P& fp);
    void link_ex(const std::string filename, FUNCP_2P& fp);
    void link_ex(const std::string filename, FUNCP_CUBA& fp);

The complete filename (including the suffix .so) of the object file has to be supplied.

The function

    void unlink_ex(const std::string filename);

is supplied for the rare cases when one wishes to close the dynamically linked object files directly and have the intermediate files (only if filename has not been given) deleted. Normally one doesn’t need this function, because all the clean-up will be done automatically upon (regular) program termination.

All the described functions will throw an exception in case they cannot perform correctly, like for example when writing the file or starting the compiler fails. Since internally the same printing methods as described in section csrc printing are used, only functions and objects that are available in standard C will compile successfully (that excludes polylogarithms for example at the moment). Another precondition for success is, of course, that it must be possible to evaluate the expression numerically. No free variables despite the ones supplied to compile_ex should appear in the expression.

compile_ex uses the shell script 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 $LIBEXECDIR (typically $PREFIX/libexec or $PREFIX/lib/ginac). You can also export additional compiler flags via the $CXXFLAGS variable:

setenv("CXXFLAGS", "-O3 -fomit-frame-pointer -ffast-math", 1);
compile_ex(...);

5.15.4 Archiving

GiNaC allows creating archives of expressions which can be stored to or retrieved from files. To create an archive, you declare an object of class archive and archive expressions in it, giving each expression a unique name:

#include <fstream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main()
{
    symbol x("x"), y("y"), z("z");

    ex foo = sin(x + 2*y) + 3*z + 41;
    ex bar = foo + 1;

    archive a;
    a.archive_ex(foo, "foo");
    a.archive_ex(bar, "the second one");
    // ...

The archive can then be written to a file:

    // ...
    ofstream out("foobar.gar", ios::binary);
    out << a;
    out.close();
    // ...

The file foobar.gar contains all information that is needed to reconstruct the expressions foo and bar. The flag ios::binary prevents locales setting of your OS tampers the archive file structure.

The tool viewgar that comes with GiNaC can be used to view the contents of GiNaC archive files:

$ viewgar foobar.gar
foo = 41+sin(x+2*y)+3*z
the second one = 42+sin(x+2*y)+3*z

The point of writing archive files is of course that they can later be read in again:

    // ...
    archive a2;
    ifstream in("foobar.gar", ios::binary);
    in >> a2;
    // ...

And the stored expressions can be retrieved by their name:

    // ...
    lst syms = {x, y};

    ex ex1 = a2.unarchive_ex(syms, "foo");
    ex ex2 = a2.unarchive_ex(syms, "the second one");

    cout << ex1 << endl;              // prints "41+sin(x+2*y)+3*z"
    cout << ex2 << endl;              // prints "42+sin(x+2*y)+3*z"
    cout << ex1.subs(x == 2) << endl; // prints "41+sin(2+2*y)+3*z"
}

Note that you have to supply a list of the symbols which are to be inserted in the expressions. Symbols in archives are stored by their name only and if you don’t specify which symbols you have, unarchiving the expression will create new symbols with that name. E.g. if you hadn’t included x in the syms list above, the ex1.subs(x == 2) statement would have had no effect because the x in ex1 would have been a different symbol than the x which was defined at the beginning of the program, although both would appear as ‘x’ when printed.

You can also use the information stored in an archive object to output expressions in a format suitable for exact reconstruction. The archive and archive_node classes have a couple of member functions that let you access the stored properties:

static void my_print2(const archive_node & n)
{
    string class_name;
    n.find_string("class", class_name);
    cout << class_name << "(";

    archive_node::propinfovector p;
    n.get_properties(p);

    size_t num = p.size();
    for (size_t i=0; i<num; i++) {
        const string &name = p[i].name;
        if (name == "class")
            continue;
        cout << name << "=";

        unsigned count = p[i].count;
        if (count > 1)
            cout << "{";

        for (unsigned j=0; j<count; j++) {
            switch (p[i].type) {
                case archive_node::PTYPE_BOOL: {
                    bool x;
                    n.find_bool(name, x, j);
                    cout << (x ? "true" : "false");
                    break;
                }
                case archive_node::PTYPE_UNSIGNED: {
                    unsigned x;
                    n.find_unsigned(name, x, j);
                    cout << x;
                    break;
                }
                case archive_node::PTYPE_STRING: {
                    string x;
                    n.find_string(name, x, j);
                    cout << '\"' << x << '\"';
                    break;
                }
                case archive_node::PTYPE_NODE: {
                    const archive_node &x = n.find_ex_node(name, j);
                    my_print2(x);
                    break;
                }
            }

            if (j != count-1)
                cout << ",";
        }

        if (count > 1)
            cout << "}";

        if (i != num-1)
            cout << ",";
    }

    cout << ")";
}

int main()
{
    ex e = pow(2, x) - y;
    archive ar(e, "e");
    my_print2(ar.get_top_node(0)); cout << endl;
    return 0;
}

This will produce:

add(rest={power(basis=numeric(number="2"),exponent=symbol(name="x")),
symbol(name="y")},coeff={numeric(number="1"),numeric(number="-1")},
overall_coeff=numeric(number="0"))

Be warned, however, that the set of properties and their meaning for each class may change between GiNaC versions.

6.2.1 A minimal example

Here is an example for the implementation of a function with two arguments that is not further evaluated:

DECLARE_FUNCTION_2P(myfcn)

REGISTER_FUNCTION(myfcn, dummy())

Any code that has seen the DECLARE_FUNCTION line can use myfcn() in algebraic expressions:

{
    ...
    symbol x("x");
    ex e = 2*myfcn(42, 1+3*x) - x;
    cout << e << endl;
     // prints '2*myfcn(42,1+3*x)-x'
    ...
}

The dummy() option in the REGISTER_FUNCTION line signifies "no options". A function with no options specified merely acts as a kind of container for its arguments. It is a pure "dummy" function with no associated logic (which is, however, sometimes perfectly sufficient).

Let’s now have a look at the implementation of GiNaC’s cosine function for an example of how to make an "intelligent" function.

6.2.2 The cosine function

The GiNaC header file inifcns.h contains the line

DECLARE_FUNCTION_1P(cos)

which declares to all programs using GiNaC that there is a function ‘cos’ that takes one 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 inifcns_trans.cpp. Here is its REGISTER_FUNCTION line:

REGISTER_FUNCTION(cos, eval_func(cos_eval).
                       evalf_func(cos_evalf).
                       derivative_func(cos_deriv).
                       latex_name("\\cos"));

There are four options defined for the cosine function. One of them (latex_name) gives the function a proper name for LaTeX output; the other three indicate the C++ functions in which the "brains" of the cosine function are defined.

The eval_func() option specifies the C++ function that implements the eval() method, GiNaC’s anonymous evaluator. This function takes the same number of arguments as the associated symbolic function (one in this case) and returns the (possibly transformed or in some way simplified) symbolically evaluated function (See Automatic evaluation and canonicalization of expressions, for a description of the automatic evaluation process). If no (further) evaluation is to take place, the eval_func() function must return the original function with .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 .hold() somewhere.

The eval_func() function for the cosine looks something like this (actually, it doesn’t look like this at all, but it should give you an idea what is going on):

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();
}

This function is called every time the cosine is used in a symbolic expression:

{
    ...
    e = cos(Pi);
     // this calls cos_eval(Pi), and inserts its return value into
     // the actual expression
    cout << e << endl;
     // prints '-1'
    ...
}

In this way, cos(4*Pi) automatically becomes 1, cos(asin(a+b)) becomes sqrt(1-(a+b)^2), etc. If no reasonable symbolic transformation can be done, the unmodified function is returned with .hold().

GiNaC doesn’t automatically transform cos(2) to ‘-0.416146...’. The user has to call evalf() for that. This is implemented in a different function:

static ex cos_evalf(const ex & x)
{
    if (is_a<numeric>(x))
        return cos(ex_to<numeric>(x));
    else
        return cos(x).hold();
}

Since we are lazy we defer the problem of numeric evaluation to somebody else, in this case the cos() function for numeric objects, which in turn hands it over to the cos() function in CLN. The .hold() isn’t really needed here, but reminds us that the corresponding eval() function would require it in this place.

Differentiation will surely turn up and so we need to tell cos what its first derivative is (higher derivatives, .diff(x,3) for instance, are then handled automatically by basic::diff and ex::diff):

static ex cos_deriv(const ex & x, unsigned diff_param)
{
    return -sin(x);
}

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.

Derivatives of some functions, for example abs() and Order(), could not be evaluated through the chain rule. In such cases the full derivative may be specified as shown for Order():

static ex Order_expl_derivative(const ex & arg, const symbol & s)
{
	return Order(arg.diff(s));
}

That is, we need to supply a procedure, which returns the expression of derivative with respect to the variable s for the argument arg. This procedure need to be registered with the function through the option expl_derivative_func (see the next Subsection). In contrast, a partial derivative, e.g. as was defined for cos() above, needs to be registered through the option derivative_func.

An implementation of the series expansion is not needed for 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 cos() is). tan(), on the other hand, does have poles and may need to do Laurent expansion:

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);
}

The series() implementation of a function must return a pseries object, otherwise your code will crash.

6.2.3 Function options

GiNaC functions understand several more options which are always specified as .option(params). None of them are required, but you need to specify at least one option to REGISTER_FUNCTION(). There is a do-nothing option called dummy() which you can use to define functions without any special options.

eval_func(<C++ function>)
evalf_func(<C++ function>)
derivative_func(<C++ function>)
expl_derivative_func(<C++ function>)
series_func(<C++ function>)
conjugate_func(<C++ function>)

These specify the C++ functions that implement symbolic evaluation, numeric evaluation, partial derivatives, explicit derivative, and series expansion, respectively. They correspond to the GiNaC methods eval(), evalf(), diff() and series().

The eval_func() function needs to use .hold() if no further automatic evaluation is desired or possible.

If no 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 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.

latex_name(const string & n)

specifies the LaTeX code that represents the name of the function in LaTeX output. The default is to put the function name in an \mbox{}.

do_not_evalf_params()

This tells evalf() to not recursively evaluate the parameters of the function before calling the evalf_func().

set_return_type(unsigned return_type, const return_type_t * return_type_tinfo)

This allows you to explicitly specify the commutation properties of the function (See Non-commutative objects, for an explanation of (non)commutativity in GiNaC). For example, with an object of type return_type_t created like

return_type_t my_type = make_return_type_t<matrix>();

you can use set_return_type(return_types::noncommutative, &my_type) to make GiNaC treat your function like a matrix. By default, functions inherit the commutation properties of their first argument. The utilized template function make_return_type_t<>()

template<typename T> inline return_type_t make_return_type_t(const unsigned rl = 0)

can also be called with an argument specifying the representation label of the non-commutative function (see section on dirac gamma matrices for more details).

set_symmetry(const symmetry & s)

specifies the symmetry properties of the function with respect to its arguments. See Indexed objects, for an explanation of symmetry specifications. GiNaC will automatically rearrange the arguments of symmetric functions into a canonical order.

Sometimes you may want to have finer control over how functions are displayed in the output. For example, the abs() function prints itself as ‘abs(x)’ in the default output format, but as ‘|x|’ in LaTeX mode, and fabs(x) in C source output. This is achieved with the

print_func<C>(<C++ function>)

option which is explained in the next section.

6.2.4 Functions with a variable number of arguments

The DECLARE_FUNCTION and REGISTER_FUNCTION macros define functions with a fixed number of arguments. Sometimes, though, you may need to have a function that accepts a variable number of expressions. One way to accomplish this is to pass variable-length lists as arguments. The Li() function uses this method for multiple polylogarithms.

It is also possible to define functions that accept a different number of parameters under the same function name, such as the psi() function which can be called either as psi(z) (the digamma function) or as psi(n, z) (polygamma functions). These are actually two different functions in GiNaC that, however, have the same name. Defining such functions is not possible with the macros but requires manually fiddling with GiNaC internals. If you are interested, please consult the GiNaC source code for the psi() function (inifcns.h and inifcns_gamma.cpp).

6.3.1 Print methods for classes

The method table for a class is set up either in the definition of the class, by passing the appropriate print_func<C>() option to GINAC_IMPLEMENT_REGISTERED_CLASS_OPT() (See Adding classes, for an example), or at run-time using set_print_func<T, C>(). The latter can also be used to override existing methods dynamically.

The argument to print_func<C>() and set_print_func<T, C>() can be a member function of the class (or one of its parent classes), a static member function, or an ordinary (global) C++ function. The C template parameter specifies the appropriate print_context type for which the method should be invoked, while, in the case of set_print_func<>(), the T parameter specifies the algebraic class (for print_func<>(), the class is the one being implemented by GINAC_IMPLEMENT_REGISTERED_CLASS_OPT).

For print methods that are member functions, their first argument must be of a type convertible to a const C &, and the second argument must be an unsigned.

For static members and global functions, the first argument must be of a type convertible to a const T &, the second argument must be of a type convertible to a const C &, and the third argument must be an unsigned. A global function will, of course, not have access to private and protected members of T.

The unsigned argument of the print methods (and of ex::print() and basic::print()) is used for proper parenthesizing of the output (and by print_tree for proper indentation). It can be used for similar purposes if you write your own output formats.

The explanations given above may seem complicated, but in practice it’s really simple, as shown in the following example. Suppose that we want to display exponents in LaTeX output not as superscripts but with little upwards-pointing arrows. This can be achieved in the following way:

void my_print_power_as_latex(const power & p,
                             const print_latex & c,
                             unsigned level)
{
    // get the precedence of the 'power' class
    unsigned power_prec = p.precedence();

    // if the parent operator has the same or a higher precedence
    // we need parentheses around the power
    if (level >= power_prec)
        c.s << '(';

    // print the basis and exponent, each enclosed in braces, and
    // separated by an uparrow
    c.s << '{';
    p.op(0).print(c, power_prec);
    c.s << "}\\uparrow{";
    p.op(1).print(c, power_prec);
    c.s << '}';

    // don't forget the closing parenthesis
    if (level >= power_prec)
        c.s << ')';
}
                                                                                
int main()
{
    // a sample expression
    symbol x("x"), y("y");
    ex e = -3*pow(x, 3)*pow(y, -2) + pow(x+y, 2) - 1;

    // switch to LaTeX mode
    cout << latex;

    // this prints "-1+{(y+x)}^{2}-3 \frac{x^{3}}{y^{2}}"
    cout << e << endl;

    // now we replace the method for the LaTeX output of powers with
    // our own one
    set_print_func<power, print_latex>(my_print_power_as_latex);

    // this prints "-1+{{(y+x)}}\uparrow{2}-3 \frac{{x}\uparrow{3}}{{y}
    //              \uparrow{2}}"
    cout << e << endl;
}

Some notes:

• The first argument of my_print_power_as_latex could also have been a const basic &, the second one a const print_context &.
• The above code depends on mul objects converting their operands to power objects for the purpose of printing.
• The output of products including negative powers as fractions is also controlled by the mul class.
• The power/print_latex method provided by GiNaC prints square roots using \sqrt, but the above code doesn’t.

It’s not possible to restore a method table entry to its previous or default value. Once you have called set_print_func(), you can only override it with another call to set_print_func(), but you can’t easily go back to the default behavior again (you can, of course, dig around in the GiNaC sources, find the method that is installed at startup (power::do_print_latex in this case), and set_print_func that one; that is, after you circumvent the C++ member access control…).

6.3.2 Print methods for functions

Symbolic functions employ a print method dispatch mechanism similar to the one used for classes. The methods are specified with print_func<C>() function options. If you don’t specify any special print methods, the function will be printed with its name (or LaTeX name, if supplied), followed by a comma-separated list of arguments enclosed in parentheses.

For example, this is what GiNaC’s ‘abs()’ function is defined like:

static ex abs_eval(const ex & arg) { ... }
static ex abs_evalf(const ex & arg) { ... }
                                                                                
static void abs_print_latex(const ex & arg, const print_context & c)
{
    c.s << "{|"; arg.print(c); c.s << "|}";
}
                                                                                
static void abs_print_csrc_float(const ex & arg, const print_context & c)
{
    c.s << "fabs("; arg.print(c); c.s << ")";
}
                                                                                
REGISTER_FUNCTION(abs, eval_func(abs_eval).
                       evalf_func(abs_evalf).
                       print_func<print_latex>(abs_print_latex).
                       print_func<print_csrc_float>(abs_print_csrc_float).
                       print_func<print_csrc_double>(abs_print_csrc_float));

This will display ‘abs(x)’ as ‘|x|’ in LaTeX mode and fabs(x) in non-CLN C source output, but as abs(x) in all other formats.

There is currently no equivalent of set_print_func() for functions.

6.3.3 Adding new output formats

Creating a new output format involves subclassing print_context, which is somewhat similar to adding a new algebraic class (see Adding classes). There is a macro GINAC_DECLARE_PRINT_CONTEXT that needs to go into the class definition, and a corresponding macro GINAC_IMPLEMENT_PRINT_CONTEXT that has to appear at global scope. Every print_context class needs to provide a default constructor and a constructor from an std::ostream and an unsigned options value.

Here is an example for a user-defined print_context class:

class print_myformat : public print_dflt
{
    GINAC_DECLARE_PRINT_CONTEXT(print_myformat, print_dflt)
public:
    print_myformat(std::ostream & os, unsigned opt = 0)
     : print_dflt(os, opt) {}
};

print_myformat::print_myformat() : print_dflt(std::cout) {}

GINAC_IMPLEMENT_PRINT_CONTEXT(print_myformat, print_dflt)

That’s all there is to it. None of the actual expression output logic is implemented in this class. It merely serves as a selector for choosing a particular format. The algorithms for printing expressions in the new format are implemented as print methods, as described above.

print_myformat is a subclass of print_dflt, so it behaves exactly like GiNaC’s default output format:

{
    symbol x("x");
    ex e = pow(x, 2) + 1;

    // this prints "1+x^2"
    cout << e << endl;
    
    // this also prints "1+x^2"
    e.print(print_myformat()); cout << endl;

    ...
}

To fill print_myformat with life, we need to supply appropriate print methods with set_print_func(), like this:

// This prints powers with '**' instead of '^'. See the LaTeX output
// example above for explanations.
void print_power_as_myformat(const power & p,
                             const print_myformat & c,
                             unsigned level)
{
    unsigned power_prec = p.precedence();
    if (level >= power_prec)
        c.s << '(';
    p.op(0).print(c, power_prec);
    c.s << "**";
    p.op(1).print(c, power_prec);
    if (level >= power_prec)
        c.s << ')';
}

{
    ...
    // install a new print method for power objects
    set_print_func<power, print_myformat>(print_power_as_myformat);

    // now this prints "1+x**2"
    e.print(print_myformat()); cout << endl;

    // but the default format is still "1+x^2"
    cout << e << endl;
}

6.4.1 Example: scalar products

Let’s suppose that we need a way to handle some kind of abstract scalar product of the form ‘<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++ struct:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

struct sprod_s {
    ex left, right;

    sprod_s() {}
    sprod_s(ex l, ex r) : left(l), right(r) {}
};

The default constructor is required. Now, to make a GiNaC class out of this data structure, we need only one line:

typedef structure<sprod_s> sprod;

That’s it. This line constructs an algebraic class sprod which contains objects of type sprod_s. We can now use sprod in expressions like any other GiNaC class:

...
    symbol a("a"), b("b");
    ex e = sprod(sprod_s(a, b));
...

Note the difference between sprod which is the algebraic class, and sprod_s which is the unadorned C++ structure containing the left and right data members. As shown above, an sprod can be constructed from an sprod_s object.

If you find the nested sprod(sprod_s()) constructor too unwieldy, you could define a little wrapper function like this:

inline ex make_sprod(ex left, ex right)
{
    return sprod(sprod_s(left, right));
}

The sprod_s object contained in sprod can be accessed with the GiNaC ex_to<>() function followed by the -> operator or get_struct():

...
    cout << ex_to<sprod>(e)->left << endl;
     // -> a
    cout << ex_to<sprod>(e).get_struct().right << endl;
     // -> b
...

You only have read access to the members of sprod_s.

The type definition of sprod is enough to write your own algorithms that deal with scalar products, for 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>
...

6.4.2 Structure output

While the sprod type is useable it still leaves something to be desired, most notably proper output:

...
    cout << e << endl;
     // -> [structure object]
...

By default, any structure types you define will be printed as ‘[structure object]’. To override this you can either specialize the template’s print() member function, or specify print methods with set_print_func<>(), as described in GiNaC’s expression output system. Unfortunately, it’s not possible to supply class options like print_func<>() to structures, so for a self-contained structure type you need to resort to overriding the print() function, which is also what we will do here.

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:

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 << ">";
}

Now we can print expressions containing scalar products:

...
    cout << e << endl;
     // -> <a|b>
    cout << swap_sprod(e) << endl;
     // -> <b|a>
...

6.4.3 Comparing structures

The 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 sprod_s. This can lead to some confusing and undesired behavior:

...
    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)
...

To remedy this, we first need to define the operators == and < for objects of type sprod_s:

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;
}

The ordering established by the < operator doesn’t have to make any algebraic sense, but it needs to be well defined. Note that we can’t use expressions like lhs.left == rhs.left or lhs.left < rhs.left in the implementation of these operators because they would construct GiNaC relational objects which in the case of < 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 sprod type to let GiNaC know that an ordering relation exists for the embedded objects:

typedef structure<sprod_s, compare_std_less> sprod;

sprod objects then behave as expected:

...
    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>
...

The compare_std_less policy parameter tells GiNaC to use the std::less and std::equal_to functors to compare objects of type sprod_s. By default, these functors forward their work to the standard < and == operators, which we have overloaded. Alternatively, we could have specialized std::less and std::equal_to for class sprod_s.

GiNaC provides two other comparison policies for structure<T> objects: the default compare_all_equal, and compare_bitwise which does a bit-wise comparison of the contained 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 T class might have.

6.4.4 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 nops() and op() methods:

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");
    }
}

Implementing nops() and op() for container types such as sprod has two other nice side effects:

• has() works as expected
• GiNaC generates better hash keys for the objects (the default implementation of calchash() takes subexpressions into account)

There is a non-const variant of op() called let_op() that allows replacing subexpressions:

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");
    }
}

Once we have provided let_op() we also get subs() and map() for free. In fact, every container class that returns a non-null nops() value must either implement let_op() or provide custom implementations of subs() and map().

In turn, the availability of map() enables the recursive behavior of a couple of other default method implementations, in particular evalf(), evalm(), normal(), diff() and expand(). Although we probably want to provide our own version of expand() for scalar products that turns expressions like ‘<a+b|c>’ into ‘<a|c>+<b|c>’. This is left as an exercise for the reader.

The 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 eval()). There is, however, one topic that shall be addressed here, as it demonstrates one peculiarity of the structure<T> template: archiving.

6.4.5 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 archive() member function and the unarchiving constructor (the 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 class_name() member function) is ‘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 get_class_name() member function. Here is a sample implementation for enabling archiving of the scalar product type defined above:

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);
}

Note that the unarchiving constructor is sprod::structure and not sprod::sprod, and that we don’t need to supply an sprod::unarchive() function.

6.5.1 Hierarchy of algebraic classes.

All algebraic classes (that is, all classes that can appear in expressions) in GiNaC are direct or indirect subclasses of the class basic. So a basic * represents a generic pointer to an algebraic class. Working with such pointers directly is cumbersome (think of memory management), hence GiNaC wraps them into ex (see Expressions are reference counted). To make such wrapping possible every algebraic class has to implement several methods. Visitors (see Visitors and tree traversal), printing, and (un)archiving (see Input and output of expressions) require helper methods too. But don’t worry, most of the work is simplified by the following macros (defined in registrar.h):

• GINAC_DECLARE_REGISTERED_CLASS
• GINAC_IMPLEMENT_REGISTERED_CLASS
• GINAC_IMPLEMENT_REGISTERED_CLASS_OPT

The GINAC_DECLARE_REGISTERED_CLASS macro inserts declarations required for memory management, visitors, printing, and (un)archiving. It takes the name of the class and its direct superclass as arguments. The GINAC_DECLARE_REGISTERED_CLASS should be the first line after the opening brace of the class definition.

GINAC_IMPLEMENT_REGISTERED_CLASS takes the same arguments as GINAC_DECLARE_REGISTERED_CLASS. It initializes certain static members of a class so that printing and (un)archiving works. The GINAC_IMPLEMENT_REGISTERED_CLASS may appear anywhere else in the source (at global scope, of course, not inside a function).

GINAC_IMPLEMENT_REGISTERED_CLASS_OPT is a variant of GINAC_IMPLEMENT_REGISTERED_CLASS. It allows specifying additional options, such as custom printing functions.

6.5.2 A minimalistic example

Now we will start implementing a new class mystring that allows placing character strings in algebraic expressions (this is not very useful, but it’s just an example). This class will be a direct subclass of basic. You can use this sample implementation as a starting point for your own classes ³.

The code snippets given here assume that you have included some header files as follows:

#include <iostream>
#include <string>   
#include <stdexcept>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

Now we can write down the class declaration. The class stores a C++ string and the user shall be able to construct a mystring object from a string:

class mystring : public basic
{
    GINAC_DECLARE_REGISTERED_CLASS(mystring, basic)
  
public:
    mystring(const string & s);

private:
    string str;
};

GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)

The GINAC_DECLARE_REGISTERED_CLASS macro insert declarations required for memory management, visitors, printing, and (un)archiving. GINAC_IMPLEMENT_REGISTERED_CLASS initializes certain static members of a class so that printing and (un)archiving works.

Now there are three member functions we have to implement to get a working class:

• mystring(), the default constructor.
• 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 other object. If it returns 0, the objects are considered equal. Please notice: This has nothing to do with the (numeric) ordering relationship expressed by <, >= etc (which cannot be defined for non-numeric classes). For example, numeric(1).compare_same_type(numeric(2)) may return +1 even though 1 is clearly smaller than 2. Every GiNaC class must provide a compare_same_type() function, even those representing objects for which no reasonable algebraic ordering relationship can be defined.
• And, of course, mystring(const string& s) which is the constructor we declared.

Let’s proceed step-by-step. The default constructor looks like this:

mystring::mystring() { }

In the default constructor you should set all other member variables to reasonable default values (we don’t need that here since our str member gets set to an empty string automatically).

Our compare_same_type() function uses a provided function to compare the string members:

int mystring::compare_same_type(const basic & other) const
{
    const mystring &o = static_cast<const mystring &>(other);
    int cmpval = str.compare(o.str);
    if (cmpval == 0)
        return 0;
    else if (cmpval < 0)
        return -1;
    else
        return 1;
}

Although this function takes a basic &, it will always be a reference to an object of exactly the same class (objects of different classes are not comparable), so the cast is safe. If this function returns 0, the two objects are considered equal (in the sense that A-B=0), so you should compare all relevant member variables.

Now the only thing missing is our constructor:

mystring::mystring(const string& s) : str(s) { }

No surprises here. We set the str member from the argument.

That’s it! We now have a minimal working GiNaC class that can store strings in algebraic expressions. Let’s confirm that the RTTI works:

ex e = mystring("Hello, world!");
cout << is_a<mystring>(e) << endl;
 // -> 1 (true)

cout << ex_to<basic>(e).class_name() << endl;
 // -> mystring

Obviously it does. Let’s see what the expression e looks like:

cout << e << endl;
 // -> [mystring object]

Hm, not exactly what we expect, but of course the mystring class doesn’t yet know how to print itself. This can be done either by implementing the print() member function, or, preferably, by specifying a print_func<>() class option. Let’s say that we want to print the string surrounded by double quotes:

class mystring : public basic
{
    ...
protected:
    void do_print(const print_context & c, unsigned level = 0) const;
    ...
};

void mystring::do_print(const print_context & c, unsigned level) const
{
    // print_context::s is a reference to an ostream
    c.s << '\"' << str << '\"';
}

The level argument is only required for container classes to correctly parenthesize the output.

Now we need to tell GiNaC that mystring objects should use the do_print() member function for printing themselves. For this, we replace the line

GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)

with

GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mystring, basic,
  print_func<print_context>(&mystring::do_print))

Let’s try again to print the expression:

cout << e << endl;
 // -> "Hello, world!"

Much better. If we wanted to have mystring objects displayed in a different way depending on the output format (default, LaTeX, etc.), we would have supplied multiple print_func<>() options with different template parameters (print_dflt, print_latex, etc.), separated by dots. This is similar to the way options are specified for symbolic functions. See GiNaC’s expression output system, for a more in-depth description of the way expression output is implemented in GiNaC.

The mystring class can be used in arbitrary expressions:

e += mystring("GiNaC rulez"); 
cout << e << endl;
 // -> "GiNaC rulez"+"Hello, world!"

(GiNaC’s automatic term reordering is in effect here), or even

e = pow(mystring("One string"), 2*sin(Pi-mystring("Another string")));
cout << e << endl;
 // -> "One string"^(2*sin(-"Another string"+Pi))

Whether this makes sense is debatable but remember that this is only an example. At least it allows you to implement your own symbolic algorithms for your objects.

Note that GiNaC’s algebraic rules remain unchanged:

e = mystring("Wow") * mystring("Wow");
cout << e << endl;
 // -> "Wow"^2

e = pow(mystring("First")-mystring("Second"), 2);
cout << e.expand() << endl;
 // -> -2*"First"*"Second"+"First"^2+"Second"^2

There’s no way to, for example, make GiNaC’s add class perform string concatenation. You would have to implement this yourself.

6.5.3 Automatic evaluation

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. This is done in the evaluation member function eval(). Let’s say that we wanted all strings automatically converted to lowercase with non-alphabetic characters stripped, and empty strings removed:

class mystring : public basic
{
    ...
public:
    ex eval() const override;
    ...
};

ex mystring::eval() const
{
    string new_str;
    for (size_t i=0; i<str.length(); i++) {
        char c = str[i];
        if (c >= 'A' && c <= 'Z') 
            new_str += tolower(c);
        else if (c >= 'a' && c <= 'z')
            new_str += c;
    }

    if (new_str.length() == 0)
        return 0;

    return mystring(new_str).hold();
}

The 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 return this->hold();.

If our class had subobjects, we would have to evaluate them first (unless they are all of type ex, which are automatically evaluated). We don’t have any subexpressions in the mystring class, so we are not concerned with this.

Let’s confirm that it works:

ex e = mystring("Hello, world!") + mystring("!?#");
cout << e << endl;
 // -> "helloworld"

e = mystring("Wow!") + mystring("WOW") + mystring(" W ** o ** W");  
cout << e << endl;
 // -> 3*"wow"

6.5.4 Optional member functions

We have implemented only a small set of member functions to make the class 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:

unsigned calchash() const override;
bool is_equal_same_type(const basic & other) const override;

The calchash() method returns an 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 calchash() calculates a hash value out of the tinfo_key of the class and all subexpressions that are accessible via op().

is_equal_same_type() works like compare_same_type() but only tests for equality without establishing an ordering relation, which is often faster. The default implementation of is_equal_same_type() just calls compare_same_type() and tests its result for zero.

6.5.5 Other member functions

For a real algebraic class, there are probably some more functions that you might want to provide:

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;

If your class stores sub-expressions (see the scalar product example in the previous section) you will probably want to override

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;

let_op() is a variant of op() that allows write access. The default implementations of subs() and map() use it, so you have to implement either let_op(), or subs() and 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 ex by stripping the ex off using the ex_to<mystring>(e) function when that should become a need.

That’s it. May the source be with you!

6.5.6 Upgrading extension classes from older version of GiNaC

GiNaC used to use a custom run time type information system (RTTI). It was removed from GiNaC. Thus, one needs to rewrite constructors which set tinfo_key (which does not exist any more). For example,

myclass::myclass() : inherited(&myclass::tinfo_static) {}

needs to be rewritten as

myclass::myclass() {}