*** empty log message ***

2 files changed, 168 insertions, 0 deletions

diff --git a/ChangeLog b/ChangeLog
index f0e203994..7e083ddc9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -19,6 +19,16 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA.
 
 
+2002-05-08  Torbjorn Granlund  <tege@swox.com>
+
+	* mpn/generic/rootrem.c: New file.
+
+	* mpz/perfpow.c: Amend allocations for mpn_rootrem requirements.
+
+	* gmp-impl.h (mpn_rootrem): Add declaration.
+
+	* mpz/root.c: Rewrite to use mpn_rootrem.
+
 2002-05-08  Kevin Ryde  <kevin@swox.se>
 
 	* gmp-impl.h (MUL_KARATSUBA_THRESHOLD etc): Remove forced nail values.
diff --git a/mpn/generic/rootrem.c b/mpn/generic/rootrem.c
new file mode 100644
index 000000000..7c2e50b0d
--- /dev/null
+++ b/mpn/generic/rootrem.c
@@ -0,0 +1,158 @@
+/* mpn_rootrem(rootp,remp,ap,an,nth) -- Compute the nth root of {ap,an}, and
+   store the truncated integer part at rootp and the remainder at remp.
+
+   THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
+   INTERFACES.  IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
+   IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
+   FUTURE GNU MP RELEASE.
+
+
+Copyright 2002 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+/*
+  We use Newton's method to compute the root of a:
+
+           n
+  f(x) := x  - a
+
+
+            n - 1
+  f'(x) := x      n
+
+
+                                       n-1            n-1           n-1
+                                x - a/x            a/x   - x     a/x   + (n-1)x
+  new x = x - f(x)/f'(x) =  x - ----------  =  x + ---------  =  --------------
+                                     n                 n                n
+
+*/
+
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+mp_size_t
+mpn_rootrem (mp_ptr rootp, mp_ptr remp,
+	     mp_srcptr up, mp_size_t un, mp_limb_t nth)
+{
+  mp_ptr pp, qp;
+  mp_size_t pn, xn, qn;
+  unsigned long int unb, xnb, bit;
+  unsigned int cnt;
+  mp_size_t i;
+  unsigned long int n_valid_bits, adj;
+  TMP_DECL (marker);
+
+  TMP_MARK (marker);
+  pp = TMP_ALLOC_LIMBS (un + 2);
+  qp = TMP_ALLOC_LIMBS (un);	/* should try and reduce this allocation */
+
+  count_leading_zeros (cnt, up[un - 1]);
+  unb = un * GMP_NUMB_BITS - cnt + GMP_NAIL_BITS;
+
+  xnb = (unb + nth - 1) / nth;
+  if (xnb == 1)
+    {
+      rootp[0] = 1;
+      if (remp == NULL)
+	remp = pp;
+      mpn_sub_1 (remp, up, un, (mp_limb_t) 1);
+      MPN_NORMALIZE (remp, un);
+      TMP_FREE (marker);
+      return un;
+    }
+
+  xn = (xnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
+
+  /* Set initial root to only ones.  This is an overestimate of the actual root
+     by less than a factor of 2.  */
+  for (i = 0; i < xn; i++)
+    rootp[i] = GMP_NUMB_MAX;
+  rootp[xnb / GMP_NUMB_BITS] = ((mp_limb_t) 1 << (xnb % GMP_NUMB_BITS)) - 1;
+
+  /* Improve the initial approximation, one bit at a time.  Keep the
+     approximations >= root(U,nth).  */
+  bit = xnb - 2;
+  n_valid_bits = 0;
+  for (i = 0; (nth >> i) != 0; i++)
+    {
+      mp_limb_t xl = rootp[bit / GMP_NUMB_BITS];
+      rootp[bit / GMP_NUMB_BITS] = xl ^ (mp_limb_t) 1 << bit % GMP_NUMB_BITS;
+      pn = mpn_pow_1 (pp, rootp, xn, nth, qp);
+      /* If the new root approximation is too small, restore old value.  */
+      if (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0)))
+	rootp[bit / GMP_NUMB_BITS] = xl;		/* restore old value */
+      n_valid_bits += 1;
+      if (bit == 0)
+	goto done;
+      bit--;
+    }
+
+  adj = n_valid_bits - 1;
+
+  /* Newton loop.  Converges downwards towards root(U,nth).  Currently we use
+   full precision from iteration 1.  Clearly, we should use just n_valid_bits
+   of precision in each step, and thus save most of the computations.  */
+  while (n_valid_bits <= xnb)
+    {
+      mp_limb_t cy;
+
+      pn = mpn_pow_1 (pp, rootp, xn, nth - 1, qp);
+      mpn_tdiv_qr (qp, pp, (mp_size_t) 0, up, un, pp, pn); /* junk remainder */
+      qn = un - pn + 1;
+      cy = mpn_addmul_1 (qp, rootp, xn, nth - 1);
+      if (qn == xn + 1)
+	{
+	  cy += qp[xn];
+	  if (cy == nth)
+	    {
+	      for (i = xn - 1; i >= 0; i--)
+		qp[i] = GMP_NUMB_MAX;
+	      cy = nth - 1;
+	    }
+	}
+
+      qp[xn] = cy;
+      qn = xn + (cy != 0);
+
+      mpn_divrem_1 (rootp, (mp_size_t) 0, qp, qn, nth);
+      n_valid_bits = n_valid_bits * 2 - adj;
+    }
+
+  /* The computed result might be one unit too large.  Adjust as necessary.  */
+ done:
+  pn = mpn_pow_1 (pp, rootp, xn, nth, qp);
+  if (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0))
+    {
+      mpn_decr_u (rootp, 1);
+      pn = mpn_pow_1 (pp, rootp, xn, nth, qp);
+
+      ASSERT_ALWAYS (! (un < pn || (un == pn && mpn_cmp (up, pp, pn) < 0)));
+    }
+
+  if (remp == NULL)
+    remp = pp;
+  mpn_sub (remp, up, un, pp, pn);
+  MPN_NORMALIZE (remp, un);
+
+  TMP_FREE (marker);
+  return un;
+}