New Version 6.0!

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

Download now!

QuickStart Samples

Generic Algorithms QuickStart Sample (C#)

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

Visual Basic code F# code Back to QuickStart Samples

using System;

namespace Extreme.Numerics.QuickStart.CSharp
{
    // Basic generic types live in Extreme.Mathematics.Generics.
    using Extreme.Mathematics.Generic;
    // We'll also need the big number types.
    using Extreme.Mathematics;

    /// <summary>
    /// Illustrates writing generic algorithms that can be 
    /// applied to different operand types using the types in the 
    /// Extreme.Mathematics.Generic namespace.
    /// </summary>
    class GenericAlgorithms
    {
        static void Main(string[] args)
        {
            // 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.

            // First, let's compute pi to 100 digits
            // by solving the equation sin(x) == 0 with
            // an initual guess of 3.
            Console.WriteLine("Computing pi by solving sin(x) == 0 with x0 = 3 using BigFloat:");
            // Create the solver object.
            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));
            // Print the results...
            Console.WriteLine("Computed value: {0:F100}", pi);
            // and verify:
            Console.WriteLine("Known value:    {0:F100}", 
                BigFloat.GetPi(AccuracyGoal.Absolute(100)));
            Console.WriteLine();

            // Next, we will use rational numbers to compute
            // an approximation to the square root of 2.
            Console.WriteLine("Computing sqrt(2) by solving x^2 == 2 using BigRational:");
            // Create the solver...
            Solver<BigRational> bigRationalSolver = new Solver<BigRational>();
            // Set properties...
            bigRationalSolver.TargetFunction = delegate(BigRational x) { return x * x - 2; };
            bigRationalSolver.DerivativeOfTargetFunction = delegate(BigRational x) { return 2 * x; };
            // Compute the solution...
            BigRational sqrt2 = bigRationalSolver.Solve(1, BigRational.Pow(10, -100));
            // And print the result.
            Console.WriteLine("Rational approximation: {0}", sqrt2);
            // To verify, we convert the BigRational to a BigFloat:
            Console.WriteLine("As real number: {0:F100}", 
                new BigFloat(sqrt2, AccuracyGoal.Absolute(100), RoundingMode.TowardsNearest));
            Console.WriteLine("Known value:    {0:F100}", 
                BigFloat.Sqrt(2, AccuracyGoal.Absolute(100), RoundingMode.TowardsNearest));
            Console.WriteLine();

            // Now, we compute the Lambert W function at x = 3.
            Console.WriteLine("Computing Lambert's W at x = 3 by solving x*exp(x) == 3 using double solver:");
            // Create the solver...
            Solver<double> doubleSolver = new Solver<double>();
            // Set properties...
            doubleSolver.TargetFunction = delegate(double x) { return x * Math.Exp(x) - 3; };
            doubleSolver.DerivativeOfTargetFunction = delegate(double x) { return Math.Exp(x) * (1 + x); };
            // Compute the solution...
            double W3 = doubleSolver.Solve(1.0, 1e-15);
            // And print the result.
            Console.WriteLine("Solution:    {0}", W3);
            Console.WriteLine("Known value: {0}", 
                Extreme.Mathematics.Elementary.LambertW(3.0));

            // Finally, we use generic functions:
            Console.WriteLine("Using generic function delegates:");
            // The delegates are slightly more complicated.
            doubleSolver.TargetFunction = fGeneric<double>;
            doubleSolver.DerivativeOfTargetFunction = dfGeneric<double>;
            double genericW3 = doubleSolver.Solve(1, 1e-15);
            Console.WriteLine("Double:      {0}", genericW3);
            bigFloatSolver.TargetFunction = fGeneric<BigFloat>;
            bigFloatSolver.DerivativeOfTargetFunction = dfGeneric<BigFloat>;
            BigFloat bigW3 = bigFloatSolver.Solve(1, BigFloat.Pow(10, -100));
            Console.WriteLine("BigFloat:    {0:F100}", bigW3);

            Console.Write("Press Enter key to exit...");
            Console.ReadLine();
        }

        // Generic versions of the above
        static T fGeneric<T>(T x)
        {
            IRealOperations<T> ops =
                (IRealOperations<T>)TypeAssociationRegistry.GetInstance(
                    typeof(T), TypeAssociationRegistry.ArithmeticKey);
            return ops.Subtract(
                ops.Multiply(x, ops.Exp(x)),
                ops.FromInt32(3));
        }
        static T dfGeneric<T>(T x)
        {
            IRealOperations<T> ops =
                (IRealOperations<T>)TypeAssociationRegistry.GetInstance(
                    typeof(T), TypeAssociationRegistry.ArithmeticKey);
            return ops.Multiply(
                ops.Exp(x),
                ops.Add(x, ops.One));
        }
    }

    /// <summary>
    /// Class that contains the generic Newton-Raphson algorithm.
    /// </summary>
    /// <typeparam name="T">The operand type.</typeparam>
    class Solver<T>
    {
        // Use a static variable to hold the arithmetic instance.
        // We get the instance from the global TypeAssociationRegistry:
        static IFieldOperations<T> ops = 
            (IFieldOperations<T>)TypeAssociationRegistry.GetInstance(
                typeof(T), TypeAssociationRegistry.ArithmeticKey);

        // Member fields:
        Func<T,T> f, df;
        int maxIterations = 100;

        // The function to solve:
        public Func<T,T> TargetFunction
        {
            get { return f; }
            set { f = value; }
        }
        // The derivative of the function to solve.
        public Func<T,T> DerivativeOfTargetFunction
        {
            get { return df; }
            set { df = value; }
        }
        // The maximum number of iterations.
        public int MaxIterations
        {
            get { return maxIterations; }
            set { maxIterations = value; }
        }

        // The core algorithm.
        // Arithmetic operations are replaced by calls to
        // methods on the arithmetic object (ops).
        public T Solve(T initialGuess, T tolerance)
        {
            int iterations = 0;

            T x = initialGuess;
            T dx = ops.Zero;
            do
            {
                iterations++;
                // Compute the denominator of the correction term.
                T dfx = df(x);
                // Relational operators map to the Compare method.
                // We also use the value of zero for the operand type.
                // if (dfx == 0)
                if (ops.Compare(dfx, ops.Zero) == 0)
                {
                    // Change value by 2x tolerance.
                    // When multiplying by a power of two, it's more efficient 
                    // to use the ScaleByPowerOfTwo method.
                    dx = ops.ScaleByPowerOfTwo(tolerance, 1);
                }
                else
                {
                    // dx = f(x) / df(x)
                    dx = ops.Divide(f(x), dfx);
                }
                // x -= dx;
                x = ops.Subtract(x, 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)
                    return x;
            }
            while (iterations < MaxIterations);
            return x;
        }
    }
}