BigNumbers in IronPython QuickStart Sample

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

This sample is also available in: C#, Visual Basic, F#.

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

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 Extreme 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)