From d9383a70ff0c0a2f25f8e526d275ac650f82390f Mon Sep 17 00:00:00 2001 From: Mark Wooding Date: Mon, 12 Aug 2019 14:44:35 +0100 Subject: [PATCH] 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.) --- math/mp-sqrt.c | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) 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 (;;) { -- 2.11.0