I'm trying to find the best way to solve for rational approximations of the square root of a number, given some pretty serious constraints on the operations I can use to calculate it. My criteria for the algorithm are that it must:
- be solvable through recursion
- operate solely on ratios of integers
- use only addition, subtraction, multiplication, division or equality
(The criteria are important because I'm trying to use c++ to approximate the root at compile-time using template-meta programming, and the available operations are extremely limited)
The approach that I'm considering is supposedly based on an ancient Babylonian method and involves iteratively solving:
$$
k_{n+1} = \frac{(k_{n} + N / k_{n})}{2}
$$
Where:
- N is the number whose root we are looking for
- K is the approximation of the root
- K[0] is chosen such that the value of k^2 is less than N
So, it seems I could pretty trivially implement this algorithm if:
- I fix K{0} to be 1, which probably hurts convergence but will work for any N > 1
- If I choose a fixed recursion depth
The problems that I would like to solve are:
- Is there in general a well known method that is faster/more accurate/etc that I should use instead?
- Is there a smarter way that I can choose the initial value of K if N is known?
- Is it possible to easily choose a recursion depth that will converge to 10-14 decimal places, a priori, if N is known? OR
- Can I test for convergence somehow using only the arithmetic mentioned above and equality? (No type of comparison operations are available, but I could for example test that the difference between two ratios is exactly equal to some other ratio)
If it's not already obvious, I'm a programmer, not a mathematician, and have no formal experience with diophantine approximation.
Best Answer
Continuing from previous answers and comments, solving for $x$ equation $f(x)=x^k-a$ (for any power $k$), can be achieved with your criteria using Newton method or, faster, with its variants such as Halley, Householder or even higher order methods they do not have a name).
For example, using Halley method, $$x_{n+1} = x_n - \frac {2 f(x_n) f'(x_n)} {2 {[f'(x_n)]}^2 - f(x_n) f''(x_n)} $$ would give $$x_{n+1}=x_n\,\frac{ a (k+1)+(k-1) x_n^k}{a (k-1)+(k+1) x_n^k}$$ If $k$ is a whole number, if $a$ is at least rational and $x_0$ rational, all iterates will be rational numbers.
For illustration, let me use $k=5$ and $a=\frac {123456}{789}$ and use $x_0=2$. The very first iterates will be $$x_1=\frac{8768}{3361}$$ $$x_2=\frac{5494130652143802990362144}{2000466367873608426046147}$$