I found a picture of Evan O'Dorney's winning project that gained him first place in the Intel Science talent search. He proposed a numerical method to find the square root, that gained him $100,000 USD.

Below are some links of pictures of the poster displaying the method.

How does this numerical method work and what is the proof?

His method makes use of Moebius Transformation.

  • 1
    $\begingroup$ Since both answers missed it (and my shoulder hurts too much to answer myself), I'd like to point out that the picture suggests that rather than applying $f$ iteratively to an initial value, you instead compute the sequence of compositions $f^{(2^n)}$. $\endgroup$ – user14972 Oct 28 '12 at 2:37
  • 1
    $\begingroup$ Is this method better than the previously known methods in some way? $\endgroup$ – user66307 Jun 19 '14 at 15:47
  • $\begingroup$ @Lembik - That reward winner claim that this method is better, but i am not sure. $\endgroup$ – Victor Jun 22 '14 at 18:37
  • $\begingroup$ @Lembik: Of course not. It is obvious (and thus long known) that the continued fraction convergents follow a recurrence relation, and clearly that can be computed by repeated squaring of the recurrence matrix. I found this by myself when I was about 18, so it's nothing great if Evan found it on his own. However, I'm sure that he knew it was well known but chose to portray otherwise to the judges, and I find it fishy that the professor who gave him the problem and obviously knew the solution let him do that. $\endgroup$ – user21820 Jul 24 '15 at 6:36

Essentially if you are interesting in evaluating $\sqrt{a}$, the idea is to first find the greatest perfect square less than or equal to $a$. Say this is $b^2$ i.e. $b = \lfloor \sqrt{a} \rfloor \implies b^2 \leq a < (b+1)^2$. Then consider the function $$f(x) = b + \dfrac{a-b^2}{x+b}$$ $$f(b) = b + \underbrace{\dfrac{a-b^2}{2b}}_{\in [0,1]} \in [b,b+1]$$ $$f(f(b)) = b + \underbrace{\dfrac{a-b^2}{f(b) + b}}_{\in [0,1]} \in [b,b+1]$$ In general $$f^{(n)}(b) = \underbrace{f \circ f \circ f \circ \cdots f}_{n \text{times}}(b) = b + \dfrac{a-b^2}{f^{(n-1)}(b)+b}$$ Hence, $f^{(n)}(b) \in [b,b+1]$ always.

If $\lim\limits_{n \to \infty}f^{(n)}(b) = \tilde{f}$ exists, then $$\tilde{f} = b + \dfrac{a-b^2}{\tilde{f}+b}$$ Hence, $$\tilde{f}^2 + b \tilde{f} = b \tilde{f} + b^2 + a - b^2 \implies \tilde{f}^2 = a$$

To prove the existence of the limit look at $$(f^{(n)}(b))^2 - a = \left(b + \dfrac{a-b^2}{f^{(n-1)}(b)+b} \right)^2 - a = \dfrac{(a-b^2)(a-(f^{(n-1)}(b))^2)}{(b+f^{(n-1)}(b))^2} = k_{n-1}(a,b)((f^{(n-1)}(b))^2-a) $$ where $\vert k_{n-1}(a,b) \vert \lt1$. Hence, convergence is also guaranteed.


Note that $k_{n-1}(a,b) = \dfrac{(a-b^2)}{(b+f^{(n-1)}(b))^2} \leq \dfrac{(b+1)^2 - 1 - b^2}{(b+b)^2} = \dfrac{2b}{(2b)^2} = \dfrac1{2b}$. This can be interpreted as larger the number, faster the convergence.

Comment: This method works only when you want to find the square of a number $\geq 1$.


To complete the answer, I am adding @Hurkyl's comment. Functions of the form $$g(z) = \dfrac{c_1z+c_2}{c_3z+c_4}$$are termed Möbius transformations. With each of these Möbius transformations, we can associate a matrix $$M = \begin{bmatrix} c_1 & c_2\\ c_3 & c_4\end{bmatrix}$$ Note that the function, $$f(x) = b + \dfrac{a-b^2}{x+b} = \dfrac{bx + a}{x+b}$$ is a Möbius transformation.

Of the many advantages of the associated matrix, one major advantage is that the associate matrix for the Möbius transformation $$g^{(n)}(z) = \underbrace{g \circ g \circ \cdots \circ g}_{n \text{ times}} = \dfrac{c_1^{(n)} z + c_2^{(n)}}{c_3^{(n)} z + c_4^{(n)}}$$ is nothing but the matrix $$M^n = \begin{bmatrix}c_1 & c_2\\ c_3 & c_4 \end{bmatrix}^n = \begin{bmatrix}c_1^{(n)} & c_2^{(n)}\\ c_3^{(n)} & c_4^{(n)} \end{bmatrix}$$ (Note that $c_k^{(n)}$ is to denote the coefficient $c_k$ at the $n^{th}$ level and is not the $n^{th}$ power of $c_k$.)

Hence, the function composition is nothing but raising the matrix $M$ to the appropriate power. This can be done in a fast way since $M^n$ can be computed in $\mathcal{O}(\log_2(n))$ operations. Thereby we can compute $g^{(2^n)}(b)$ in $\mathcal{O}(n)$ operations.

  • $\begingroup$ To Marvis -Does this method garuantee that it at least gain one extra correct decimal place when it iterate once? $\endgroup$ – Victor Oct 28 '12 at 1:37
  • $\begingroup$ This converges linearly, which can be interpreted as saying that roughly each iteration gives one more digit accuracy. $\endgroup$ – user17762 Oct 28 '12 at 1:45
  • $\begingroup$ Marvis - Apreciate for your answer! $\endgroup$ – Victor Oct 28 '12 at 1:48
  • 1
    $\begingroup$ @Victor Consecutive approximations say $f^{(n)}$ and $f^{(n+1)}$ decreases the error by a factor if $\dfrac1{2b}$ at each step. So if $b \geq 5$, we can recover one digit with every recursion. If $b \geq 50$, we can recover two digits with every recursion and so on. $\endgroup$ – user17762 Oct 28 '12 at 1:54
  • $\begingroup$ Marvis - Appreciate for your edit, it will take me a extra few hour without your edit! $\endgroup$ – Victor Oct 28 '12 at 1:59

The iteration to find $\sqrt k$ is $f(x) = \frac{d x+k}{x+d}$ where $d = \lfloor \sqrt k \rfloor$. The iterations start with $x = d$.

If $x$ is a fixed point of this, $x = \frac{d x+k}{x+d}$, or $x(x+d) = dx + k$ or $x^2 = k$, so any fixed point must be the square root.

Now wee see if the iteration increases or decreases. If $y = \frac{d x+k}{x+d}$, $$y - x = \frac{d x+k}{x+d} - x = \frac{d x+k - x(x+d)}{x+d} = \frac{k - x^2}{x+d} $$ so if $x^2 < k$, $y > x$ and if $x^2 > k$, $y < x$.

Also, proceeding like analyses of Newton's method, $y^2-k = \frac{(d x+k)^2}{(x+d)^2} - k = \frac{d^2 x^2 +2 d x k + k^2 - k(x+d)^2}{(x+d)^2} = \frac{d^2 x^2 +2 d x k + k^2 - k(x^2 + 2dx +d^2)}{(x+d)^2} = \frac{d^2 x^2 +2 d x k + k^2 - kx^2 - 2dkx -kd^2)}{(x+d)^2} = \frac{d^2 x^2 + k^2 - kx^2 -kd^2)}{(x+d)^2} = \frac{d^2 (x^2-k) + k^2 - kx^2)}{(x+d)^2} = \frac{d^2 (x^2-k) - k(x^2-k))}{(x+d)^2} = \frac{(d^2-k) (x^2-k)}{(x+d)^2} = (x^2-k)\frac{d^2-k}{(x+d)^2} $.

Since $d = \lfloor \sqrt k \rfloor$, $d < \sqrt k < d+1$ or $d^2 < k < d^2 + 2d +1$ or $-2d - 1 < d^2 - k < 0$, so $|d^2-k| < 2d+1$. Using this, $|y^2-k| < |x^2-k|\frac{2d+1}{(x+d)^2}| = |x^2-k|\frac{2d+1}{x^2+2dx+d^2} $, so $|y^2-k|< |x^2-k|$, and each iteration gets closer to the square root.

Since the starting iterate is $d$, all following iterates exceed $d$ so $|y^2-k| < |x^2-k|\frac{2d+1}{(d+d)^2}| < |x^2-k|\frac{2d+1}{4d^2}| < |x^2-k|\frac{1+1/(2d)}{2d}| \le 3|x^2-k|/4$ since $d \ge 1$.

This show that the iteration converges. However, this does not show that it converges quadratically like Newton's, only that it converges linearly.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.