5. Methods and Functions

In this chapter the most important algorithms provided by GiNaC will be described. Some of them are implemented as functions on expressions, others are implemented as methods provided by expression objects. If they are methods, there exists a wrapper function around it, so you can alternatively call it in a functional way as shown in the simple example:

    ...
    cout << "As method:   " << sin(1).evalf() << endl;
    cout << "As function: " << evalf(sin(1)) << endl;
    ...

The general rule is that wherever methods accept one or more parameters (arg1, arg2, ...) the order of arguments the function wrapper accepts is the same but preceded by the object to act on (object, arg1, arg2, ...). This approach is the most natural one in an OO model but it may lead to confusion for MapleV users because where they would type A:=x+1; subs(x=2,A); GiNaC would require A=x+1; subs(A,x==2); (after proper declaration of A and x). On the other hand, since MapleV returns 3 on A:=x^2+3; coeff(A,x,0); (GiNaC: A=pow(x,2)+3; coeff(A,x,0);) it is clear that MapleV is not trying to be consistent here. Also, users of MuPAD will in most cases feel more comfortable with GiNaC's convention. All function wrappers are implemented as simple inline functions which just call the corresponding method and are only provided for users uncomfortable with OO who are dead set to avoid method invocations. Generally, nested function wrappers are much harder to read than a sequence of methods and should therefore be avoided if possible. On the other hand, not everything in GiNaC is a method on class ex and sometimes calling a function cannot be avoided.

5.1 Getting information about expressions

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;
unsigned 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 section 4.4 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 section 4.4 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 section 4.15 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), 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 section 4.11 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.

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 std::binary_function<ex, ex, bool> {
public:
    bool operator()(const ex &lh, const ex &rh) const;
};

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

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(),
                                  std::bind2nd(ex_is_equal(), 1));

The implementation of ex_is_less uses the member function

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

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

5.2 Numerical Evaluation

GiNaC keeps algebraic expressions, numbers and constants in their exact form. To evaluate them using floating-point arithmetic you need to call

ex ex::evalf(int level = 0) const;

The accuracy of the evaluation is controlled by the global object Digits which can be assigned an integer value. The default value of Digits is 17. See section 4.6 Numbers, for more information and examples.

To evaluate an expression to a double floating-point number you can call evalf() followed by numeric::to_double(), like this:

{
    // Approximate sin(x/Pi)
    symbol x("x");
    ex e = series(sin(x/Pi), x == 0, 6);

    // Evaluate numerically at x=0.1
    ex f = evalf(e.subs(x == 0.1));

    // ex_to<numeric> is an unsafe cast, so check the type first
    if (is_a<numeric>(f)) {
        double d = ex_to<numeric>(f).to_double();
        cout << d << endl;
         // -> 0.0318256
    } else
        // error
}

5.3 Substituting expressions

Algebraic objects inside expressions can be replaced with arbitrary expressions via the .subs() method:

ex ex::subs(const ex & e, unsigned options = 0);
ex ex::subs(const exmap & m, unsigned options = 0);
ex ex::subs(const lst & syms, const lst & repls, unsigned options = 0);

In the first form, subs() accepts a relational of the form `object == expression' or a lst of such relationals:

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

    ex e1 = 2*x^2-4*x+3;
    cout << "e1(7) = " << e1.subs(x == 7) << endl;
     // -> 73

    ex e2 = x*y + x;
    cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
     // -> -10
}

If you specify multiple substitutions, they are performed in parallel, so e.g. subs(lst(x == y, y == x)) exchanges `x' and `y'.

The second form of subs() takes an exmap object which is a pair associative container that maps expressions to expressions (currently implemented as a std::map). This is the most efficient one of the three subs() forms and should be used when the number of objects to be substituted is large or unknown.

Using this form, the second example from above would look like this:

{
    symbol x("x"), y("y");
    ex e2 = x*y + x;

    exmap m;
    m[x] = -2;
    m[y] = 4;
    cout << "e2(-2, 4) = " << e2.subs(m) << endl;
}

The third form of subs() takes two lists, one for the objects to be replaced and one for the expressions to be substituted (both lists must contain the same number of elements). Using this form, you would write

{
    symbol x("x"), y("y");
    ex e2 = x*y + x;

    cout << "e2(-2, 4) = " << e2.subs(lst(x, y), lst(-2, 4)) << endl;
}

The optional last argument to subs() is a combination of subs_options flags. There are two options available: subs_options::no_pattern disables pattern matching, which makes large subs() operations significantly faster if you are not using patterns. The second option, subs_options::algebraic enables algebraic substitutions in products and powers. 5.4 Pattern matching and advanced substitutions, for more information about patterns and algebraic substitutions.

subs() performs syntactic substitution of any complete algebraic object; it does not try to match sub-expressions as is demonstrated by the following example:

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

    ex e1 = pow(x+y, 2);
    cout << e1.subs(x+y == 4) << endl;
     // -> 16

    ex e2 = sin(x)*sin(y)*cos(x);
    cout << e2.subs(sin(x) == cos(x)) << endl;
     // -> cos(x)^2*sin(y)

    ex e3 = x+y+z;
    cout << e3.subs(x+y == 4) << endl;
     // -> x+y+z
     // (and not 4+z as one might expect)
}

A more powerful form of substitution using wildcards is described in the next section.

5.4 Pattern matching and advanced substitutions

GiNaC allows the use of patterns for checking whether an expression is of a certain form or contains subexpressions of a certain form, and for substituting expressions in a more general way.

A pattern is an algebraic expression that optionally contains wildcards. A wildcard is a special kind of object (of class wildcard) that represents an arbitrary expression. Every wildcard has a label which is an unsigned integer number to allow having multiple different wildcards in a pattern. Wildcards are printed as `$label' (this is also the way they are specified in ginsh). In C++ code, wildcard objects are created with the call

ex wild(unsigned label = 0);

which is simply a wrapper for the wildcard() constructor with a shorter name.

Some examples for patterns:

Notes:

• Wildcards behave like symbols and are subject to the same algebraic rules. E.g., `$0+2*$0' is automatically transformed to `3*$0'.
• As shown in the last example, to use wildcards for indices you have to use them as the value of an idx object. This is because indices must always be of class idx (or a subclass).
• Wildcards only represent expressions or subexpressions. It is not possible to use them as placeholders for other properties like index dimension or variance, representation labels, symmetry of indexed objects etc.
• Because wildcards are commutative, it is not possible to use wildcards as part of noncommutative products.
• A pattern does not have to contain wildcards. `x' and `x+y' are also valid patterns.

5.4.1 Matching expressions

bool ex::match(const ex & pattern);
bool ex::match(const ex & pattern, lst & 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 repls object as a list of relations of the form `wildcard == expression'. If match() returns false, the state of repls is undefined. For reproducible results, the list should be empty when passed to match(), but it is also possible to find similarities in multiple expressions by passing in the result of a previous match.

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==c,$2==b}
  (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

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, lst & 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

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 Algebraic substitutions

An example clarifies it all (hopefully):

cout << (a*a*a*a+b*b*b*b+pow(x+y,4)).subs(wild()*wild()==pow(wild(),3),
                                        subs_options::algebraic) << endl;
// --> (y+x)^6+b^6+a^6

cout << ((a+b+c)*(a+b+c)).subs(a+b==x,subs_options::algebraic) << endl;
// --> (c+b+a)^2
// Powers and products are smart, but addition is just the same.

cout << ((a+b+c)*(a+b+c)).subs(a+b+wild()==x+wild(), subs_options::algebraic)
                                                                      << endl;
// --> (x+c)^2
// As I said: addition is just the same.

cout << (pow(a,5)*pow(b,7)+2*b).subs(b*b*a==x,subs_options::algebraic) << endl;
// --> x^3*b*a^2+2*b

cout << (pow(a,-5)*pow(b,-7)+2*b).subs(1/(b*b*a)==x,subs_options::algebraic)
                                                                       << endl;
// --> 2*b+x^3*b^(-1)*a^(-2)

cout << (4*x*x*x-2*x*x+5*x-1).subs(x==a,subs_options::algebraic) << endl;
// --> -1-2*a^2+4*a^3+5*a

cout << (4*x*x*x-2*x*x+5*x-1).subs(pow(x,wild())==pow(a,wild()),
                                subs_options::algebraic) << endl;
// --> -1+5*x+4*x^3-2*x^2
// You should not really need this kind of patterns very often now.
// But perhaps this it's-not-a-bug-it's-a-feature (c/sh)ould still change.

cout << ex(sin(1+sin(x))).subs(sin(wild())==cos(wild()),
                                subs_options::algebraic) << endl;
// --> cos(1+cos(x))

cout << expand((a*sin(x+y)*sin(x+y)+a*cos(x+y)*cos(x+y)+b)
        .subs((pow(cos(wild()),2)==1-pow(sin(wild()),2)),
                                subs_options::algebraic)) << endl;
// --> b+a

5.5 Applying a Function on Subexpressions

Sometimes you may want to perform an operation on specific parts of an expression while leaving the general structure of it intact. An example of this would be a matrix trace operation: the trace of a sum is the sum of the traces of the individual terms. That is, the trace should map on the sum, by applying itself to each of the sum's operands. It is possible to do this manually which usually results in code like this:

ex calc_trace(ex e)
{
    if (is_a<matrix>(e))
        return ex_to<matrix>(e).trace();
    else if (is_a<add>(e)) {
        ex sum = 0;
        for (size_t i=0; i<e.nops(); i++)
            sum += calc_trace(e.op(i));
        return sum;
    } else if (is_a<mul>)(e)) {
        ...
    } else {
        ...
    }
}

This is, however, slightly inefficient (if the sum is very large it can take a long time to add the terms one-by-one), and its applicability is limited to a rather small class of expressions. If calc_trace() is called with a relation or a list as its argument, you will probably want the trace to be taken on both sides of the relation or of all elements of the list.

GiNaC offers the map() method to aid in the implementation of such operations:

ex ex::map(map_function & f) const;
ex ex::map(ex (*f)(const ex & e)) const;

In the first (preferred) form, map() takes a function object that is subclassed from the map_function class. In the second form, it takes a pointer to a function that accepts and returns an expression. map() constructs a new expression of the same type, applying the specified function on all subexpressions (in the sense of op()), non-recursively.

The use of a function object makes it possible to supply more arguments to the function that is being mapped, or to keep local state information. The map_function class declares a virtual function call operator that you can overload. Here is a sample implementation of calc_trace() that uses map() in a recursive fashion:

struct calc_trace : public map_function {
    ex operator()(const ex &e)
    {
        if (is_a<matrix>(e))
            return ex_to<matrix>(e).trace();
        else if (is_a<mul>(e)) {
            ...
        } else
            return e.map(*this);
    }
};

This function object could then be used like this:

{
    ex M = ... // expression with matrices
    calc_trace do_trace;
    ex tr = do_trace(M);
}

Here is another example for you to meditate over. It removes quadratic terms in a variable from an expanded polynomial:

struct map_rem_quad : public map_function {
    ex var;
    map_rem_quad(const ex & var_) : var(var_) {}

    ex operator()(const ex & e)
    {
        if (is_a<add>(e) || is_a<mul>(e))
     	    return e.map(*this);
        else if (is_a<power>(e) && 
                 e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
            return 0;
        else
            return e;
    }
};

...

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

    ex e;
    for (int i=0; i<8; i++)
        e += pow(x, i) * pow(y, 8-i) * (i+1);
    cout << e << endl;
     // -> 4*y^5*x^3+5*y^4*x^4+8*y*x^7+7*y^2*x^6+2*y^7*x+6*y^3*x^5+3*y^6*x^2+y^8

    map_rem_quad rem_quad(x);
    cout << rem_quad(e) << endl;
     // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
}

ginsh offers a slightly different implementation of map() that allows applying algebraic functions to operands. The second argument to map() is an expression containing the wildcard `$0' which acts as the placeholder for the operands:

> map(a*b,sin($0));
sin(a)*sin(b)
> map(a+2*b,sin($0));
sin(a)+sin(2*b)
> map({a,b,c},$0^2+$0);
{a^2+a,b^2+b,c^2+c}

Note that it is only possible to use algebraic functions in the second argument. You can not use functions like `diff()', `op()', `subs()' etc. because these are evaluated immediately:

> map({a,b,c},diff($0,a));
{0,0,0}
  This is because "diff($0,a)" evaluates to "0", so the command is equivalent
  to "map({a,b,c},0)".

5.6 Visitors and Tree Traversal

Suppose that you need a function that returns a list of all indices appearing in an arbitrary expression. The indices can have any dimension, and for indices with variance you always want the covariant version returned.

You can't use get_free_indices() because you also want to include dummy indices in the list, and you can't use find() as it needs specific index dimensions (and it would require two passes: one for indices with variance, one for plain ones).

The obvious solution to this problem is a tree traversal with a type switch, such as the following:

void gather_indices_helper(const ex & e, lst & l)
{
    if (is_a<varidx>(e)) {
        const varidx & vi = ex_to<varidx>(e);
        l.append(vi.is_covariant() ? vi : vi.toggle_variance());
    } else if (is_a<idx>(e)) {
        l.append(e);
    } else {
        size_t n = e.nops();
        for (size_t i = 0; i < n; ++i)
            gather_indices_helper(e.op(i), l);
    }
}

lst gather_indices(const ex & e)
{
    lst l;
    gather_indices_helper(e, l);
    l.sort();
    l.unique();
    return l;
}

This works fine but fans of object-oriented programming will feel uncomfortable with the type switch. One reason is that there is a possibility for subtle bugs regarding derived classes. If we had, for example, written

    if (is_a<idx>(e)) {
      ...
    } else if (is_a<varidx>(e)) {
      ...

in gather_indices_helper, the code wouldn't have worked because the first line "absorbs" all classes derived from idx, including varidx, so the special case for varidx would never have been executed.

Also, for a large number of classes, a type switch like the above can get unwieldy and inefficient (it's a linear search, after all). gather_indices_helper only checks for two classes, but if you had to write a function that required a different implementation for nearly every GiNaC class, the result would be very hard to maintain and extend.

The cleanest approach to the problem would be to add a new virtual function to GiNaC's class hierarchy. In our example, there would be specializations for idx and varidx while the default implementation in basic performed the tree traversal. Unfortunately, in C++ it's impossible to add virtual member functions to existing classes without changing their source and recompiling everything. GiNaC comes with source, so you could actually do this, but for a small algorithm like the one presented this would be impractical.

One solution to this dilemma is the Visitor design pattern, which is implemented in GiNaC (actually, Robert Martin's Acyclic Visitor variation, described in detail in http://objectmentor.com/publications/acv.pdf). Instead of adding virtual functions to the class hierarchy to implement operations, GiNaC provides a single "bouncing" method accept() that takes an instance of a special visitor class and redirects execution to the one visit() virtual function of the visitor that matches the type of object that accept() was being invoked on.

Visitors in GiNaC must derive from the global visitor class as well as from the class T::visitor of each class T they want to visit, and implement the member functions void visit(const T &) for each class.

A call of

void ex::accept(visitor & v) const;

will then dispatch to the correct visit() member function of the specified visitor v for the type of GiNaC object at the root of the expression tree (e.g. a symbol, an idx or a mul).

Here is an example of a visitor:

class my_visitor
 : public visitor,          // this is required
   public add::visitor,     // visit add objects
   public numeric::visitor, // visit numeric objects
   public basic::visitor    // visit basic objects
{
    void visit(const add & x)
    { cout << "called with an add object" << endl; }

    void visit(const numeric & x)
    { cout << "called with a numeric object" << endl; }

    void visit(const basic & x)
    { cout << "called with a basic object" << endl; }
};

which can be used as follows:

...
    symbol x("x");
    ex e1 = 42;
    ex e2 = 4*x-3;
    ex e3 = 8*x;

    my_visitor v;
    e1.accept(v);
     // prints "called with a numeric object"
    e2.accept(v);
     // prints "called with an add object"
    e3.accept(v);
     // prints "called with a basic object"
...

The visit(const basic &) method gets called for all objects that are not numeric or add and acts as an (optional) default.

From a conceptual point of view, the visit() methods of the visitor behave like a newly added virtual function of the visited hierarchy. In addition, visitors can store state in member variables, and they can be extended by deriving a new visitor from an existing one, thus building hierarchies of visitors.

We can now rewrite our index example from above with a visitor:

class gather_indices_visitor
 : public visitor, public idx::visitor, public varidx::visitor
{
    lst l;

    void visit(const idx & i)
    {
        l.append(i);
    }

    void visit(const varidx & vi)
    {
        l.append(vi.is_covariant() ? vi : vi.toggle_variance());
    }

public:
    const lst & get_result() // utility function
    {
        l.sort();
        l.unique();
        return l;
    }
};

What's missing is the tree traversal. We could implement it in visit(const basic &), but GiNaC has predefined methods for this:

void ex::traverse_preorder(visitor & v) const;
void ex::traverse_postorder(visitor & v) const;
void ex::traverse(visitor & v) const;

traverse_preorder() visits a node before visiting its subexpressions, while traverse_postorder() visits a node after visiting its subexpressions. traverse() is a synonym for traverse_preorder().

Here is a new implementation of gather_indices() that uses the visitor and traverse():

lst gather_indices(const ex & e)
{
    gather_indices_visitor v;
    e.traverse(v);
    return v.get_result();
}

Alternatively, you could use pre- or postorder iterators for the tree traversal:

lst gather_indices(const ex & e)
{
    gather_indices_visitor v;
    for (const_preorder_iterator i = e.preorder_begin();
         i != e.preorder_end(); ++i) {
        i->accept(v);
    }
    return v.get_result();
}

5.7 Polynomial arithmetic

5.7.1 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 (written down here in output-style). It is equivalent to the factorized polynomial . 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 and .

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 . 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 , 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.2 Degree and coefficients

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

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

which also work reliably on non-expanded input polynomials (they even work on rational functions, returning the asymptotic degree). By definition, the degree of zero is zero. To extract a coefficient with a certain power from an expanded polynomial you use

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

The additional function

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

computes the pseudo-remainder of `a' and `b' which satisfies , where .

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.4 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.5 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.6 Square-free decomposition

GiNaC still lacks proper factorization support. Some form of factorization is, however, easily implemented by noting that factors appearing in a polynomial with power two or more also appear in the derivative and hence can easily be found by computing the GCD of the original polynomial and its derivatives. Any decent system has an interface for this so called square-free factorization. So we provide one, too:

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

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

5.8 Rational expressions

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, calling numer_denom() is faster than using numer() and denom() separately.

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);
ex ex::to_polynomial(lst & l);

ex ex::to_rational(exmap & m);
ex ex::to_rational(lst & l);

on the expression to be converted. The supplied exmap or lst 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;

    lst lp;
    ex p = a.to_polynomial(lp);
    cout << " = " << p << "\n   with " << lp << endl;
     // = symbol3*symbol2*y+2*symbol2*x
     //   with {symbol2==sin(x)^(-1),symbol3==-1/3}

    lst lr;
    ex r = a.to_rational(lr);
    cout << " = " << r << "\n   with " << lr << endl;
     // = -1/3*symbol4^(-1)*y+2*symbol4^(-1)*x
     //   with {symbol4==sin(x)}
}

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.9 Symbolic differentiation

GiNaC's objects know how to differentiate themselves. Thus, a polynomial (class add) knows that its derivative is the sum of the derivatives of all the monomials:

{
    symbol x("x"), y("y"), z("z");
    ex P = pow(x, 5) + pow(x, 2) + y;

    cout << P.diff(x,2) << endl;
     // -> 20*x^3 + 2
    cout << P.diff(y) << endl;    // 1
     // -> 1
    cout << P.diff(z) << endl;    // 0
     // -> 0
}

If a second integer parameter n is given, the diff method returns the nth derivative.

If every object and every function is told what its derivative is, all derivatives of composed objects can be calculated using the chain rule and the product rule. Consider, for instance the expression 1/cosh(x). Since the derivative of cosh(x) is sinh(x) and the derivative of pow(x,-1) is -pow(x,-2), GiNaC can readily compute the composition. It turns out that the composition is the generating function for Euler Numbers, i.e. the so called nth Euler number is the coefficient of x^n/n! in the expansion of 1/cosh(x). We may use this identity to code a function that generates Euler numbers in just three lines:

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

ex EulerNumber(unsigned n)
{
    symbol x;
    const ex generator = pow(cosh(x),-1);
    return generator.diff(x,n).subs(x==0);
}

int main()
{
    for (unsigned i=0; i<11; i+=2)
        std::cout << EulerNumber(i) << std::endl;
    return 0;
}

When you run it, it produces the sequence 1, -1, 5, -61, 1385, -50521. We increment the loop variable i by two since all odd Euler numbers vanish anyways.

5.10 Series expansion

Expressions know how to expand themselves as a Taylor series or (more generally) a Laurent series. As in most conventional Computer Algebra Systems, no distinction is made between those two. There is a class of its own for storing such series (class pseries) and a built-in function (called Order) for storing the order term of the series. As a consequence, if you want to work with series, i.e. multiply two series, you need to call the method ex::series again to convert it to a series object with the usual structure (expansion plus order term). A sample application from special relativity could read:

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

int main()
{
    symbol v("v"), c("c");
    
    ex gamma = 1/sqrt(1 - pow(v/c,2));
    ex mass_nonrel = gamma.series(v==0, 10);
    
    cout << "the relativistic mass increase with v is " << endl
         << mass_nonrel << endl;
    
    cout << "the inverse square of this series is " << endl
         << pow(mass_nonrel,-2).series(v==0, 10) << endl;
}

Only calling the series method makes the last output simplify to , without that call we would just have a long series raised to the power .

As another instructive application, let us calculate the numerical value of Archimedes' constant (for which there already exists the built-in constant Pi) using John Machin's amazing formula . This equation (and similar ones) were used for over 200 years for computing digits of pi (see Pi Unleashed). We may expand the arcus tangent around 0 and insert the fractions 1/5 and 1/239. However, as we have seen, a series in GiNaC carries an order term with it and the question arises what the system is supposed to do when the fractions are plugged into that order term. The solution is to use the function series_to_poly() to simply strip the order term off:

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

ex machin_pi(int degr)
{
    symbol x;
    ex pi_expansion = series_to_poly(atan(x).series(x,degr));
    ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
                   -4*pi_expansion.subs(x==numeric(1,239));
    return pi_approx;
}

int main()
{
    using std::cout;  // just for fun, another way of...
    using std::endl;  // ...dealing with this namespace std.
    ex pi_frac;
    for (int i=2; i<12; i+=2) {
        pi_frac = machin_pi(i);
        cout << i << ":\t" << pi_frac << endl
             << "\t" << pi_frac.evalf() << endl;
    }
    return 0;
}

Note how we just called .series(x,degr) instead of .series(x==0,degr). This is a simple shortcut for ex's method series(): if the first argument is a symbol the expression is expanded in that symbol around point 0. When you run this program, it will type out:

2:      3804/1195
        3.1832635983263598326
4:      5359397032/1706489875
        3.1405970293260603143
6:      38279241713339684/12184551018734375
        3.141621029325034425
8:      76528487109180192540976/24359780855939418203125
        3.141591772182177295
10:     327853873402258685803048818236/104359128170408663038552734375
        3.1415926824043995174

5.11 Symmetrization

The three methods

ex ex::symmetrize(const lst & l);
ex ex::antisymmetrize(const lst & l);
ex ex::symmetrize_cyclic(const lst & l);

symmetrize an expression by returning the sum over all symmetric, antisymmetric or cyclic permutations of the specified list of objects, weighted by the number of permutations.

The three additional methods

ex ex::symmetrize();
ex ex::antisymmetrize();
ex ex::symmetrize_cyclic();

symmetrize or antisymmetrize an expression over its free indices.

Symmetrization is most useful with indexed expressions but can be used with almost any kind of object (anything that is subs()able):

{
    idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
    symbol A("A"), B("B"), a("a"), b("b"), c("c");
                                           
    cout << indexed(A, i, j).symmetrize() << endl;
     // -> 1/2*A.j.i+1/2*A.i.j
    cout << indexed(A, i, j, k).antisymmetrize(lst(i, j)) << endl;
     // -> -1/2*A.j.i.k+1/2*A.i.j.k
    cout << lst(a, b, c).symmetrize_cyclic(lst(a, b, c)) << endl;
     // -> 1/3*{a,b,c}+1/3*{b,c,a}+1/3*{c,a,b}
}

5.12 Predefined mathematical functions

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 for C++ as defined in the ANSI standard as far as possible. In particular: the natural logarithm (log) and the square root (sqrt) both have their branch cuts running along the negative real axis where the points on the axis itself belong to the upper part (i.e. continuous with quadrant II). The inverse trigonometric and hyperbolic functions are not defined for complex arguments by the C++ standard, however. In GiNaC we follow the conventions used by CLN, which in turn follow the carefully designed definitions in the Common Lisp standard. It should be noted that this convention is identical to the one used by the C99 standard and by most serious CAS. It is to be expected that future revisions of the C++ standard incorporate these functions in the complex domain in a manner compatible with C99.

5.12.2 Multiple polylogarithms

The multiple polylogarithm is the most generic member of a family of functions, to which others like the harmonic polylogarithm, Nielsen's generalized polylogarithm and the multiple zeta value belong. Everyone of these functions can also be written as a multiple polylogarithm with specific 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 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. 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 So in comparison to 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 and G(lst(a,b), lst(-1,1), c) means 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. 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.13 Complex Conjugation

The method

ex ex::conjugate();

returns the complex conjugate of the expression. For all built-in functions and objects the conjugation gives the expected results:

{
    varidx a(symbol("a"), 4), b(symbol("b"), 4);
    symbol x("x");
    realsymbol y("y");
                                           
    cout << (3*I*x*y + sin(2*Pi*I*y)).conjugate() << endl;
     // -> -3*I*conjugate(x)*y+sin(-2*I*Pi*y)
    cout << (dirac_gamma(a)*dirac_gamma(b)*dirac_gamma5()).conjugate() << endl;
     // -> -gamma5*gamma~b*gamma~a
}

For symbols in the complex domain the conjugation can not be evaluated and the GiNaC function conjugate is returned. GiNaC functions conjugate by applying the conjugation to their arguments. This is the default strategy. If you want to define your own functions and want to change this behavior, you have to supply a specialized conjugation method for your function (see 6.2 Symbolic functions and the GiNaC source-code for abs as an example).

5.14 Solving Linear Systems of Equations

The function lsolve() provides a convenient wrapper around some matrix operations that comes in handy when a system of linear equations needs to be solved:

ex lsolve(const ex & eqns, const ex & symbols,
          unsigned options = solve_algo::automatic);

Here, eqns is a lst of equalities (i.e. class relational) while symbols is a lst of indeterminates. (See section 4.4 The Class Hierarchy, for an exposition of class lst).

It returns the lst of solutions as an expression. As an example, let us solve the two equations a*x+b*y==3 and x-y==b:

{
    symbol a("a"), b("b"), x("x"), y("y");
    lst eqns, vars;
    eqns = a*x+b*y==3, x-y==b;
    vars = x, y;
    cout << lsolve(eqns, vars) << endl;
     // -> {x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)}

When the linear equations eqns are underdetermined, the solution will contain one or more tautological entries like x==x, depending on the rank of the system. When they are overdetermined, the solution will be an empty lst. Note the third optional parameter to lsolve(): it accepts the same parameters as matrix::solve(). This is because lsolve is just a wrapper around that method.

5.15 Input and output of expressions

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 construct an expression from a string, specifying the list of symbols to be used:

{
    symbol x("x"), y("y");
    ex e("2*x+sin(y)", lst(x, y));
}

The input syntax is the same as that used by ginsh and the stream output operator <<. The symbols in the string are matched by name to the symbols in the list and if GiNaC encounters a symbol not specified in the list it will throw an exception.

With this constructor, it's also easy to implement interactive GiNaC programs:

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

int main()
{
    symbol x("x");
    string s;

    cout << "Enter an expression containing 'x': ";
    getline(cin, s);

    try {
        ex e(s, lst(x));
        cout << "The derivative of " << e << " with respect to x is ";
        cout << e.diff(x) << ".\n";
    } catch (exception &p) {
        cerr << p.what() << endl;
    }
}

5.15.3 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>
using namespace std;
#include <ginac/ginac.h>
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");
    out << a;
    out.close();
    // ...

The file `foobar.gar' contains all information that is needed to reconstruct the expressions foo and bar.

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");
    in >> a2;
    // ...

And the stored expressions can be retrieved by their name:

    // ...
    lst syms;
    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.