Using 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, 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<T> class. This gives us:
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:
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.
int n = 20;
Matrix<BigRational> h = Matrix.Create(n, n,
(i, j) => new BigRational(1, i + j + 1));
Matrix<BigRational> hInv = h.GetInverse();Operand Types and Arithmetic Types
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<T>, which means it supports addition, multiplication and related operations, but not division.
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.
Generic calculations
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<T, TResult> with its generic equivalent, Func<T, TResult>. This gives us the following:
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<T>. The different generic arithmetic interfaces are covered in the next section. This gives us the following final result:
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);