BigNumbers in IronPython QuickStart Sample
Illustrates the basic use of the arbitrary precision classes: BigInteger, BigRational, BigFloat in IronPython.
View this sample in: C# Visual Basic F#
```Python
import numerics
from math import log
# The arbitrary precision types reside in the Extreme.Mathematics
# namespace.
from Extreme.Mathematics import *
# Illustrates the use of the arbitrary precision number
# classes in Numerics.NET.
# 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.
# The number of decimal digits.
digits = 10000
# The equivalent number of binary digits, to account for round-off error:
binaryDigits = (int)(8 + digits * log(10, 2))
# The number of digits in the last correction, if applicable.
# correctionDigits
# First, create an AccuracyGoal for the number of digits we want.
# We'll add 5 extra digits to account for round-off error.
goal = AccuracyGoal.Absolute(digits + 5)
print "Calculating {0} digits of pi:".format(digits)
# Create a stopwatch so we can time the results.
from System.Diagnostics import Stopwatch
sw = 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.
# Helper function to compute Arctan(1/p)
# <param name="p">The reciprocal of the argument.</param>
# <param name="binaryDigits">The number of binary digits in the result.</param>
# <returns>Arctan(1/<paramref name="p"/>) to <paramref name="binaryDigits"/> binary digits.</returns>
def Arctan(p, binaryDigits):
# We scale the result by a factor of 2^binaryDigits.
# The first term is 1/p.
power = BigInteger.Pow(2, binaryDigits) / p
# We store the sum in result.
result = power
subtract = True
k = 0
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 /= p ** 2
# And we alternately add and subtract
if subtract:
result -= power / (2 * k + 1)
else:
result += power / (2 * k + 1)
subtract = not subtract
# Scale the result.
return BigFloat.ScaleByPowerOfTwo(BigFloat(result), -binaryDigits)
print "Arctan formula using integer arithmetic:"
sw.Start()
coefficients = [ 88, 51, 32, 44, 68 ]
arguments = [ 172, 239, 682, 5357, 12943 ]
pi = BigFloat.Zero
for k in range(0, 5):
pi += coefficients[k] * Arctan(arguments[k], binaryDigits)
print "Step {0}: ({1:.3f} seconds)".format(k + 1, sw.Elapsed.TotalSeconds)
# The ScaleByPowerOfTwo is the fastest way to multiply
# or divide by a power of two:
pi = BigFloat.ScaleByPowerOfTwo(pi, 2)
sw.Stop()
print "Total time: {0:.3f} seconds.".format(sw.Elapsed.TotalSeconds, pi)
print
#
# 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.
print "Rational approximation using rational arithmetic."
print "This is very inefficient, so we only do up to 10000 digits."
sw.Start()
an = BigRational.Zero
n0 = binaryDigits if digits <= 10000 else (int)(8 + 10000 * log(10, 2))
for n in range(n0,0,-1):
an = BigRational(n, 2 * n + 1) * an + 1
pi = BigFloat(2 * an, goal, RoundingMode.TowardsNearest)
sw.Stop()
print "Total time: {0:.3f} seconds.".format(sw.Elapsed.TotalSeconds, pi)
print
#
# 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
print "Arithmetic-Geometric Mean:"
sw.Reset()
sw.Start()
x1 = BigFloat.Sqrt(2, goal, RoundingMode.TowardsNearest)
x2 = BigFloat.One
S = BigFloat.Zero
c = BigFloat.One
for k in xrange(0, int.MaxValue):
S += BigFloat.ScaleByPowerOfTwo(c, k - 1)
aMean = BigFloat.ScaleByPowerOfTwo(x1 + x2, -1)
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()
print "Iteration {0}: {1:F1} digits ({2:.3f} seconds)".format(k, correctionDigits, sw.Elapsed.TotalSeconds)
if correctionDigits >= digits:
break
pi = x1 * x1 / (1 - S)
sw.Stop()
print "Total time: {0:.3f} seconds.".format(sw.Elapsed.TotalSeconds, pi)
print
#
# 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
print "Quartic formula:"
sw.Reset()
sw.Start()
sqrt2 = BigFloat.Sqrt(2, goal, RoundingMode.TowardsNearest)
y = sqrt2 - BigFloat.One
a = BigFloat(6, goal) - BigFloat.ScaleByPowerOfTwo(sqrt2, 2)
y2 = y * y
y4 = y2 * y2
correctionDigits = 0
k=1
while 4 * correctionDigits < digits:
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()
print "Iteration {0}: {1:F1} digits ({2:.3f} seconds)".format(k, correctionDigits, sw.Elapsed.TotalSeconds)
k = k+1
pi = BigFloat.Inverse(a)
sw.Stop()
print "Total time: {0:.3f} seconds.".format(sw.Elapsed.TotalSeconds, pi)
print
#
# Method 5: The built-in method
#
# The method used to compute pi internally is an order of magnitude
# faster than any of the above.
print "Built-in function:"
sw.Reset()
sw.Start()
pi = BigFloat.GetPi(goal)
sw.Stop()
print "Total time: {0:.3f} seconds.".format(sw.Elapsed.TotalSeconds, pi)
# The highest precision value of pi is cached, so
# getting pi to any precision up to that is super fast.
print "Built-in function (cached):"
sw.Reset()
sw.Start()
pi = BigFloat.GetPi(goal)
sw.Stop()
print "Total time: {0:.3f} seconds.".format(sw.Elapsed.TotalSeconds, pi)
```