New Version 6.0!

Try it for free with our fully functional 60-day trial version.

Download now!

QuickStart Samples

Generic Algorithms QuickStart Sample (F#)

Illustrates how to write algorithms that are generic over the numerical type of the arguments in F#.

C# code Visual Basic code Back to QuickStart Samples

#light

open System

// Basic generic types live in Extreme.Mathematics.Generics.
open Extreme.Mathematics.Generic
// We'll also need the big number types.
open Extreme.Mathematics

// Illustrates writing generic algorithms that can be 
// applied to different operand types using the types in the 
// Extreme.Mathematics.Generic namespace.

// We will implement a simple Newton-Raphson solver class.
// The code for the solver is below.

// Here we will call the generic solver with three
// different operand types: BigFloat, BigRational and Double.

let f a b = a - b

// Class that contains the generic Newton-Raphson algorithm.
type Solver<'a>(f : 'a -> 'a, df : 'a -> 'a) =

    // Use a static variable to hold the arithmetic instance.
    // We get the instance from the global TypeAssociationRegistry:
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IFieldOperations<'a>

    let mutable maxIterations = 100
    
    // The core algorithm.
    // Arithmetic operations are replaced by calls to
    // methods on the arithmetic object (ops).
    let rec solveRec x0 tolerance maxIterations iterations =
        // Compute the denominator of the correction term.
        let dfx = df(x0)
        // Relational operators map to the Compare method.
        // We also use the value of zero for the operand type.
        // if (dfx == 0)
        let dx =
            if (ops.Compare(dfx, ops.Zero) = 0) then
                // Change value by 2x tolerance.
                // When multiplying by a power of two, it's more efficient 
                // to use the ScaleByPowerOfTwo method.
                ops.ScaleByPowerOfTwo(tolerance, 1)
            else
                // dx = f(x) / df(x)
                let fx = f(x0)
                ops.Divide(fx, dfx)

        // x -= dx
        let x = ops.Subtract(x0, dx)
        
        // if |dx|^2<tolerance
        // Convergence is quadratic (in most cases), so we should be good here:
        if (ops.Compare(ops.Multiply(dx,dx), tolerance) < 0 || iterations > maxIterations) then
            x
        else 
            solveRec x tolerance maxIterations (iterations+1)

    let solve initialGuess tolerance maxIterations =
        solveRec initialGuess tolerance maxIterations 0

    // The maximum number of iterations.
    member this.MaxIterations
        with get() = maxIterations
        and set(value) = maxIterations <- value

    member this.Solve initialGuess tolerance = solve initialGuess tolerance this.MaxIterations

// First, let's compute pi to 100 digits
// by solving the equation sin(x) == 0 with
// an initual guess of 3.
printfn "Computing pi by solving sin(x) == 0 with x0 = 3 using BigFloat:"
// Create the solver object.
let bigFloatSolver = 
    let f = fun (x : BigFloat) -> sin(x)
    let df = fun (x : BigFloat) -> cos(x)
    Solver(f, df)
// Now solve to within a tolerance of 10^-100.
let pi = bigFloatSolver.Solve (BigFloat 3) (BigFloat.Pow(BigFloat 10, -100))
// Print the results...
printfn "Computed value: %s" (pi.ToString("F100"))
// and verify:
printfn "Known value:    %s" (BigFloat.GetPi(AccuracyGoal.Absolute(100.0)).ToString("F100"))
printfn ""

// Next, we will use rational numbers to compute
// an approximation to the square root of 2.
printfn "Computing sqrt(2) by solving x^2 == 2 using BigRational:"
// Create the solver...
let bigRationalSolver = 
    let f = fun (x : BigRational) -> x ** 2 - 2
    let df = fun (x : BigRational) -> 2 * x
    Solver(f, df)
// Compute the solution...
let sqrt2 = bigRationalSolver.Solve BigRational.One (BigRational.Pow(BigRational(10,1), -100))
// And print the result.
printfn "Rational approximation: %A" sqrt2
// To verify, we convert the BigRational to a BigFloat:
printfn "As real number: %s" (BigFloat(sqrt2, AccuracyGoal.Absolute(100.0), RoundingMode.TowardsNearest).ToString("F100"))
printfn "Known value:    %s" (BigFloat.Sqrt(BigFloat 2, AccuracyGoal.Absolute(100.0), RoundingMode.TowardsNearest).ToString("F100"))
printfn ""

// Now, we compute the Lambert W function at x = 3.
printfn "Computing Lambert's W at x = 3 by solving x*exp(x) == 3 using double solver:"
// Create the solver...
let doubleSolver = 
    let f = fun x -> x * exp(x) - 3.0
    let df = fun x -> (1.0 + x) * exp(x)
    Solver(f, df)
// Compute the solution...
let W3 = doubleSolver.Solve 1.0 1e-15
// And print the result.
printfn "Solution:    %A" W3
printfn "Known value: %A" (Elementary.LambertW(3.0))



// Finally, we use generic functions:
printfn "Using generic function delegates:"

// We can define some inline helper functions:
let inline (+) (a : 'a) b =
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IGroupOperations<'a>
    ops.Add(a, b)

let inline (-) (a : 'a) b =
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IGroupOperations<'a>
    ops.Subtract(a, b)

let inline (*) (a : 'a) b =
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IRingOperations<'a>
    ops.Multiply(a, b)

let inline one<'a> =
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IRingOperations<'a>
    ops.One

let inline zero<'a> =
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IGroupOperations<'a>
    ops.Zero

let inline exp (x : 'a) =
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IRealOperations<'a>
    ops.Exp(x)

let inline ofint n =
    let ops = 
        TypeAssociationRegistry.GetInstance(typeof<'a>, TypeAssociationRegistry.ArithmeticKey) 
            :?> IRingOperations<'a>
    ops.FromInt32(n)

// Using these definitions, we can use standard notation:
let fGeneric x = x * exp x - ofint 3
let dfGeneric x = exp x * (x + one<_>)

let genericDoubleSolver = Solver<float>(fGeneric, dfGeneric)
let genericW3 = doubleSolver.Solve 1.0 1e-15
printfn "Double:      %A" genericW3

let genericBigFloatSolver = Solver<BigFloat>(fGeneric, dfGeneric)
let bigW3 = genericBigFloatSolver.Solve BigFloat.One (BigFloat.Pow(BigFloat 10, -100))
printfn "BigFloat:    %s" (bigW3.ToString("F100"))

printf "Press Enter key to exit..."
Console.ReadLine() |> ignore