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

diff --git a/math/mp-sqrt.c b/math/mp-sqrt.c
index 2dbe4189..9530ff99 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 (;;) {
@@ -106,6 +126,24 @@ mp *mp_sqrt(mp *d, mp *a)
   return (d);
 }
 
+/* --- @mp_squarep@ --- *
+ *
+ * Arguments:	@mp *n@ = an integer
+ *
+ * Returns:	Nonzero if and only if @n@ is a perfect square, i.e.,
+ *		%$n = a^2$% for some rational integer %$a$%.
+ */
+
+int mp_squarep(mp *n)
+{
+  mp *t = MP_NEW;
+  int rc;
+
+  if (MP_NEGP(n)) return (0);
+  t = mp_sqrt(t, n); t = mp_sqr(t, t);
+  rc = MP_EQ(t, n); mp_drop(t); return (rc);
+}
+
 /*----- Test rig ----------------------------------------------------------*/
 
 #ifdef TEST_RIG