X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/0a80c8cefb56e80ccb95277b250dab0c10e99d9d..d9383a70ff0c0a2f25f8e526d275ac650f82390f:/math/mp-sqrt.c

diff --git a/math/mp-sqrt.c b/math/mp-sqrt.c
index 2dbe4189..bf19fe8e 100644
--- 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 (;;) {