I have had a recent interest in rational numbers. These numbers are those which can be expressed as a ratio of integers. Compared to real numbers they have many attractivie qualities, conceptual ease is one, but also many problems which are non-computable over real numbers become easily computable when restricted to rational numbers (e.g. given any two rational numbers I can tell you in a finite amount of time whether one is less than the other.) Yet for all their advantages, they are also deeply problematic. Calculus simply will not work on rational numbers, at least not in its traditional formulation. As a computational scientist this is a heavy blow because real analysis is a bedrock tool we use in order to prove correctness of many algorithms. Is this a fixable situation? Can we keep the nice decidable properties of rational numbers while also doing useful analytic arguments with them?
I have decided to investigate this for my own personal curiosity. I show below a result for using Newton’s algorithm to compute the square root of a number.
This algorithm computes the square root of a number \( c \) by taking any starting guess value \( x0 \) and repeatedly iterating the sequence:
\[ x _ {n+1} = \frac{ x _ {n} }{2} + \frac{c}{2x _ {n}}. \]
The way we know this computes the right thing is usually by appealing to real analysis, which allows us to write a theorem such as:
\[ \lim _{n \to \infty} x_n = \sqrt{c} \]
A convergence rate should accompany this result to make it more useful, but this is the basic statement that our algorithm is correct. There are a few problems with this approach however when we are limited to rational numbers.
There are a few things wrong with the correctness result stated above as a limit. One issue is that the vast majority of rational numbers do not have a square root (google “square root of 2 is irrational” and pick your favorite proof), so if we are to limit ourselves to rational numbers then the symbol \( \sqrt{c} \) has no meaning.
Another issue is the use of the limit. If we appeal back to the foundational definitions of real analysis, the very idea of a limit heavily depends on completeness properties of the real number system. Thus if we are to reproduce a correctness result for rational numbers we can not use the limit either.
Fortunately we may still make sense of the square root if we understand it as an approximation rather than as an exact number. If instead of asking exactly for the real number \( \sqrt{c} \) I ask for a rational number \( r \) such that \( r^2 - c \) is below some error threshold \( \epsilon \), then I may use \(r \) very much like I would use its likely-irrational cousin \( \sqrt{c} \) just I would have to track the error \(\epsilon \) in any resulting expression using it.
The question then is how do we know if \( r^2 - c \) is beneath our chosen error threshold? I give below a fully computable error estimate that does not rely on calculus.
A simple algebraic trick reveals that no matter what starting point we use in the Newton iteration, the resulting error \( x _ n ^ 2 - c \) is always nonnegative (for \( n>=1 \) ). This follows from
This lemma shows that the error is bounded below by \( 0 \), the next lemma shows how we may bound the error from above
Now that we know the Newton iteration error does not careen off into negative numbers, we must show that it is bounded above as well. Ideally this bound will get smaller as we increase \(n \), indicating successively better approximations. Indeed it does, asusming that \(n>1\) we have
The inequality above arises from the fact that \(x _ n^2 - c \) is nonnegative, proved in Lemma 1.
A simple mathematical induction argument combined with Lemma 2 and Lemma 1 shows that
\[ 0 \leq x _ {n+1}^2 - c \leq \left (\frac{1}{4} \right) ^n \left (x _ {1} ^ 2 - c \right) \]
I have formalized the above argument in Coq using QArith, Coq’s formalized rational number library. I won’t go into great detail here how I formalized it, but you may find the coq file on my github.
I will just end with three simple computations done in Coq using my formalization of the above arguments. The input to the routine is the rational number \( c \) you wish to have an approximate square root for, and the error tolerance \( \epsilon \) you want to satisfy. The smaller your tolerance, the more expensive the square root becomes to compute.
Eval compute in (sqrt_with_err (2#1) (1#10)).
(*
result: 544/384
*)
The result is
\[ \frac{544}{384} \]
Eval compute in (sqrt_with_err (2#1) (1#100)).
(*
5585613357056 / 3949625081856
*)
The result is
\[ \frac{5585613357056}{3949625081856} \]
Eval compute in (sqrt_with_err (2#1) (1#1000)).
(*
124796306297948855091593216 / 88244314450313083820703744
*)
The result is
\[ \frac{124796306297948855091593216}{88244314450313083820703744} \]
Eval compute in (sqrt_with_err (2#1) (1#10000)).
(*
15523401825382753748015991456494612479890246964696135663503014851104858287082550715570437825607166429495296/
10976702697811775276479974021712309192862272555841355215716474568827946263816058771853758825719716730896384
*)
The result is
\[ \frac{15523401825382753748015991456494612479890246964696135663503014851104858287082550715570437825607166429495296}{10976702697811775276479974021712309192862272555841355215716474568827946263816058771853758825719716730896384}\]
I started this with the thought of bringing real-analytic techniques to rational numbers so as to make them more useful for proving properties about iterative algorithms, but it is quite a bit more challenging than with real numbers and I am still on the fence about whether the extra work is worth it. It is something I have to think a bit more about.
Big thanks to my good friend Tiago Cogumbreiro for teaching me Coq and guiding me through the world of functional programming in general. Tiago helped me with the square root Coq formalization (but any unpleasantness you may find in that code likely has its origin in me rather than Tiago).