BigNumbers QuickStart Sample (IronPython)

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 the Extreme Optimization Mathematics Library for .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)