# Approximating constant π² to within error

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?

• Probably you want return prev (8 / (n * n) <= error < 8 / ((n - 1) * (n - 1)) ) – vaeta Mar 6 '18 at 1:56
• You probably already know this, but just in case: if you want a good-enough approximation of pi-squared rather than computing that particular formula, then hardcoding the result, exact to the limits of floating point precision, is better by almost any metric, and certainly has better performance. – Emilio M Bumachar Mar 6 '18 at 16:03
• You're going to run into rounding errors quite quickly with this one. – David Foerster Mar 7 '18 at 14:17

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

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

• this is by far the best way to do it. – Ev. Kounis Mar 6 '18 at 12:19
• @Ev.Kounis, no, it isn't. It addresses the first point raised in my (earlier, FWIW) answer, but not the second one, which is about numerical analysis and actually getting the correct answer. – Peter Taylor Mar 6 '18 at 12:36
• @PeterTaylor: The "correct answer" is to use a better formula because this one converges too slowly. With error=1e-14, 14 million terms are needed, sum(l) and sum(l[::-1]) are 6.5e-11 apart but they both are 1.5e-7 away from π². Somehow, the sum with the largest floats first is a better approximation than the other one. – Eric Duminil Mar 6 '18 at 13:01

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

• This solution seems less understandable at a glance than OP's. Does anyone agree or is it because I am not well versed/experienced with python? – Real Mar 6 '18 at 18:34
• Directly calling a "rich comparison" magic method such as __lt__() is a bad habit to get into. They are allowed to return NotImplemented, in which case you're expected to swap the operands and try the reverse magic method. operator.lt() will do the whole thing correctly in the face of arbitrary type mismatches (though in this case you'd also need functools.partial() to pass tolerance). – Kevin Mar 6 '18 at 19:32
• Just bear in mind that this is not a style issue - it's a correctness issue. So if you're doing something with NumPy or another of the "interesting" libraries, this difference could really bite you. – Kevin Mar 6 '18 at 19:37
• @Kevin I could also have used lambda diff: tolerance < diff which both corrected and readable. But, lambdas... Meh. – Mathias Ettinger Mar 6 '18 at 19:43
• @Real To me, the solution in this post is much more understandable, in the following sense: if I were given both pieces of code with no explanation, it's a lot easier for me to work out the "big picture" of what the answer's code is doing (summing a series), compared to doing the same for the code in the question. But to many other people, the answer's code will be less understandable. It makes sense that the amount of Python experience you have will influence that. – David Z Mar 7 '18 at 10:19
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):
...


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

• instead of break, just return new – hjpotter92 Mar 6 '18 at 1:24
• Point 4 is for the most part invalid. Just because most people use i, j, k for counting, doesn't mean that other variable names are any less valid. For example, if you're looping over $(x, y)$ pairs, x and y make more sense than i and j. Likewise, if some argument is called $n$ in a mathematical operation, it makes sense to use n to help people familiar with the maths understand your code. – Daniel Mar 6 '18 at 7:16
• From the math perspective, n is most often the maximum value of the counter and the current value is more likely to be k. – allo Mar 6 '18 at 10:47
• You calculate new = prev + something and then new - prev? – Eric Duminil Mar 6 '18 at 11:10
• @SolomonUcko Speaking of unnecessary parentheses, the new = (prev + ... has quite a few. You can write new = prev + 8 / i**2 – Ev. Kounis Mar 6 '18 at 12:17

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

• Good catch! For smaller numbers of decimals (<10) it seems to be a 50-100% speedup. – maxb Mar 6 '18 at 9:37
• It really isn't too hard to calculate the exact n, your method should be called approx_pi_squared_fast BTW ;) Other than that, it looks good and it's much faster than the alternatives. – Eric Duminil Mar 6 '18 at 11:51
• It is possible to find a more accurate value for n, but I mainly wanted to be in the same ballpark, and have it be larger than the number of iterations in the original question for any number of decimals (including floating point). – maxb Mar 6 '18 at 11:57

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.

• Reversing the terms doesn't seem to have much of an influence in this case. Do you have an example where it leads to a significant difference? – Eric Duminil Mar 6 '18 at 16:36
• @EricDuminil : it will make a difference if you set error small enough. Specifically if it is less than about 1E-15, adding successively smaller terms will make no difference to the output value. – Martin Bonner Mar 7 '18 at 14:46
• @MartinBonner: With error=1e-15, the formula will still be off by 4.5e-8. It might be an important point, I just don't think it applies to this case. – Eric Duminil Mar 7 '18 at 15:02
• @EricDuminil, despite the title of the question, the actual spec as presented asks for a specific sum, and it's the error in calculating that sum which I consider relevant, not the difference relative to $\pi^2$. – Peter Taylor Mar 7 '18 at 17:22
• Please don't use the name sum. It's a python builtin. – Snakes and Coffee Mar 8 '18 at 6:42

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.

• Question after your suggestion is, what would be the value of diff in the first iteration then – A. Romeu Mar 6 '18 at 12:37
• @A.Romeu whatever the hell you want it to be, as long as it's bigger than error. error + 1 would work. – Nico Mar 6 '18 at 13:51

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.

• Are you aware of the recursion limit in Python? You won't be able to go much below than error = 0.00001 – Mathias Ettinger Mar 6 '18 at 13:41
• Yes, I'm aware of it. I was giving my answer as another way to think about it. I'm not sure of OP's intent or if they even need to go lower. – TheStrangeQuark Mar 6 '18 at 13:46
• Recursion could be interesting for series which converge very fast. It doesn't make much sense for OPs case but it could work fine for Ramanujan's formula for example. – Eric Duminil Mar 6 '18 at 14:34
• Yeah, I agree. I didn't investigate how quickly OP's formula converges, or even that it does converge to something meaningful. – TheStrangeQuark Mar 6 '18 at 15:57

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

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