One Million Binary Digits of π

A file containing the first 2²⁰=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.)