BigNumbers in Visual Basic QuickStart Sample
Illustrates the basic use of the arbitrary precision classes: BigInteger, BigRational, BigFloat in Visual Basic.
View this sample in: C# F# IronPython
Option Infer On
' We'll also need the big number types.
Imports Extreme.Mathematics
' Illustrates the use of the arbitrary precision number
' classes in Extreme Numerics.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 New 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 = BigRational.Zero
Dim n0 = binaryDigits
If digits <= 10000 Then n0 = CInt(8 + 10000 * Math.Log(10, 2))
For n = 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 = BigFloat.Sqrt(2, goal, RoundingMode.TowardsNearest)
Dim x2 = BigFloat.One
Dim S = BigFloat.Zero
Dim c = BigFloat.One
For k = 0 To Integer.MaxValue
S = S + BigFloat.ScaleByPowerOfTwo(c, k - 1)
Dim aMean = BigFloat.ScaleByPowerOfTwo(x1 + x2, -1)
Dim gMean = 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 = sqrt2 - BigFloat.One
Dim a = New BigFloat(6, goal) - BigFloat.ScaleByPowerOfTwo(sqrt2, 2)
Dim y2 = y * y
Dim y3 As BigFloat
Dim y4 = 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 = 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(p As Integer, binaryDigits As Integer) As BigFloat
' We scale the result by a factor of 2^binaryDigits.
' The first term is 1/p.
Dim power = BigInteger.Pow(2, binaryDigits) / p
' We store the sum in result.
Dim result = power
Dim subtract = True
Dim k = 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