I am trying to define a function which will approximate pi
in python using one of Euler's methods. His formula is as follows:
My code so far is this:
def pi_euler1(n):
numerator = list(range(2 , n))
for i in numerator:
j = 2
while i * j <= numerator[-1]:
if i * j in numerator:
numerator.remove(i * j)
j += 1
for k in numerator:
if (k + 1) % 4 == 0:
denominator = k + 1
else:
denominator = k - 1
#Because all primes are odd, both numbers inbetween them are divisible by 2,
#and by extension 1 of the 2 numbers is divisible by 4
term = numerator / denominator
I know this is wrong, and also incomplete. I'm just not quite sure what the TypeError that I mentioned earlier actually means. I'm just quite stuck with it, I want to create a list of the terms and then find their products. Am I on the right lines?
Update: I have worked ways around this, fixing the clearly obvious errors that were prevalent thanks to msconi and Johanc, now with the following code:
import math
def pi_euler1(n):
numerator = list(range(2 , 13 + math.ceil(n*(math.log(n)+math.log(math.log(n))))))
denominator=[]
for i in numerator:
j = 2
while i * j <= numerator[-1]:
if (i * j) in numerator:
numerator.remove(i * j)
j += 1
numerator.remove(2)
for k in numerator:
if (k + 1) % 4 == 0:
denominator.append(k+1)
else:
denominator.append(k-1)
a=1
for i in range(n):
a *= numerator[i] / denominator[i]
return 4*a
This seems to work, when I tried to plot a graph of the errors from pi in a semilogy axes scale, I was getting a domain error, but i needed to change the upper bound of the range to n+1 because log(0) is undefined. Thank you guys
Here is the code with some small modifications to get it working:
import math
def pi_euler1(n):
lim = n * n + 4
numerator = list(range(3, lim, 2))
for i in numerator:
j = 3
while i * j <= numerator[-1]:
if i * j in numerator:
numerator.remove(i * j)
j += 2
euler_product = 1
for k in numerator[:n]:
if (k + 1) % 4 == 0:
denominator = k + 1
else:
denominator = k - 1
factor = k / denominator
euler_product *= factor
return euler_product * 4
print(pi_euler1(3))
print(pi_euler1(10000))
print(math.pi)
Output:
3.28125
3.148427801913721
3.141592653589793
Remarks:
j
can start with 3
and increment in steps of 2. In fact, j
can start at i
because all the multiples of i
smaller than i*i
are already removed earlier.n
as the limit for highest prime, while in fact n
should be the number of primes. I adapted the code to rectify that. In the above version with a crude limit; the version below tries a tighter approximation for the limit.Here is a reworked version, without removing from the list you're iterating. Instead of removing elements, it just marks them. This is much faster, so a larger n
can be used in a reasonable time:
import math
def pi_euler_v3(n):
if n < 3:
lim = 6
else:
lim = n*n
while lim / math.log(lim) / 2 > n:
lim //= 2
print(n, lim)
numerator = list(range(3, lim, 2))
odd_primes = []
for i in numerator:
if i is not None:
odd_primes.append(i)
if len(odd_primes) >= n:
break
j = i
while i * j < lim:
numerator[(i*j-3) // 2] = None
j += 2
if len(odd_primes) != n:
print(f"Wrong limit calculation, only {len(odd_primes)} primes instead of {n}")
euler_product = 1
for k in odd_primes:
denominator = k + 1 if k % 4 == 3 else k - 1
euler_product *= k / denominator
return euler_product * 4
print(pi_euler_v2(100000))
print(math.pi)
Output:
3.141752253548891
3.141592653589793