Python π = 1 + (1/2) + (1/3) + (1/4) - (1/5) + (1/6) + (1/7) + (1/8) + (1/9) - (1/10) ...1748 Euler 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 &lt;= n:
        if n % i:
            i += 1
        else:
            n //= i
            factors.append(i)
    if n &gt; 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)
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 


How can we prove \$\pi =1+\frac{1}{2}+\frac{1}{3}+\frac{1}{4}-\frac{1}{5}+\frac{1}{6}+\frac{1}{7}+\cdots\,\$?


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 &lt;= n:
        if n % i:
            i += 1
        else:
            n //= i
            if i % 4 == 1:
                s = -s
    if n &gt; 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)