Search code examples
pythonpiapproximation

Trying to define one of Euler's approximations to pi, getting unsupported operand type(s) for 'list and 'int'

I am trying to define a function which will approximate pi in python using one of Euler's methods. His formula is as follows: formula

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

Solution

  • 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:

    • You only want the odd primes, so you can start with a list of odd numbers.
    • 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.
    • In general it is very bad practise to remove elements from the list over which you are iterating. See e.g. this post. Internally, Python uses an index into the list over which it iterates. Coincidently, this is not a problem in this specific case, because only numbers larger than the current are removed.
    • Also, removing elements from a very long list is very slow, as each time the complete list needs to be moved to fill the gap. Therefore, it is better to work with two separate lists.
    • You didn't calculate the resulting product, nor did you return it.
    • As you notice, this formula converges very slowly.
    • As mentioned in the comments, the previous version interpreted 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