pi.py

"a linear RNN that receives ones as input and gives increasingly better approximations to pi as output"

import numpy as np
import math

def binary(digits: int):
    "Make a basis of powers of two of dimension `digits`, lowest bits first"
    return 1 << np.arange(digits)

def leibniz(n):
    "pi = 4*sum(leibniz(n) for n in range(infinity))"
    sign = (-1)**n
    deno = 2*n+1
    return sign/deno

T = 32768
D = int(math.log2(T)) # D: dimension of the hidden state depends on the maximum sequence length

if False:
    # leibniz approximation to pi, supervision for code below
    s = 0
    for i in range(0, T):
        s += leibniz(i)
        print(s*4)

# input layer: just ones
u = np.ones(T, dtype=np.int32)
# first linear recurrence: x_t = x_{t-1} + u_t, A = 1, B = 1
x = np.cumsum(u)
# first linear recurrence: C = 2, D = 1
y = 2 * x + 1
# grow
bits = x[:, np.newaxis] & binary(D) # move into basis of powers of two of dimension D
# shrink
x_mod_2 = bits @ np.eye(D)[0]  # x_mod_2 = x % 2, test the lowest bit

sign = -2 * x_mod_2 + 1
z = sign/y   # nonlinearity

# final recurrence: A = 1, B = 1
output = np.cumsum(z)
# C = 4, D = 4
pi = 4 + output * 4
print(pi)