[Math] Rational approximation of square roots

approximationdiophantine-approximation

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

Related Question