I have the following datastructure in Haskell for computable numbers $\mathbb{M}$:
data Computable = Exact Rational | Inexact (Word -> Integer)
Exact labels numbers in $\mathbb{Q}$. Inexact labels Cauchy sequences whose limit is in $\mathbb{M}$, like this (using lambda calculus):
$$ \text{Inexact}^{-1}(x) = \lambda k. \lfloor x × 10^k \rceil $$
Note that numbers in $\mathbb{Q}$ can be labeled Inexact. In other words, numbers labeled Exact has guarantee to be in $\mathbb{Q}$, while those labeled Inexact doesn't.
I want to declare instance Floating Computable, which requires a definition of $\pi$. I decided to implement it using Gauss-Legendre algorithm, like this:
approximateQ :: Rational -> Word -> Integer
approximateQ q n = round (q * 10 ^ n)
inexactToReal :: (Word -> Computable) -> Computable
inexactToReal f = Inexact $ \n -> case f n of
Exact x -> approximateQ x n
Inexact g -> g n
instance Floating Computable where
pi = inexactToReal $ \n -> go n 1 (sqrt 0.5) 0.25 1 where
go n a b t x = if a - b < 10 ^^ negate n
then (a + b)^2 / (4 * t)
else let
newA = (a + b) / 2
newB = sqrt (b * a)
in go n newA newB (t - x * (newA - a)^2) (2 * x)
$\text{approximateQ}(q,n) = \lfloor q × 10^n \rceil$ where $q \in \mathbb{Q}$ and $n \in \mathbb{N}$. inexactToReal flattens nested Inexacts.
But there is a problem. The implementation I refered to is in this paper, but it looks like the paper assumes the comparison to end iterations (a - b < 10 ^^ negate n) will be done in $O(1)$ time. But my implementation won't:
inexactify :: Computable -> Computable
inexactify (Exact x) = Inexact (approximateQ x)
inexactify x = x
instance Ord Computable where
compare (Exact x) (Exact y) = compare x y
compare (Inexact f) (Inexact g) = go 0 where
go n = compare (f n) (g n) <> go (n + 1)
compare x y = compare (inexactify x) (inexactify y)
inexactify removes the guarantee that the number is in $\mathbb{Q}$. Let $n$ be the location below the decimal point where the numbers start to differ. Since it takes $O(m)$ time to compare $m$-digit integers, assuming it takes $O(1)$ to evaluate both operands, it will take $O(n^2)$ time to compare numbers in $\mathbb{M}$. If the numbers are equal, the method won't halt. This is inevitable, for it is proven that there doesn't exist an algorithm that compares numbers in $\mathbb{M}$.
Furthermore, it actually doesn't take $O(1)$ to evaluate both operands. Gauss-Legendre algorithm recursively stacks sqrts. Given that there is $m$ sqrts stacked and I want to evaluate it, it takes $\Omega(2^m)$ time. If I were to evaluate $\pi$ to $n$ precisions, though I didn't do exact analysis, I persume it will take $\Omega(n!)$ time.
TL;DR: So I have a different plan. Given that I have to evaluate $\pi$ to $n$ precisions, rather than doing comparison to decide whether to end iterations, I will estimate the required number of iterations at the first place. But how? Since Gauss-Legendre algorithm has quadratic convergence, in theory, the number of precision must be doubled each iteration, but it actually seems to converge faster: $$ 3.14\cdots \quad \text{(precision 2)} $$ $$ 3.1415926\cdots \quad \text{(precision 7)} $$ $$ 3.141592653589793238\cdots \quad \text{(precision 18)} $$
How can I estimate the numbers of precision given that there has been $m$ iterations? Not asymptotically, I need a formula in closed form.