BigNumbers in Visual Basic QuickStart Sample
Illustrates the basic use of the arbitrary precision classes: BigInteger, BigRational, BigFloat in Visual Basic.
This sample is also available in: C#, F#, IronPython.
Overview
This QuickStart sample demonstrates how to use Numerics.NET’s arbitrary precision number types to perform high-precision mathematical calculations, specifically computing pi to thousands of digits.
The sample implements and compares five different methods for computing pi:
- Using an arctan formula with BigInteger arithmetic
- Using a rational approximation with BigRational arithmetic
- Using the arithmetic-geometric mean algorithm with BigFloat arithmetic
- Using Borweins’ quartic formula with BigFloat arithmetic
- Using the built-in BigFloat.GetPi method
Each method demonstrates different aspects of arbitrary precision arithmetic and shows how to work with the various number types. The sample includes performance timing for each approach and shows how to control precision and rounding modes. It also demonstrates important concepts like:
- Converting between different arbitrary precision types
- Controlling precision and accuracy goals
- Handling performance considerations
- Using mathematical functions with arbitrary precision
- Working with rational numbers and continued fractions
- Understanding convergence rates of different algorithms
The code
Option Infer On
' We'll also need the big number types.
Imports Numerics.NET
' Illustrates the use of the arbitrary precision number
' classes in 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 {digits} digits of pi:")
' 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