Approximating constant π² to within error

Question

This function takes as input float error and approximates constant π² to within error by computing this sum, term by term, until the difference between the new and the previous sum is less than error. The function should return the new sum.

$$ \pi^2 = 8 + \dfrac{8}{3^2} + \dfrac{8}{5^2} + \dfrac{8}{7^2} + \dfrac{8}{9^2} + \cdots $$

Example:

approxPIsquared(0.0001)

Result is 9.855519952254232

My working solution is shown below:

def approxPIsquared(error):
    prev = 8
    new =0
    n = 3
    while (True):
        new = (prev + (8 / (n * n)))
        diff = new - prev
        if (diff <= error):
            break
        prev = new
        n = n + 2

    return new

Is this a good solution or is there a better way to do this?

Answer 1

You don't need to compare prev and new during each iteration.

The difference between the new and the previous sum is simply the current term : $$\frac{8}{(2i+1)^2}$$

If you want this term to be smaller than error, you can solve:

$$\mathrm{error} > \frac{8}{(2i+1)^2}\\ \iff (2i+1)^2 > \frac{8}{error}\\ \iff 2i+1 > \sqrt{\frac{8}{error}}\\ \iff i > \frac{\sqrt{\frac{8}{error}} - 1}{2}\\ $$

Now that you know how many terms your series should have, you can return the result directly:

def approx_pi_squared(error):
    n = int(((8 / error)**0.5 - 1) / 2) + 1
    return sum(8 / (2 * i + 1)**2 for i in range(n))

Better formulas

Adding a delta

Note that error represents how small the terms are, not how close approx_pi_squared is from π²:

>>> import math
>>> approx_pi_squared(1e-10)
9.869590258918535
>>> math.pi**2 - approx_pi_squared(1e-7)
0.0004472271895057389
>>> math.pi**2 - approx_pi_squared(1e-10)
1.414217082285063e-05

Even with more than 140 000 terms, the series only gives the 3 first digits of π². This formula is very simple but converges too slowly.

What's very interesting, though, is that the difference between math.pi**2 and approx_pi_squared(error) seems very close to \$\sqrt{2\mathrm{error}}\$. It seems to hold true for any error, so we could update the function:

def approx_pi_squared(error):
    n = int(((8 / error)**0.5 - 1) / 2) + 1
    delta = (2 * error)**0.5
    return sum(8 / (2 * i + 1)**2 for i in range(n)) + delta

approx_pi_squared(1e-10) now returns 10 correct digits for π².

This new formula is unproven, so use at your own risk!

BBP-Type Formula

There are many π² formulas, so feel free to pick another one. For example:

def approx_pi_squared(error):
    n = int(((12 / error)**0.5 - 1) / 2) + 1
    return 12 * sum((-1)**i / (i + 1)**2 for i in range(n))

error seems to have the same order of magnitude as math.pi**2 - approx_pi_squared(error) now:

>>> math.pi**2 - approx_pi_squared(1e-9)
2.0001476030984122e-09
>>> math.pi**2 - approx_pi_squared(1e-10)
-1.9977974829998857e-10

delta looks like (-1)**n * 2 * error.

With Sympy

You can delegate the job to sympy and be sure you'll get a correct result with arbitrary precision:

>>> from sympy import Sum, N, pi, oo, init_printing
>>> from sympy.abc import i, n
>>> init_printing()
>>> Sum(8/(2*i+1)**2, (i, 0, oo))
  ∞             
 ____           
 ╲              
  ╲       8     
   ╲  ──────────
   ╱           2
  ╱   (2⋅i + 1) 
 ╱              
 ‾‾‾‾           
i = 0  
>>> N(Sum(8/(2*i+1)**2, (i, 0, oo)), 100)
9.869604401089358618834490999876151135313699407240790626413349376220044822419205243001773403718552232
>>> N(pi**2, 100)
9.869604401089358618834490999876151135313699407240790626413349376220044822419205243001773403718552232

Answer 2

So, you have a polynomial sequence and want to sum its terms while they are greater than a provided tolerance (i.e.: the tolerance is lower than the current computed term).

This is easily expressed using basic python constructs / builtins and the itertools module:

1. An infinite polynomial sequence can be described using a generator expression and the infinite counter itertools.count:

   polynomial_sequence = (8 / (n * n) for n in itertools.count(1, 2))
   

2. You can extract terms out of an iterable (and generator expressions are iterables) while they comply to a condition using itertools.takewhile:

   approximate_finite_sequence = itertools.takewhile(tolerance.__lt__, polynomial_sequence)
   

   Here __lt__ is the magic method that is called when tolerance < … is written, so the terms of the sequence will be kept while tolerance < term. But contrary to your implementation, the first term that is lower than the tolerance will not be kept and thus not added to the sum for the approximation.

3. You can sum all the generated terms using the builtin sum:

   pi_squared = sum(approximate_finite_sequence)
   

---

Putting it all togeter:

import itertools


def approximate_pi_squared(tolerance=0.0001):
    polynomial_sequence = (8 / (n * n) for n in itertools.count(1, 2))
    return sum(itertools.takewhile(tolerance.__lt__, polynomial_sequence))

Answer 3

1. Please fix the indentation (probably just a copy-pasting issue).
2. In Python, you don't need parenthesis around the expressions for block statements like if and while. You also don't need it around your expressions (thanks Ev. Kounis).
3. You can use i ** 2 instead of i * i.
4. You can use the increment operator: n += 2 instead of n = n + 2.
5. You should probably use i instead of n for counter variables (in most cases [thanks Coal_]).

6. You can use (make sure to import itertools; I personally prefer to use from itertools import count):

   for n in itertools.count(3, 2):
       ...
   

   instead of

   n = 3
   while (True):
       ...
       n = n + 2
   

7. There's no need to have a temporary diff variable.
8. Use snake_case for everything but classes and 'constants'.
9. You can return instead of breaking (thanks @hjpotter92).

Result:

from itertools import count

def approx_pi_squared(error):
    prev = 8
    new = 0
    for i in count(3, 2):
        new = prev + 8 / i ** 2
        if new - prev <= error:
            return new
        prev = new

Answer 4

Since Solomon Ucko provided excellent feedback on the coding style, and gave an appropriate modification given your original algorithm, I thought I'd focus on the performance tag you added.

There are many more efficient formulas for calculating pi, but let's use the one you provided, and make it faster. By running your function for different inputs on the form 10^(-n), I found out that your program performed about 10**(n/2+.5) iterations (it was about 5-10% less).

Once I had bounded the number of iterations, I could resort to using numpy, which is blazingly fast for these kinds of operations. This is the script I ended up using:

import numpy as np

def approx_pi_squared_fast(decimals):
    n = int(10**(decimals/2+.5))
    denominators = np.arange(1, 2*n, 2)**2
    pi_squared = np.sum(8/denominators)
    return pi_squared

I changed the input from error to decimals, so the new program doesn't return the exact same values as the ones you posted. However, it does return more exact values for all inputs between 1-15 (after that your script takes >10s to test). It also returns the answer 4-6 times faster than your original script.

Edit: edited function name

Answer 5

Ockham's razor (do not unnecessarily multiply entities) can be applied to variables as well.

    prev = 8
    new = 0
    ...
        new = (prev + (8 / (n * n)))
        diff = new - prev
        ...
        prev = new

is unnecessarily complicated. You know diff, so you can reduce prev and new to one variable, with a more informative name (e.g. sum):

    sum = 8
    ...
        diff = 8 / (n * n)
        sum += diff
        ...

---

Having made it explicit what the difference is, we can address the biggest problem: correctness. To correctly sum a list of floating point numbers, you should start with the smallest ones, not the largest ones. But since we now have a simple expression for diff as a function of n, we can invert it to find the first value of n for which the term is less than the desired error, using a sqrt.

Answer 6

All of the existing answers provide valuable feedback, but there's one point I think is lacking, a discussion of the loop operator choice. In your initial algorithm you used while True: along with a if... break. Why not put the condition you evaluate in your if in your while, i.e. something like while diff > error:?

This has the advantage of making immediately clear what the stopping condition is, hence making your code easier to understand. This is the one default in Solomon Ucko's otherwise great answer, at least in my opinion. His for n in itertools.count(3, 2): does not say anything about termination condition.

Answer 7

Another way to rewrite this problem is using recursion. You can supply the function back into the function if it does not fit the error bounds by doing something like this:

def approx_pi_squared(error, curr=8, i=3):
    new = curr + (8 / (i * i))
    if new - curr <= error:
        return new
    else:
        return approx_pi_squared(error, new, i+2)

So if the current iteration does not fit the error bounds, then you will just do the function again, but use the previous value and i+2 in the function. In my opinion, it looks a little cleaner even though it doesn't save any time.

Since you're only approximating up to a certain decimal, you may also want to truncate your value only up the decimal place that your error is bounded to, since everything after is not certain. You can do that using the round() function like this:

def approx_pi_squared(error, curr=8, i=3):
    new = curr + (8 / (i * i))
    if new - curr <= error:
        decimal_index = str(error).index('.') + 1  # after decimal
        num_decimal = len(str(error)[decimal_index:])
        return round(new, num_decimal)
    else:
        return approx_pi_squared(error, new, i+2)

All this is doing is finding where the decimal is in the error input and counts the number of digits after that so the output matches the number of decimals as the error input. So if you run this you get approx_pi_squared(0.0001) you are returned 9.8555. Looking up pi * pi, I find that it is closer to 9.8696. Where did your formula come from? I'm trying to find it, but am not able to.

Answer 8

When you add lots of floating point numbers, have a look at Kahan summation, which looks complicated first, but it reduces the rounding error.

Answer 9

By the integral test, the distance of the partial sum of the first \$N\$ terms to the limit is about $$\sum_{n=N}^\infty\frac{8}{(2n+1)^2}\simeq\int_{N-0.5}^\infty\frac2{(x+0.5)^2}dx=\frac2{N}.$$ This implies the often reported fact that to get an accuracy of about 1e-14 you need about 1e14 terms in the partial sum.

The error of this error approximation is \begin{align}\sum_{n=N}^\infty\frac{8}{(2n+1)^2}-\frac2N &= \sum_{n=N}^\infty\frac{8}{(2n+1)^2}-\left(\frac2n-\frac2{n+1}\right) \\ &= -\sum_{n=N}^\infty\frac{1}{2n(n+0.5)^2(n+1)}\simeq-\frac1{6N^3}\end{align}

Now use the first error estimate to correct the partial sum to \$\sum_{i=0}^{N-1}\frac{8}{(2i+1)^2}+\frac2N\$. Then the error estimate of this formula is the second error estimate above. Thus to get an error of \$ε\$ with the corrected formula one needs \$N=1/\sqrt[3]{6ε}\$ terms. For example, to get \$ε=10^{-14}\$ one would have to sum up about \$N=2.6\cdot 10^4\$ terms.

Implementing this using summation from smallest to largest to mitigate floating point noise accumulation

def approx_pi_squared(error):
    n = int( 1 / (6*error)**(1/3) ) + 1
    return sum(8 / (2 * (n - i) - 1)**2 for i in range(n)) + 2 / n

for k in range(4,17):
     pi2 = approx_pi_squared(10**-k);
     print "%2d   %20.15f   %10.5e"%(k, pi2, pi2-pi**2)

gives the expected results

k       approx pi squared     error

 4      9.869700618552868   9.62175e-05
 5      9.869613878813270   9.47772e-06
 6      9.869605350023797   9.48934e-07
 7      9.869604499989549   9.89002e-08
 8      9.869604411023413   9.93406e-09
 9      9.869604402085667   9.96309e-10
10      9.869604401189266   9.99076e-11
11      9.869604401099350   9.99201e-12
12      9.869604401090358   1.00009e-12
13      9.869604401089459   1.01252e-13
14      9.869604401089369   1.06581e-14
15      9.869604401089360   1.77636e-15
16      9.869604401089360   1.77636e-15