BigNumbers QuickStart Sample (Visual Basic)

Illustrates the basic use of the arbitrary precision classes: BigInteger, BigRational, BigFloat in Visual Basic.

C# code F# code IronPython code Back to QuickStart Samples

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


Namespace Extreme.Numerics.QuickStart.VB
    ' Illustrates the use of the arbitrary precision number
    ' classes in the Extreme Optimization Mathematics Library for .NET.
    Module BigNumbers

        Sub main()
            ' In this QuickStart sample, we'll use 5 different methods to compute 
            ' the value of pi, the ratio of the circumference to the diameter of
            ' a circle.
            Dim pi As BigFloat

            ' The number of decimal digits.
            Dim digits As Integer = 10000
            ' The equivalent number of binary digits, to account for round-off error:
            Dim binaryDigits As Integer = CInt(8 + digits * Math.Log(10, 2))
            ' The number of digits in the last correction, if applicable.
            Dim correctionDigits As Double

            ' First, create an AccuracyGoal for the number of digits we want.
            ' We'll add 5 extra digits to account for round-off error.
            Dim goal As AccuracyGoal = AccuracyGoal.Absolute(digits + 5)
            Console.WriteLine("Calculating {0} digits of pi:", digits)

            ' Create a stopwatch so we can time the results.
            Dim sw As System.Diagnostics.Stopwatch = New System.Diagnostics.Stopwatch()

            '
            ' Method 1: Arctan formula
            ' 
            ' pi/4 = 88L(172) + 51L(239) + 32L(682) + 44L(5357) + 68L(12943)
            ' Where L(p) = Arctan(1/p)
            ' We will use big integer arithmetic for this.
            ' See the helper function Arctan later in this file.
            Console.WriteLine("Arctan formula using integer arithmetic:")
            sw.Start()
            Dim coefficients As Integer() = {88, 51, 32, 44, 68}
            Dim arguments As Integer() = {172, 239, 682, 5357, 12943}
            pi = BigFloat.Zero
            For k As Integer = 0 To 4
                pi = pi + coefficients(k) * Arctan(arguments(k), binaryDigits)
                Console.WriteLine("Step {0}: ({1:F3} seconds)", k + 1, sw.Elapsed.TotalSeconds)
            Next
            ' The ScaleByPowerOfTwo is the fastest way to multiply
            ' or divide by a power of two:
            pi = BigFloat.ScaleByPowerOfTwo(pi, 2)
            sw.Stop()
            Console.WriteLine("Total time: {0:F3} seconds.", sw.Elapsed.TotalSeconds, pi)
            Console.WriteLine()

            '
            ' Method 2: Rational approximation
            ' 
            ' pi/2 = 1 + 1/3 + (1*2)/(3*5) + (1*2*3)/(3*5*7) + ...
            '      = 1 + 1/3 * (1 + 2/5 * (1 + 3/7 * (1 + ...)))
            ' We gain 1 bit per iteration, so we know where to start.
            Console.WriteLine("Rational approximation using rational arithmetic.")
            Console.WriteLine("This is very inefficient, so we only do up to 10000 digits.")
            sw.Start()
            Dim an As BigRational = BigRational.Zero
            Dim n0 As Integer = binaryDigits
            If digits <= 10000 Then n0 = CInt(8 + 10000 * Math.Log(10, 2))
            For n As Integer = n0 To 1 Step -1
                an = New BigRational(n, 2 * n + 1) * an + 1
            Next
            pi = New BigFloat(2 * an, goal, RoundingMode.TowardsNearest)
            sw.Stop()
            Console.WriteLine("Total time: {0:F3} seconds.", sw.Elapsed.TotalSeconds, pi)
            Console.WriteLine()

            '
            ' Method 3: Arithmetic-Geometric mean
            '
            ' By Salamin & Brent, based on discoveries by C.F.Gauss.
            ' See http:'www.cs.miami.edu/~burt/manuscripts/gaussagm/agmagain-hyperref.pdf
            Console.WriteLine("Arithmetic-Geometric Mean:")
            sw.Reset()
            sw.Start()
            Dim x1 As BigFloat = BigFloat.Sqrt(2, goal, RoundingMode.TowardsNearest)
            Dim x2 As BigFloat = BigFloat.One
            Dim S As BigFloat = BigFloat.Zero
            Dim c As BigFloat = BigFloat.One
            For k As Integer = 0 To Integer.MaxValue
                S = S + BigFloat.ScaleByPowerOfTwo(c, k - 1)
                Dim aMean As BigFloat = BigFloat.ScaleByPowerOfTwo(x1 + x2, -1)
                Dim gMean As BigFloat = BigFloat.Sqrt(x1 * x2)
                x1 = aMean
                x2 = gMean
                c = (x1 + x2) * (x1 - x2)
                ' GetDecimalDigits returns the approximate number of digits in a number.
                ' A negative return value means the number is less than 1.
                correctionDigits = -c.GetDecimalDigits()
                Console.WriteLine("Iteration {0}: {1:F1} digits ({2:F3} seconds)", k, correctionDigits, sw.Elapsed.TotalSeconds)
                If (correctionDigits >= digits) Then Exit For
            Next
            pi = x1 * x1 / (1 - S)
            sw.Stop()
            Console.WriteLine("Total time: {0:F3} seconds.", sw.Elapsed.TotalSeconds, pi)
            Console.WriteLine()

            '
            ' Method 4: Borweins' quartic formula
            '
            ' This algorithm quadruples the number of correct digits
            ' in each iteration.
            ' See http:'en.wikipedia.org/wiki/Borwein's_algorithm
            Console.WriteLine("Quartic formula:")
            sw.Reset()
            sw.Start()
            Dim sqrt2 As BigFloat = BigFloat.Sqrt(2, goal, RoundingMode.TowardsNearest)
            Dim y As BigFloat = sqrt2 - BigFloat.One
            Dim a As BigFloat = New BigFloat(6, goal) - BigFloat.ScaleByPowerOfTwo(sqrt2, 2)
            Dim y2 As BigFloat = y * y
            Dim y3 As BigFloat
            Dim y4 As BigFloat = y2 * y2
            Dim da As BigFloat
            correctionDigits = 0
            For k As Integer = 0 To Integer.MaxValue
                If 4 * correctionDigits > digits Then Exit For
                Dim qrt As BigFloat = BigFloat.Root(1 - y4, 4)
                y = BigFloat.Subtract(1, qrt, goal, RoundingMode.TowardsNearest) _
                    / BigFloat.Add(1, qrt, goal, RoundingMode.TowardsNearest)
                ' y = BigFloat.Divide(1 - qrt, 1 + qrt, AccuracyGoal.InheritAbsolute, RoundingMode.TowardsNearest)
                y2 = y * y
                y3 = y * y2
                y4 = y2 * y2
                da = (BigFloat.ScaleByPowerOfTwo(y + y3, 2) + (6 * y2 + y4)) * a _
                    - BigFloat.ScaleByPowerOfTwo(y + y2 + y3, 2 * k + 1)
                da = da.RestrictPrecision(goal, RoundingMode.TowardsNearest)
                a += da
                correctionDigits = -da.GetDecimalDigits()
                Console.WriteLine("Iteration {0}: {1:F1} digits ({2:F3} seconds)", k, correctionDigits, sw.Elapsed.TotalSeconds)
            Next
            pi = BigFloat.Inverse(a)
            sw.Stop()
            Console.WriteLine("Total time: {0:F3} seconds.", sw.Elapsed.TotalSeconds, pi)
            Console.WriteLine()

            ' 
            ' Method 5: The built-in method
            '
            ' The method used to compute pi internally is an order of magnitude
            ' faster than any of the above.
            Console.WriteLine("Built-in function:")
            sw.Reset()
            sw.Start()
            pi = BigFloat.GetPi(goal)
            sw.Stop()
            Console.WriteLine("Total time: {0:F3} seconds.", sw.Elapsed.TotalSeconds, pi)
            ' The highest precision value of pi is cached, so
            ' getting pi to any precision up to that is super fast.
            Console.WriteLine("Built-in function (cached):")
            sw.Reset()
            sw.Start()
            pi = BigFloat.GetPi(goal)
            sw.Stop()
            Console.WriteLine("Total time: {0:F3} seconds.", sw.Elapsed.TotalSeconds, pi)

            Console.Write("Press Enter key to exit...")
            Console.ReadLine()
        End Sub

        ' Helper function to compute Arctan(1/p)
        '    p: The reciprocal of the argument.
        '    binaryDigits: The number of binary digits in the result.
        ' Returns; Arctan(1/p) to "binaryDigits" binary digits.
        Function Arctan(ByVal p As Integer, ByVal binaryDigits As Integer) As BigFloat

            ' We scale the result by a factor of 2^binaryDigits.
            ' The first term is 1/p.
            Dim power As BigInteger = BigInteger.Pow(2, binaryDigits) / p
            ' We store the sum in result.
            Dim result As BigInteger = power
            Dim subtract As Boolean = True
            Dim k As Integer = 0
            Do While (Not power.IsZero)
                k = k + 1
                ' Expressions involving big integers look exactly like any other arithmetic expression:

                ' The kth term is (-1)^k 1/(2k+1) 1/p^2k.
                ' So the power is 1/p^2 times the previous power.
                power = power / (p * p)
                ' And we alternately add and subtract
                If (subtract) Then
                    result -= power / (2 * k + 1)
                Else
                    result += power / (2 * k + 1)
                End If
                subtract = Not subtract
            Loop
            ' Scale the result.
            Return BigFloat.ScaleByPowerOfTwo(New BigFloat(result), -binaryDigits)
        End Function

    End Module

End Namespace