Why is the computing of the value of pi using the Machin Formula giving a wrong value?
<p>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).    </p>

<p>I wrote the following code in python:</p>

<pre><code>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) &lt;= 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) &lt;= 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)
</code></pre>

<p>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.</p>

<p>For example:</p>

<p>at 11 repetitions of the first loop</p>

<p>pi = 3.141592653589793408632493</p>

<p>at 100 repetitions of the first loop</p>

<p>pi = 3.141592653589793410703296</p>

<p>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.</p>

<p>EXTRA INFORMATION:</p>

<p>The Machin Formula states that:</p>

<pre><code>   Ï = (16 * Summation of (((-1)^n) / 2n+1) * ((1/5)^(2n+1))) - (4 * Summation of (((-1)^n) / 2n+1) * ((1/239)^(2n+1)))    
</code></pre>

<p>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).</p>

<p>You can fix that by changing </p>

<pre><code>c = pow(decimal.Decimal(1/5), decimal.Decimal(b))
</code></pre>

<p>to </p>

<pre><code>c = pow(1 / decimal.Decimal(5), decimal.Decimal(b))
</code></pre>

<p>or</p>

<pre><code>c = pow(decimal.Decimal(5), decimal.Decimal(-b))
</code></pre>

<p>Obviously, a similar change needs to be made to </p>

<pre><code>c = pow(decimal.Decimal(1/239), decimal.Decimal(b))
</code></pre>

<hr>

<p>You could make your code a <em>lot</em> 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).</p>

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

<pre><code>a = (-1) ** count
</code></pre>

<p>or you could just set <code>a</code> to 1 outside the loop and negate it each time through the loop.</p>

<p>Here's a more compact version of your code.</p>

<pre><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)
</code></pre>

<p><strong>output</strong></p>

<pre><code>3.14159265358979323846264338327950288419716939937510582094048
</code></pre>

<p>The last 4 digits are rubbish, but the rest are fine.</p>

python math pi