math/mp-sqrt.c: Explain the algorithm, and particularly the end condition.

The gnomic remark that `Increasing x is pointless when -q < 2 x + 1' was
enough to lead me back on the right lines, but is hardly adequate.  I
ended up wasting quite a lot of time with a whiteboard, because the
`... + 1' makes the termination condition look just enough like `the
fraction comes out zero'.  But no: that's (not quite) a
coincidence.  (Thinking more carefully, it's no surprise that the delta
looks similar to the derivative, but it's definitely the former we're
interested in here, rather than the latter.)

diff --git a/math/mp-sqrt.c b/math/mp-sqrt.c
index 2dbe418..bf19fe8 100644 (file)
--- a/math/mp-sqrt.c
+++ b/math/mp-sqrt.c
@@ -71,12 +71,32 @@ mp *mp_sqrt(mp *d, mp *a)
 
   /* --- Main approximation --- *
    *
-   * We use the Newton-Raphson recurrence relation
+   * The Newton--Raphson method finds approximate zeroes of a function by
+   * starting with a guess and repeatedly refining the guess by approximating
+   * the function near the guess by its tangent at the guess
+   * %$x$%-coordinate, using where the tangent cuts the %$x$%-axis as the new
+   * guess.
    *
-   *   %$x_{i+1} = x_i - \frac{x_i^2 - a}{2 x_i}$%
+   * Given a function %$f(x)$% and a guess %$x_i$%, the tangent has the
+   * equation %$y = f(x_i) + f'(x_i) (x - x_i)$%, which is zero when
    *
-   * We inspect the term %$q = x^2 - a$% to see when to stop.  Increasing
-   * %$x$% is pointless when %$-q < 2 x + 1$%.
+   *   %$\displaystyle x = x_i - \frac{f(x_i)}{f'(x_i)}
+   *
+   * We set %$f(x) = x^2 - a$%, so our recurrence will be
+   *
+   *   %$\displaystyle x_{i+1} = x_i - \frac{x_i^2 - a}{2 x_i}$%
+   *
+   * It's possible to simplify this, but it's useful to see %$q = x_i^2 - a$%
+   * so that we know when to stop.  We want the largest integer not larger
+   * than the true square root.  If %$q > 0$% then %$x_i$% is definitely too
+   * large, and we should decrease it by at least one even if the adjustment
+   * term %$(x_i^2 - a)/2 x$% is less than one.
+   *
+   * Suppose, then, that %$q \le 0$%.  Then %$(x_i + 1)^2 - a = {}$%
+   * $%x_i^2 + 2 x_i + 1 - a = q + 2 x_i + 1$%.  Hence, if %$q \ge -2 x_i$%
+   * then %$x_i + 1$% is an overestimate and we should settle where we are.
+   * Otherwise, %$x_i + 1$% is an underestimate -- but, in this case the
+   * adjustment will always be at least one.
    */
 
   for (;;) {