#!/usr/bin/env lua

--[[
This program computes Pi by the formula:
pi / 4 = 4 * ArcTan(1/5) - ArcTan(1/239)

Adapted from some program I found on the Internet a long time
ago that computed Pi by:
pi / 4 = ArcTan(1/2) + ArcTan(1/3)

I have no idea who wrote the original program nor remember where
I found it--there wasn't any credits in the original C source.  This
one also no longer looks like the original as I've completely rewritten
it.  I got a lot of my information from
http://www.boo.net/~jasonp/pipage.html when I rewrote the program.  This
program will reliably compute upto 1,000,000 digits, that I've verified.

There are faster ways of computing Pi, and this program probably
could be optimized further, but it is sufficiently fast for computing
Pi using relatively portable and understandable routines that can
be easily converted into other languages.

Computing Pi using FFTs are much much faster, but also much more
complicated.  Using the ArcTan methods are fairly simple to implement,
but not nearly as fast as the FFT methods.  Also, the
pi / 4 = 4 * ArcTan(1/5) - ArcTan(1/239) formula is one of the faster
of the various ArcTan formulas for computing PI.

Author: M. Scott Reynolds
Date: Lua adaptation 26 October 2020
]]

--[[
 * Compute pi up to maxDigits long.
 * Return result as a string. 
]]

function debug(...)
--	print("DEBUG:", ...)
end

function write(...)
	io.write(...)
end

function writeln(...)
	io.write(...)
	io.write("\n")
end

function computePi(maxDigits)
    local SIZE = 1000
	
    local precision = (maxDigits-1) // 3 + 3
	debug(precision)
	
	-- initialize variables
	local remainder1, remainder2, remainder3, remainder4 = 0, 0, 0, 0
	local b, n, n2, carry = 0, 0, 0, 0
	local i, l = 0, 0
	local isZero = false

	-- arrays of numbers
    local p = {}
	local t = {}

    -- Initialize arrays.
    for i = 1, precision+1 do
        t[i] = 0
		p[i] = 0
    end
	
    -- Note, borrows and carries from the addition and subtraction
    -- are postponed till last.  See: http://www.boo.net/~jasonp/ord

    -- Compute arctan(1/5)

    -- t = t / 5, p = t
    t[1] = 1
    remainder1 = 0
    for i = 1, precision+1 do
        b = SIZE * remainder1 + t[i]
        t[i] = b // 5
        p[i] = t[i]
        remainder1 = b % 5
    end

    -- While t is not zero.
    n = -1
    n2 = 1
	isZero = false
    repeat
        -- t = t / 25, p = p - t / n, t = t / 25, p = p + t / (n+2)
        remainder1 = 0
        remainder2 = 0
        remainder3 = 0
        remainder4 = 0
        isZero = true
        n = n + 4
        n2 = n2 + 4
        for i = 1, precision+1 do
            b = SIZE * remainder1 + t[i]
            t[i] = b // 25
            remainder1 = b % 25

            b = SIZE * remainder2 + t[i]
            p[i] = p[i] - b // n
            remainder2 = b % n

            b = SIZE * remainder3 + t[i]
            t[i] = b // 25
            remainder3 = b % 25

            b = SIZE * remainder4 + t[i]
            p[i] = p[i] + b // n2
            remainder4 = b % n2 

            if isZero and t[i] ~= 0 then
                isZero = false
            end
        end
    until isZero

    -- p = p * 4
    carry = 0
    for i = precision+1, 1, -1 do
        b = p[i] * 4 + carry
        p[i] = b % SIZE
        carry = b // SIZE
    end

    -- Compute arctan(1/239)

    -- t = t / 239, p = p - t
    t[1] = 1
    remainder1 = 0
    for i = 1, precision+1 do
        b = SIZE * remainder1 + t[i]
        t[i] = b // 239
        p[i] = p[i] - t[i]
        remainder1 = b % 239
    end

    -- While t is not zero.
    n = -1
    n2 = 1
    repeat
        -- t = t / 57121, p = p + t / n, t = t / 57121, p = p - t / (n+2)
        remainder1 = 0
        remainder2 = 0
        remainder3 = 0
        remainder4 = 0
        isZero = true
        n = n + 4
        n2 = n2 + 4
        for i = 1, precision+1 do
            b = SIZE * remainder1 + t[i]
            t[i] = b // 57121
            remainder1 = b % 57121

            b = SIZE * remainder2 + t[i]
            p[i] = p[i] + b // n
            remainder2 = b % n

            b = SIZE * remainder3 + t[i]
            t[i] = b // 57121
            remainder3 = b % 57121

            b = SIZE * remainder4 + t[i]
            p[i] = p[i] - b // n2
            remainder4 = b % n2

            if isZero and t[i] ~= 0 then    -- is t zero?
                isZero = false
            end
        end
    until isZero

    -- p = p * 4
    carry = 0
    for i = precision+1, 1, -1 do
        b = p[i] * 4 + carry
        p[i] = b % SIZE
        carry = b // SIZE
    end

    -- Borrow and carry.
    for i = precision+1, 2, -1 do
        if p[i] < 0 then
            b = p[i] // SIZE
            p[i] = p[i] - (b - 1) * SIZE
            p[i-1] = p[i-1] + b - 1
        end
        if p[i] >= SIZE then
            b = p[i] // SIZE
            p[i] = p[i] - b * SIZE
            p[i-1] = p[i-1] + b
        end
    end

    ---[[ Store results in string buffer.
	local text = {}
	text[1] = tostring(p[1])
	text[2] = '.'
    for i = 2, precision+1 do
        if p[i] < 10 then
			text[#text + 1] = "00" .. tostring(p[i])
        elseif p[i] < 100 then
            text[#text + 1] = '0' .. tostring(p[i])
        else
			text[#text + 1] = tostring(p[i])
        end
    end
	--]]
    return table.concat(text):sub(1, maxDigits+2)
end

do
	local digits = math.tointeger(arg[1]) or 10000
	writeln("Computing ", digits, " digits of Pi.")
	
	local pi = computePi(digits)
	
	-- Format and print the output
	local line = ""
	writeln("3.")
	for i = 3, pi:len() do
		local j = i - 2
		line = line .. pi:sub(i, i)
		if j % 1000 == 0 then
			writeln(line)
			writeln()
			line = ""
		elseif j % 50 == 0 then
			writeln(line)
			line = ""
		elseif j % 10 == 0 then
			line = line .. " "
		end
    end
	if line:len() > 0 then
		writeln(line)
	end
end