# Algebraic Computation without Mathematica

The only real purpose of programming algebraic computations in a traditional language is efficiency. As noted previously, tensor operations involving large numbers of contractions do not scale well in a truly general-purpose environment such as that provided by Mathematica. We begin then by narrowing the scope of what we intend to accomplish to something that can be represented and manipulated easily and quickly.

As an example, consider the Kerr-Newman Metric. It can be parameterized as

gr r = -f / (2 h r)

gθ θ = -f / q

gφ φ = -a4 q / f + 2 a2 h q2 r / f -

2 a2 q r2 / f - q r4 / f
gφ t = a e2 q / f - 2 a M q r / f

gt t = -a2 c2 / f - e2 / f + 2 M r / f - r2 / f

with
f = a2 c2 + r2

h = (a2 + e2 + r (-2 M + r)) / (-2 r)

q = c2 - 1

By reducing the class of problems of interest to those that can be represented as polynomials with real coefficients and integer exponents, we

• simplify the storage and multiplication of terms, and parameterization functions and their derivatives can be pre-computed for easy chain-rule differentiation;
• make most computations integer, which are far faster than floating point;
• reduce memory requirements by the sparse storage of terms on linked lists (the lists are doubly-linked to facilitate insertion and removal).
This approach has several potential problems:

1. The parameterization of fractional exponents can be awkward, but is tractable:
1 / ( 1 + k / r3 ) ( 1 / n ) = 1 / q ->
qn = 1 + k / r3
so that, for instance,
n q ( n - 1 ) dq / dr = - 3 k / r4

and

dq / dr = - 3 k / ( n q( n - 1 ) r4 )

2. Constants must be handled as such: a, M, e, I need to have zero derivatives.
3. Recurring fractions are most efficiently handled through the definition of additional constants; this also reduces floating point representational error.
4. Ansatze are also awkward: if u(r) is not known at computation time, neither are its derivatives, so we parametrize them as well:
u'(r) = v(r) and

u''(r) = w(r)

Only two derivatives are needed as long as any computations only involve the Christoffel Symbols and the various curvature tensors.

### Implementation in C or C++

What follows is a series of suggestions, most of which apply to any programming project.

Sanity checks are an important part of any algorithm, especially one as general as we are attempting:

• Class members should be protected as much as possible.
• All unused or freed pointers should be set to NULL.
• Make sure classes or structures are unlinked from lists before freeing or deleting them.
• Abort on any unsuccessful malloc or new.
• Use magic numbers for each class or structure, and check them at each reference. Invert them when the structure is freed (A magic number is a unique constant identifier for each structure type).
• All indices should be checked for bounds violations before using them.
• Inter-process message tracing is necessary for debugging sequencing problems.
• Use
#ifdef debug
to allow some sanity checking to be easily turned off when the application is "stable".
A number of optimizations are also possible:

• Inlines can be used to speed procedure calls, although they can make debuggers crazy.
• Check for factors of 1 or 0 to avoid unnecessary computation.
• Free all unused memory as soon as possible.
• Linked lists are good for processing nonzero components, but lookups are expensive even with good hashing. For that reason, maintain a non-sparse array for each tensor with pointers to components to enable fast location of a component with specified indices.
• When raising or lowering flat-space lorentz indices, use the fact that the operation only reverses the sign of timelike indices and leaves the others unchanged.
• Precompute (once) and keep ε matrices in memory (make sure that their sum is zero to check consistency). ε components can be stored as
(index value) << (4 * index number) + ((value < 0) ? 1 << 50 : 0) in long longs.
• If you are only contracting ε symbols with antisymmetric indices in pairs, you will need 2D / 2 fewer components (the consistency sum should now equal (D/2)!).
• Gamma matrices (and their antisymmetric products) can likewise be precomputed, but you need to differentiate tensor and spinor indices, and allow for values of +- I:
((spinor index) ? (index value) << (4 * number of spacetime indices + 8 * spinor index number) :
(index value) << (4 * spacetime index number)) +

((value < 0) ? 1 << 50 : 0) + ((value = +- I) ? 1 << 51 : 0)

Software metrics are also a very good idea for any algorithm of this complexity:

• Record the number of entries into procedures and execution timings for long operations.

• Keep track of the maximum depths of hash overflows or searched linked lists.

### Classes used in PTAH

A good set of data structures will make any program easier to code and debug. It is for this reason that C++ was chosen as the implementation language for
PTAH (whose predecessor was written in C). The following classes are used in PTAH:

These classes are connected as follows:

### PVM - Parallel Virtual Machine

PVM is the underlying MIMD cluster software chosen for PTAH. It was chosen primarily because it provides tools to build in fault tolerance. PVM provides procedures for interrogating the cluster configuration, spawning tasks, sending and receiving messages (both data and control) between tasks and providing notification when a task dies. PVM does not guarantee message ordering, so that must be done by hand.

### Hybridization with Mathematica

PTAH provides Mathematica functions to write C++ header files for the exportation of tensors and parameters. C++ procedures in PTAH then build sums and ntuples from the information in the header files. PTAH output is formatted as Mathematica FullForm, providing easy importation into Mathematica for further analysis. PTAH also outputs the definition of a simplification function to replace parameterization functions by their definitions.

Here is a partial comparison of PTAH timings with Mathematica for the above analyses of KerrD (times are in seconds):

DMathematicaPTAHnumber of products computed by PTAH maximum PTAH working set (MB)
4356873287
5891711971511
67852118578718
715245223914428
85086*10633656644
91434830942773567
1039482*1332382770099
11953673117723598142

* = does not include computation of the Euler Class due to the excessive time required

These computations were performed on one 1533 MHz Athlon with 512 MB of DDR-266 RAM.