Rational approximation of square roots

Question

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.

Answer 1

The typical approach here is probably to use Newton's method (https://en.wikipedia.org/wiki/Newton's_method), which for this case gives:

$k_{n+1} = k_n - \frac{k_n^2-N}{2k_n}$

Since you are solving a quadratic this will always converge to the correct root (or it's negative) no matter what starting value you give, however it will be faster if you pick an initial value nearby. If you're able to, you could try combining a lookup table with linear interpolation to get a reasonable initial guess (but maybe that requires too much utility?)

Answer 2

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}$$

Answer 3

Here are some of my own, they don't require iteration; just enter in an $n,$ they are extremely accurate; by some small n it should already be greater than your calculator will display and increasing to limit of $0$ error at infinity, and they meet your criteria.

$\sqrt{n}\sqrt{(n+1)}\approx\frac{2n(n+1)}{2n+1}+\frac{2n+1}{4(2n+1)^{2}+1}.$

$\sqrt{n^{2}+1}\approx\frac{2n(n^{2}+1)}{2n^{2}+1}+\frac{2n^{2}+1}{n(4(2n^{2}+1)^{2}+1)}.$

$\sqrt{n^{2}-1}\approx\frac{2n(n^{2}-1)}{2n^{2}-1}+\frac{2n^{2}-1}{n(4(2n^{2}-1)^{2}+1)}.$

$\sqrt{n}\sqrt{(2n+1)}\approx(\frac{1}{\sqrt{2}})(\frac{4n(2n+1)}{4n+1}+\frac{4n+1}{4(4n+1)^{2}+1}).$

$\sqrt{n}\sqrt{(2n-1)}\approx(\frac{1}{\sqrt{2}})(\frac{4n(2n-1)}{4n-1}+\frac{4n-1}{4(4n-1)^{2}+1}).$

At the very least one of these may be a good start for K. As an example for $n=35,$ calculation of $\sqrt{35^{2}+1}$ using the second formula above is correct to 19 decimal places, which is already past your convergence margin of $10$-$14$ decimal places in one simple fraction. And, like I said, that increases for higher $n$.

The only downfall is none of these contain just $n.$

If you need more info on their derivation you can look at my post. "Close approximation of absolute value function"