### One Million Binary Digits of π

A file containing the first 2^{20}=1,048,576 binary digits of π after the radix point can be found at `http://weitz.de/files/pi.zip`

. The file is meant to be read byte by byte (i.e. octet by octet) and each individual byte should be read from least to most significant bit. Here's some demo code in Python to read the file bit by bit:

def bits (name):
with open(name, "rb") as infile:
byte = infile.read(1)
while byte:
num = int.from_bytes(byte, "little")
for i in range(8):
yield num % 2
num >>= 1
byte = infile.read(1)

The file was created using the Julia program below. It implements a so-called "streaming" spigot algorithm by Jeremy Gibbons. The code can in theory create an unlimited number of digits although it will become slower over time and its memory consumption will increase. (Creating one million binary digits took about two hours on my run-of-the-mill PC.)

using Printf, Dates # for Julia 1.x
function spigot(n)
V1::BigInt = 1
V2::BigInt = 0
V4::BigInt = 1
num, den, prod = 1, 3, 3
i, j, byte = -1, 0, 0
start = now()
open("pi", "w") do out
while i <= n
if 4V1+V2-V4 < prod*V4
if i >= 0
byte += prod << j
j += 1
if j == 8
write(out, byte % UInt8)
j, byte = 0, 0
end
end
i += 1
if i % 1000 == 0
@printf("%9d, %s\n", i, Dates.canonicalize(Dates.CompoundPeriod(now() - start)))
end
# see remark below
V1, V2, prod = 2V1, 2*(V2-prod*V4), div(2*(3V1+V2),V4)-2prod
else
V1, V2, V4, prod = num*V1, den*(2V1+V2), den*V4, div(7num*V1+2+den*V2,den*V4)
num += 1
den += 2
end
if num % 50 == 0
g::BigInt = gcd(V1, V2, V4)
V1, V2, V4 = div(V1, g), div(V2, g), div(V4, g)
end
end
end
end

If you want to emit decimal instead of binary digits, change the four 2's in the line below the comment to 10's. Of course, you'll then have to change the output code as well which currently collects bits until it can write one byte to the file. (That's the `if i >= 0`

part which is executed for each digit with the digit in the variable `prod`

.)