Generic Arithmetic

Writing generic algorithms that involve arithmetic is fairly straightforward. Unfortunately, the convenience of using generic arithmetic is limited by lack of support in programming languages.

Generic Operations

The Operations<T> type contains fields and methods that perform arithmetic operations and compute elementary functions for arbitrary types. The types that are supported out-of-the-box are: Double, Single, Quad, Decimal, SByte, Int16, Int32, Int64, Byte, UInt16, UInt32, UInt64, Boolean, BigInteger, BigRational, BigFloat. The details of how to implement these operations for your own types are discussed later.

Using the Operations<T> class is very straightforward. The following operations are available:

  • Constants: Zero, One, MinusOne, Epsilon, PositiveInfinity, NegativeInfinity, NaN.

  • Arithmetic operations: Add, Subtract, Negate, Abs, Multiply, Divide, Pow, Reciprocal.

  • Relational operations: Compare, LessThan, LessThanOrEqual, GreaterThan, GreaterThanOrEqual, EqualTo, NoEqualTo, IsZero, IsOne.

  • Elementary functions: Sqrt, Log, Log2, Log10, Log1PlusX, Exp, ExpM1, Exp2, Exp10, Hypot, Sin, Cos, Tan, Asin, Acos, Atan, Atan2, Sinh, Cosh, Tanh, Asinh, Acosh, Atanh.

  • Misc. operations: IsNaN, IsInfinity, FromInt32, FromDouble, Floor, Ceiling, IsReal.

To illustrate how to use these methods, we will implement the Newton-Raphson method for solving equations in generic arithmetic. The code below implements the algorithm:

C# 
public static double Solve(Func<double, double> function,
    Func<double, double> derivative, double initialGuess, double tolerance)
{
    int maxIterations = 20;
    double x = initialGuess;
    for (int i = 0; i < maxIterations; i++)
    {
        var y = function(x);
        var dy = derivative(x);
        var dx = y / dy;
        var x1 = x - dx;
        var activeTolerance = tolerance * Math.Abs(x1);
        if (Math.Abs(dx) < activeTolerance)
            return x1;
        x = x1;
    }
    throw new Exception("Max. number of iterations exceeded.");
}

To make this code generic, we add a generic type parameter to the method and replace all references to Double with this type parameter. We then replace all operators and method calls with the equivalent member of the Operations<T> class. This gives us:

C# 
public static T Solve<T>(Func<T, T> function, Func<T,T> derivative,
    T initialGuess, T tolerance) {

    int maxIterations = 20;
    T x = initialGuess;
    for (int i = 0; i < maxIterations; i++)
    {
        var y = function(x);
        var dy = derivative(x);
        var dx = Operations<T>.Divide(y, dy);
        var x1 = Operations<T>.Subtract(x, dx);
        var activeTolerance = Operations<T>.Multiply(tolerance,
            Operations<T>.Abs(x1));
        if (Operations<T>.LessThan(Operations<T>.Abs(dx), activeTolerance))
            return x1;
        x = x1;
    }
    throw new Exception("Max. number of iterations exceeded.");
}

We can now call this method with different types. We will keep it really simple here, and calculate the square root of 3 by solving the equation x2-3=0 :

C# 
var xSingle = Solve(x => x * x - 3.0f, x => 2.0f * x, 1.0f, 1e-5f);
var xDouble = Solve(x => x * x - 3.0, x => 2.0 * x, 1.0, 1e-15);
var xQuad = Solve(x => x * x - 3, x => 2 * x, Quad.One, 10 * Quad.Epsilon);
var xRational = Solve(x => x * x - 3, x => 2 * x, BigRational.One, 
    new BigRational(1, 1000000000));
Console.WriteLine("Single precision result: {0}", xSingle);
Console.WriteLine("Double precision result: {0}", xDouble);
Console.WriteLine("Quad precision result: {0}", xQuad);
Console.WriteLine("Rational result: {0} ({1})", xRational, (double)xRational);

Generic linear algebra

All types that represent vectors, matrices, or matrix decompositions are generic over the element type. This means that the same code can be used for different element types.

Not all methods are available for all operand types. For example, in order to solve a system of equations, you need at least full division. Integer division is not enough. Calling a method for which the operand type doesn't support the necessary operations will result in a GenericArithmeticException.

The example below constructs a 20x20 Hilbert matrix and computes its inverse. Hilbert matrices are infamous for the numerical difficulties involved in performing calculations with them. The condition number of a 20x20 Hilbert matrix is about 1030, which means that about 30 digits of precision are lost in the calculation. So, we use the BigRational type for the matrix elements.

C# 
int n = 20;
Matrix<BigRational> h = Matrix.CreateFromFunction(n, n,
        (i, j) => new BigRational(1, i + j + 1));
Matrix<BigRational> hInv = h.GetInverse();