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 OperationsT
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 OperationsT
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.");
}Public Shared Function Solve(f As Func(Of Double, Double),
df As Func(Of Double, Double), initialGuess As Double, tolerance As Double)
Dim maxIterations As Integer = 20
Dim x As Double = initialGuess
For i = 1 To maxIterations
Dim y = f(x)
Dim dy = df(x)
Dim dx = y / dy
Dim x1 = x - dx
Dim activeTolerance = tolerance * Math.Abs(x1)
If (Math.Abs(dx) < activeTolerance) Then
Return x1
End If
x = x1
Next
Throw New Exception("Max. number of iterations exceeded.")
End FunctionNo code example is currently available or this language may not be supported.
let SolveDouble f df initialGuess tolerance =
let maxIterations = 20
let rec iterate x i =
match i with
| _ when i < maxIterations ->
let y = f x
let dy = df x
let dx = y / dy
let x1 = x - dx
let activeTolerance = tolerance * abs(x1)
if (abs(dx) < activeTolerance) then x1 else iterate x1 (i+1)
| _ -> failwith "Max. number of iterations exceeded."
iterate initialGuess 0
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
OperationsT
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.");
}Public Shared Function Solve(Of T)(f As Func(Of T, T), df As Func(Of T, T),
initialGuess As T, tolerance As T) As T
Dim maxIterations As Integer = 20
Dim x As T = initialGuess
For i = 1 To maxIterations
Dim y = f(x)
Dim dy = df(x)
Dim dx = Operations(Of T).Divide(y, dy)
Dim x1 = Operations(Of T).Subtract(x, dx)
Dim activeTolerance = Operations(Of T).Multiply(tolerance,
Operations(Of T).Abs(x1))
If (Operations(Of T).LessThan(Operations(Of T).Abs(dx), activeTolerance)) Then
Return x1
End If
x = x1
Next
Throw New Exception("Max. number of iterations exceeded.")
End FunctionNo code example is currently available or this language may not be supported.
let Solve f df (initialGuess : 'a) tolerance =
let maxIterations = 20
let rec iterate x i =
match i with
| _ when i < maxIterations ->
let y = f(x)
let dy = df(x)
let dx = Operations<_>.Divide.Invoke(y, dy)
let x1 = Operations<_>.Subtract.Invoke(x, dx)
let activeTolerance =
Operations<_>.Multiply.Invoke(tolerance,
Operations<_>.Abs.Invoke(x1))
if (Operations<_>.LessThan.Invoke(Operations<_>.Abs.Invoke(dx), activeTolerance))
then x1
else
iterate x1 (i+1)
| _ -> failwith "Max. number of iterations exceeded."
iterate initialGuess 0
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);Dim xSingle = Solve(Function(x) x * x - 3.0F, Function(x) 2.0F * x, 1.0F, 0.00001F)
Dim xDouble = Solve(Function(x) x * x - 3.0, Function(x) 2.0 * x, 1.0, 0.000000000000001)
Dim xQuad = Solve(Function(x) x * x - 3, Function(x) 2 * x, Quad.One, 10 * Quad.Epsilon)
Dim xRational = Solve(Function(x) x * x - 3, Function(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, CType(xRational, Double))No code example is currently available or this language may not be supported.
let xSingle = Solve (fun x -> x * x - 3.0f) (fun x -> 2.0f * x) 1.0f 1e-5f
let xDouble = Solve (fun x -> x * x - 3.0) (fun x -> 2.0 * x) 1.0 1e-15
let xQuad = Solve (fun x -> x * x - Quad(3)) (fun x -> Quad(2) * x) Quad.One (Quad(10) * Quad.Epsilon)
let xRational = Solve(fun x -> x * x - 3) (fun x -> 2 * x) BigRational.One (BigRational(1, 1000000000))
printfn "Single precision result: %f" xSingle
printfn "Double precision result: %f" xDouble
printfn "Quad precision result: %A" xQuad
printfn "Rational result: %A (%f)" 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
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();Dim n As Integer = 20
Dim h As Matrix(Of BigRational) = Matrix.Create(n, n,
Function(i, j) New BigRational(1, i + j + 1))
Dim hInv = h.GetInverse()No code example is currently available or this language may not be supported.
let n = 20
let h = Matrix.Create(n, n,
fun i j -> BigRational(1, i + j + 1))
let 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
IRingOperationsT,
which means it supports addition, multiplication and related operations, but not division.
IRingOperations<T> R = TypeAssociationRegistry.GetInstance(typeof(T),
TypeAssociationRegistry.ArithmeticKey) as IRingOperations<T>;Dim R As IRingOperations(Of T) = CType(TypeAssociationRegistry.GetInstance(GetType(T),
TypeAssociationRegistry.ArithmeticKey), IRingOperations(Of T));No code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.
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;
}
}Class Solver
' Member fields:
Dim f, df As RealFunction
' The function to solve:
Public Property TargetFunction() As RealFunction
Get
Return f
End Get
Set(ByVal value As RealFunction)
f = value
End Set
End Property
' The derivative of the function to solve.
Public Property DerivativeOfTargetFunction() As RealFunction
Get
Return df
End Get
Set(ByVal value As RealFunction)
df = value
End Set
End Property
' The core algorithm.
Public Function Solve(ByVal initialGuess As Double, ByVal tolerance As Double) As Double
Dim iterations As Integer = 0
Dim x As Double = initialGuess
Dim dx As Double = 0.0
Do
iterations = iterations + 1
' Compute the denominator of the correction term.
Dim dfx As Double = df(x)
If (dfx = 0) Then
' Change value by 2x tolerance.
dx = 2 * tolerance
Else
dx = f(x) / dfx
End If
x = x - dx
If Math.Abs(dx) < tolerance Then
Return x
End If
Loop While (iterations < MaxIterations)
Return x
End Function
End ClassNo code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.
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
FuncT, TResult
with its generic equivalent,
FuncT, 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;
}
}Class Solver(Of T)
' Member fields:
Dim f, df As GenericFunction(Of T)
' The function to solve:
Public Property TargetFunction() As GenericFunction(Of T)
Get
Return f
End Get
Set(ByVal value As GenericFunction(Of T))
f = value
End Set
End Property
' The derivative of the function to solve.
Public Property DerivativeOfTargetFunction() As GenericFunction(Of T)
Get
Return df
End Get
Set(ByVal value As GenericFunction(Of T))
df = value
End Set
End Property
' The core algorithm.
Public Function Solve(ByVal initialGuess As T, ByVal tolerance As T) As T
Dim iterations As Integer = 0
Dim x As T = initialGuess
Dim dx As T = 0.0
Do
iterations = iterations + 1
Dim dfx As T = df(x)
If (dfx = 0) Then
dx = 2 * tolerance
Else
dx = f(x) / dfx
End If
x = x - dx
If Math.Abs(dx) < tolerance Then
Return x
End If
Loop While (iterations < MaxIterations)
Return x
End Function
End ClassNo code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.
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 IFieldOperationsT.
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;
}
}Class Solver(Of T)
Dim ops As IFieldOperations(Of T) = _
CType(TypeAssociationRegistry.GetInstance( _
GetType(T), TypeAssociationRegistry.ArithmeticKey), _
IFieldOperations(Of T))
' Member fields:
Dim f, df As GenericFunction(Of T)
' The function to solve:
Public Property TargetFunction() As GenericFunction(Of T)
Get
Return f
End Get
Set(ByVal value As GenericFunction(Of T))
f = value
End Set
End Property
' The derivative of the function to solve.
Public Property DerivativeOfTargetFunction() As GenericFunction(Of T)
Get
Return df
End Get
Set(ByVal value As GenericFunction(Of T))
df = value
End Set
End Property
' The core algorithm.
Public Function Solve(ByVal initialGuess As T, ByVal tolerance As T) As T
Dim iterations As Integer = 0
Dim x As T = initialGuess
Dim dx As T = ops.Zero
Do
iterations = iterations + 1
Dim dfx As T = df(x)
If (ops.Compare(dfx, ops.Zero) = 0) Then
dx = ops.ScaleByPowerOfTwo(tolerance, 1)
Else
dx = ops.Divide(f(x), dfx)
End If
x = ops.Subtract(x, dx)
If ops.Compare(ops.Abs(dx), tolerance) < 0 Then
Return x
End If
Loop While (iterations < MaxIterations)
Return x
End Function
End ClassNo code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.
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));Function fBigFloat(ByVal x As BigFloat) As BigFloat
Return BigFloat.Sin(x)
End Function
Function dfBigFloat(ByVal x As BigFloat) As BigFloat
Return BigFloat.Cos(x)
End Function
Dim bigFloatSolver As New Solver(Of BigFloat)
' Set the function to solve, and its derivative.
bigFloatSolver.TargetFunction = New GenericFunction(Of BigFloat)(AddressOf fBigFloat)
bigFloatSolver.DerivativeOfTargetFunction = New GenericFunction(Of BigFloat)(AddressOf dfBigFloat)
' Now solve to within a tolerance of 10^-100.
Dim sqrt2 As BigFloat = bigFloatSolver.Solve(1, BigFloat.Pow(10, -100))No code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.
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));Function fBigRational(ByVal x As BigRational) As BigRational
Return x * x - 2
End Function
Function dfBigRational(ByVal x As BigRational) As BigRational
Return 2 * x
End Function
Dim bigRationalSolver As New Solver(Of BigRational)()
' Set the function to solve, and its derivative.
bigRationalSolver.TargetFunction = New GenericFunction(Of BigRational)(AddressOf fBigRational)
bigRationalSolver.DerivativeOfTargetFunction = New GenericFunction(Of BigRational)(AddressOf dfBigRational)
' Now solve to within a tolerance of 10^-100.
Dim sqrt2 As BigRational = bigRationalSolver.Solve(1, BigRational.Pow(10, -100))No code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.
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);Function fDouble(ByVal x As double) As double
Return x * x - 2
End Function
Function dfDouble(ByVal x As double) As double
Return 2 * x
End Function
Dim doubleSolver As New Solver(Of Double)()
' Set the function to solve, and its derivative.
doubleSolver.TargetFunction = New GenericFunction(Of Double)(AddressOf fDouble)
doubleSolver.DerivativeOfTargetFunction = New GenericFunction(Of Double)(AddressOf dfDouble)
' Now solve to within a tolerance of 10^-10.
Dim sqrt2 As Double = doubleSolver.Solve(1.0, 0.0000000001)No code example is currently available or this language may not be supported.
No code example is currently available or this language may not be supported.