Using Generic Arithmetic | Extreme Optimization Numerical Libraries for .NET Professional |

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.

The Operations

Using the Operations

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, Log1PlusX, Exp, 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:

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

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
x^{2}-3=0:

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

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
10^{30}, which means that about 30 digits of precision
are lost in the calculation. So, we use the
BigRational type
for the matrix elements.

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

Generic arithmetic is implemented by associating with each operand type a class (the arithmetic type) that implements the operations on that operand type. The association between an operand type and its arithmetic type is declared using a TypeAssociationAttribute.

To perform generic operations, you need to first acquire an instance of the arithmetic type. This is done with a call to the GetInstance method of the TypeAssociationRegistry class. This is a static method, and takes two arguments. The first argument is the operand type. The second argument is a key that identifies the association. The key is a string. The key that identifies the arithmetic type is ArithmeticKey.

The return value is an instance of the arithmetic type. It implements one or more of
the generic interfaces defined in the next section. The example below illustrates how to
obtain the arithmetic for a generic type. The arithmetic implements
IRingOperations

IRingOperations<T> R = TypeAssociationRegistry.GetInstance(typeof(T), TypeAssociationRegistry.ArithmeticKey) as IRingOperations<T>;

Obtaining the arithmetic is a relatively expensive operation, especially compared to the work done by individual operations. It is usually a good idea to cache the arithmetic in a static variable that is available to all instances of the containing type.

Once you have an instance of the arithmetic type, you can use it to perform calculations with the operand type. Not all operand types support the same set of operations. Which operations they support is determined by the interfaces that the arithmetic type implements. A thorough overview of the different interfaces for generic arithmetic is given in the next section.

To illustrate how generic algorithms can be coded, we will implement a generic algorithm for the solution of an equation using Newton's method. Newton's method finds the solution of an equation f(x)=0 using the following iteration formula:

We start out with the non-generic version using doubles.

class Solver { // Member fields: RealFunction f, df; // The function to solve: public RealFunction TargetFunction { get { return f; } set { f = value; } } // The derivative of the function to solve. public RealFunction DerivativeOfTargetFunction { get { return df; } set { df = value; } } // The core algorithm. public double Solve(double initialGuess, double tolerance) { int iterations = 0; double x = initialGuess; double dx = 0.0; do { iterations++; // Compute the denominator of the correction term. double dfx = df(x); if (dfx == 0.0) { // Change value by 2x tolerance. dx = 2 * tolerance; } else { dx = f(x) / dfx; } x -= dx; if (Math.Abs(dx) < tolerance) return x; } while (iterations < 100); return x; } }

The first step is to add a generic type parameter to the class definition,
and replace all occurrences of the operand type
with the generic parameter type. We also want to replace all related types, like delegates
and lists, with the corresponding generic type. In this case, we replace all occurrences
of Double
with T.
We also replace
Func

class Solver<T> { // Member fields: GenericFunction<T> f, df; // The function to solve: public GenericFunction<T> TargetFunction { get { return f; } set { f = value; } } // The derivative of the function to solve. public GenericFunction<T> DerivativeOfTargetFunction { get { return df; } set { df = value; } } // The core algorithm. public T Solve(T initialGuess, T tolerance) { int iterations = 0; T x = initialGuess; T dx = 0.0; do { iterations++; // Compute the denominator of the correction term. T dfx = df(x); if (dfx == 0.0) { // Change value by 2x tolerance. dx = 2 * tolerance; } else { dx = f(x) / dfx; } x -= dx; if (Math.Abs(dx) < tolerance) return x; } while (iterations < 100); return x; } }

The last step is to create a static (Shared in Visual Basic)
variable in the class to hold an instance of the arithmetic type, and to
replace all arithmetic operations with calls to methods of this instance.
For Newton's algorithm, we need division, so the type of the arithmetic should
be IFieldOperations

class Solver<T> { static IFieldOperations<T> ops = (IFieldOperations<T>)TypeAssociationRegistry.GetInstance( typeof(T), TypeAssociationRegistry.ArithmeticKey); // Member fields: GenericFunction<T> f, df; // The function to solve: public GenericFunction<T> TargetFunction { get { return f; } set { f = value; } } // The derivative of the function to solve. public GenericFunction<T> DerivativeOfTargetFunction { get { return df; } set { df = value; } } // The core algorithm. public T Solve(T initialGuess, T tolerance) { int iterations = 0; T x = initialGuess; T dx = ops.Zero; do { iterations++; T dfx = df(x); if (ops.Compare(dfx, ops.Zero) == 0) { // Change value by 2x tolerance. dx = ops.ScaleByPowerOfTwo(tolerance, 1); } else { dx = ops.Divide(f(x), dfx); } x = ops.Subtract(x, dx); if (ops.Compare(ops.Abs(dx), tolerance) < 0) return x; } while (iterations < 100); return x; } }

Values of 0 were replaced with calls to the Zero property. Relational operators were replaced with calls to the Compare method.

We are now ready to use this class with any operand type that has an associated arithmetic with a full division operator. As a first example, we'll compute the number π to 100 digits by solving sin x = 0 using the BigFloat with a starting value of 3:

Solver<BigFloat> bigFloatSolver = new Solver<BigFloat>(); // Set the function to solve, and its derivative. bigFloatSolver.TargetFunction = delegate(BigFloat x) { return BigFloat.Sin(x); }; bigFloatSolver.DerivativeOfTargetFunction = delegate(BigFloat x) { return BigFloat.Cos(x); }; // Now solve to within a tolerance of 10^-100. BigFloat pi = bigFloatSolver.Solve(3, BigFloat.Pow(10, -100));

Next, we compute a rational approximation to the square root of 2:

Solver<BigRational> bigRationalSolver = new Solver<BigRational>(); // Set the function to solve, and its derivative. bigRationalSolver.TargetFunction = delegate(BigRational x) { return x*x - 2; }; bigRationalSolver.DerivativeOfTargetFunction = delegate(BigRational x) { return 2*x; }; // Now solve to within a tolerance of 10^-100. BigRational pi = bigRationalSolver.Solve(3, BigRational.Pow(10, -100));

Finally, we show that our solver class works just as well for the built-in types, including Double. We use the exact same code as above, and simply replace all references to BigRational with references to Double.

Solver<double> doubleSolver = new Solver<double>(); // Set the function to solve, and its derivative. doubleSolver.TargetFunction = delegate(double x) { return x*x - 2; }; doubleSolver.DerivativeOfTargetFunction = delegate(double x) { return 2*x; }; // Now solve to within a tolerance of 10^-100. double pi = doubleSolver.Solve(3.0, 1e-10);

Copyright © 2004-20116,
Extreme Optimization. All rights reserved.

*Extreme Optimization,* *Complexity made simple*, *M#*, and *M
Sharp* are trademarks of ExoAnalytics Inc.

*Microsoft*, *Visual C#, Visual Basic, Visual Studio*, *Visual
Studio.NET*, and the *Optimized for Visual Studio* logo

are
registered trademarks of Microsoft Corporation.