12
\$\begingroup\$

I wrote this code to show that my reddit post is correct.

After the first two terms, the signs are determined as follows: If the denominator is a prime of the form 4m − 1, the sign is positive; if the denominator is a prime of the form 4m + 1, the sign is negative; for composite numbers, the sign is equal the product of the signs of its factors.

Basically it's the harmonic series minus the non-Gaussian prime reciprocals and reciprocals which factors are an odd multiple of the non-Gaussian primes- beautifully embodying quadratic reciprocity.

This is a very, very inefficient way to calculate π. However, I believe it's the most beautiful. Which is a hard sell, because π algorithms are pretty much the most harmonious things.

This formula for π is my favorite because it clearly shows how a circle is related to the harmonic series, and how that series is related to the prime number theorem and quadratic reciprocity. I love all algorithms for π, this is for its 'insight' it provides in relating π, primes, and quadratics.

Something to note is that swapping the logic here from n%4==1 to n%4==3 results in the sum = π/2 , and converges a bit faster.

import decimal
iters = int(input('Number of Iterations: '))
D = decimal.Decimal
decimal.getcontext().prec = 100
def prime_factors(n):
    i = 2
    factors = []
    while i * i <= n:
        if n % i:
            i += 1
        else:
            n //= i
            factors.append(i)
    if n > 1:
        factors.append(n)
    return factors
s = D(0)
for x in range(1, iters):
    clist = [int(i) for i in prime_factors(x)]
    plist = [n for n in clist if n%4==1]
    if len(plist)%2!=0:
        s-=1/D(x)
    else:
        s += 1/D(x)
    print(s)
\$\endgroup\$
  • \$\begingroup\$ Minor mistake, but it will hurt you sooner or later if you are not cautious: you ask the number of iterations, then loop on range(1, iters), which has length iters-1 (range(a, b) contains values from a to be excluding b. So, the loop should be on range(1, iters + 1) \$\endgroup\$ – user179210 Aug 2 at 5:42
  • \$\begingroup\$ Nice catch! Merci! \$\endgroup\$ – TheHoyt Aug 2 at 5:52
  • \$\begingroup\$ I could be wrong, but I thought the series for PI alternated signs with each term. \$\endgroup\$ – Thomas Matthews Aug 2 at 16:37
  • \$\begingroup\$ You're thinking of the Leibniz series my friend :) There are alot of series for pi. \$\endgroup\$ – TheHoyt Aug 2 at 16:59
  • 1
    \$\begingroup\$ @ThomasMatthews: if you alternate signs on each term, the sum is log 2. \$\endgroup\$ – Martin Argerami Aug 2 at 18:29
13
\$\begingroup\$

What you are approximating is $$ \pi =\sum_{m=1}^{\infty}\frac{(-1)^{s(m)}}{m},$$ where \$s(m)\$ counts the number of appearances of primes of the form \$4k+1\$ in the prime decomposition of \$m\$, compare

on Mathematics Stack Exchange.

The computation of the factor \$ (-1)^{s(m)} \$ can be done more efficiently: Instead of creating a list of all prime factors of \$m \$, then filtering the list for prime factors of the form \$4k+1\$, and finally counting the filtered list, you can compute the sign while factoring the number:

def sign_for_pi_series(n):
    """ Returns sign for 1/n in the pi series.

    The sign is -1 if n has an odd number of prime factors of the form 4k+1,
    and +1 otherwise.
    """
    s = 1
    i = 2
    while i * i <= n:
        if n % i:
            i += 1
        else:
            n //= i
            if i % 4 == 1:
                s = -s
    if n > 1 and n % 4 == 1:
        s = -s
    return s

The summation then simplifies to

sum_pi = sum(D(sign_for_pi_series(x))/D(x) for x in range(1, iters + 1))
print(sum_pi)
\$\endgroup\$
6
\$\begingroup\$

Blank lines

You require 2 blank lines before def statements to be PEP 8 compliant.

decimal.getcontext().prec = 100


def prime_factors(n):

Whitespace

You have your operators written as follows:

len(plist)%2!=0

You should use whitespace around operators instead:

 if len(plist) % 2 != 0:

Refactored

import decimal
iters = int(input('Number of Iterations: '))
D = decimal.Decimal
decimal.getcontext().prec = 100


def prime_factors(n):
    i = 2
    factors = []
    while i * i <= n:
        if n % i:
            i += 1
        else:
            n //= i
            factors.append(i)
    if n > 1:
        factors.append(n)
    return factors


s = D(0)
for x in range(1, iters):
    clist = [int(i) for i in prime_factors(x)]
    plist = [n for n in clist if n % 4 == 1]
    if len(plist) % 2 != 0:
        s -= 1 / D(x)
    else:
        s += 1 / D(x)
    print (s)
\$\endgroup\$
  • 3
    \$\begingroup\$ Sorry, your review is not useful at all. \$\endgroup\$ – stefan Aug 2 at 7:30
  • 15
    \$\begingroup\$ @stefan apology accepted \$\endgroup\$ – dfhwze Aug 2 at 7:41
  • 1
    \$\begingroup\$ I would also put a blank line between import and the rest of the code. Not sure if it is somewhere in PEP 8 or no, though. Also, you introduced a space on the last line which shouldn't be there. \$\endgroup\$ – Georgy Aug 2 at 14:45
  • 5
    \$\begingroup\$ @stefan - How is this not useful? There are code standards for a reason. In any case, +1. \$\endgroup\$ – Brian Aug 2 at 21:54
  • 1
    \$\begingroup\$ I almost thought you speak python but seeing the first comment... you probably thought the same :-P \$\endgroup\$ – t3chb0t Aug 23 at 6:42
5
\$\begingroup\$

You don't need plist, only its length. Likewise, you don't really need to build clist, its values would be enough, making prime_factors a generator. You may also remove the if. Note that given the extremely slow convergence, it's really not useful to compute with 100 decimals, so you could use floats instead.

import decimal

iters = int(input('Number of Iterations: '))
D = decimal.Decimal
decimal.getcontext().prec = 100


def prime_factors(n):
    i = 2
    while i * i <= n:
        if n % i:
            i += 1
        else:
            n //= i
            yield i
    if n > 1:
        yield n

s = D(0)
for n in range(1, iters + 1):
    k = sum(1 for p in prime_factors(n) if p % 4 == 1)
    s += (-1)**k / D(n)
    print(s)
\$\endgroup\$
  • 1
    \$\begingroup\$ Good point about the convergence, I shall simplify when I share this algorithm. I had been comparing 6 different Pi algorithms when writing this one and needed the precision there :) I shall share that project here too. \$\endgroup\$ – TheHoyt Aug 2 at 6:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.