GiNaC allows you to handle expressions containing general indexed objects in arbitrary spaces. It is also able to canonicalize and simplify such expressions and perform symbolic dummy index summations. There are a number of predefined indexed objects provided, like delta and metric tensors.
There are few restrictions placed on indexed objects and their indices and it is easy to construct nonsense expressions, but our intention is to provide a general framework that allows you to implement algorithms with indexed quantities, getting in the way as little as possible.
Indexed expressions in GiNaC are constructed of two special types of objects, index objects and indexed objects.
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 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.
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 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.
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
...
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
...
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.
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:
get_free_indices() does
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’.
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).
The delta tensor takes two indices, is symmetric and has the matrix
representation diag(1, 1, 1, ...). It is constructed by the function
delta_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
}
The function metric_tensor() creates a general symmetric metric
tensor with two indices that can be used to raise/lower tensor indices. The
metric tensor is denoted as ‘g’ in the output and if its indices are of
mixed variance it is automatically replaced by a delta tensor:
{
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
}
The Minkowski metric tensor is a special metric tensor with a constant
matrix representation which is either diag(1, -1, -1, ...) (negative
signature, the default) or diag(-1, 1, 1, ...) (positive signature).
It is created with the function lorentz_g() (although it is output as
‘eta’):
{
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
}
The function spinor_metric() creates an antisymmetric tensor with
two indices that is used to raise/lower indices of 2-component spinors.
It is output as ‘eps’:
{
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
}
The matrix representation of the spinor metric is [[0, 1], [-1, 0]].
The epsilon tensor is totally antisymmetric, its number of indices is equal to the dimension of the index space (the indices must all be of the same numeric dimension), and ‘eps.1.2.3...’ (resp. ‘eps~0~1~2...’) is defined to be 1. Its behavior with indices that have a variance also depends on the signature of the metric. Epsilon tensors are output as ‘eps’.
There are three functions defined to create epsilon tensors in 2, 3 and 4 dimensions:
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);
The first two functions create an epsilon tensor in 2 or 3 Euclidean
dimensions, the last function creates an epsilon tensor in a 4-dimensional
Minkowski space (the last bool argument specifies whether the metric
has negative or positive signature, as in the case of the Minkowski metric
tensor):
{
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
}
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(2, 2), X(2, 1);
A = 1, 2,
3, 4;
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.