Search code examples
pythonmathpi

Why is the computing of the value of pi using the Machin Formula giving a wrong value?


For my school project I was trying to compute the value of using different methods. One of the formula I found was the Machin Formula that can be calculated using the Taylor expansion of arctan(x).

I wrote the following code in python:

import decimal

count = pi = a = b = c = d = val1 = val2 = decimal.Decimal(0) #Initializing the variables      
decimal.getcontext().prec = 25 #Setting percision

while (decimal.Decimal(count) <= decimal.Decimal(100)): 
    a = pow(decimal.Decimal(-1), decimal.Decimal(count))
    b = ((decimal.Decimal(2) * decimal.Decimal(count)) + decimal.Decimal(1))
    c = pow(decimal.Decimal(1/5), decimal.Decimal(b))
    d = (decimal.Decimal(a) / decimal.Decimal(b)) * decimal.Decimal(c)
    val1 = decimal.Decimal(val1) + decimal.Decimal(d)
    count = decimal.Decimal(count) + decimal.Decimal(1)
    #The series has been divided into multiple small parts to reduce confusion

count = a = b = c = d = decimal.Decimal(0) #Resetting the variables

while (decimal.Decimal(count) <= decimal.Decimal(10)):
    a = pow(decimal.Decimal(-1), decimal.Decimal(count))
    b = ((decimal.Decimal(2) * decimal.Decimal(count)) + decimal.Decimal(1))
    c = pow(decimal.Decimal(1/239), decimal.Decimal(b))
    d = (decimal.Decimal(a) / decimal.Decimal(b)) * decimal.Decimal(c)
    val2 = decimal.Decimal(val2) + decimal.Decimal(d)
    count = decimal.Decimal(count) + decimal.Decimal(1)
    #The series has been divided into multiple small parts to reduce confusion

pi = (decimal.Decimal(16) * decimal.Decimal(val1)) - (decimal.Decimal(4) * decimal.Decimal(val2))
print(pi)

The problem is that I am getting the right value of pi only till 15 decimal places, no matter the number of times the loop repeats itself.

For example:

at 11 repetitions of the first loop

pi = 3.141592653589793408632493

at 100 repetitions of the first loop

pi = 3.141592653589793410703296

I am not increasing the repetitions of the second loop as arctan(1/239) is very small and reaches an extremely small value with a few repetitions and therefore should not affect the value of pi at only 15 decimal places.

EXTRA INFORMATION:

The Machin Formula states that:

   π = (16 * Summation of (((-1)^n) / 2n+1) * ((1/5)^(2n+1))) - (4 * Summation of (((-1)^n) / 2n+1) * ((1/239)^(2n+1)))    

Solution

  • That many terms is enough to get you over 50 decimal places. The problem is that you are mixing Python floats with Decimals, so your calculations are polluted with the errors in those floats, which are only precise to 53 bits (around 15 decimal digits).

    You can fix that by changing

    c = pow(decimal.Decimal(1/5), decimal.Decimal(b))
    

    to

    c = pow(1 / decimal.Decimal(5), decimal.Decimal(b))
    

    or

    c = pow(decimal.Decimal(5), decimal.Decimal(-b))
    

    Obviously, a similar change needs to be made to

    c = pow(decimal.Decimal(1/239), decimal.Decimal(b))
    

    You could make your code a lot more readable. For starters, you should put the stuff that calculates the arctan series into a function, rather than duplicating it for arctan(1/5) and arctan(1/239).

    Also, you don't need to use Decimal for everything. You can just use simple Python integers for things like count and a. Eg, your calculation for a can be written as

    a = (-1) ** count
    

    or you could just set a to 1 outside the loop and negate it each time through the loop.

    Here's a more compact version of your code.

    import decimal
    
    decimal.getcontext().prec = 60 #Setting precision
    
    def arccot(n, terms):
        base = 1 / decimal.Decimal(n)
        result = 0
        sign = 1
        for b in range(1, 2*terms, 2):
            result += sign * (base ** b) / b
            sign = -sign
        return result
    
    pi = 16 * arccot(5, 50) - 4 * arccot(239, 11)
    print(pi)
    

    output

    3.14159265358979323846264338327950288419716939937510582094048
    

    The last 4 digits are rubbish, but the rest are fine.