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)