[Gmp-commit] /var/hg/gmp: 3 new changesets

[mercurial at gmplib.org]

details:   /var/hg/gmp/rev/ce95fa86a330
changeset: 14214:ce95fa86a330
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Wed Aug 03 21:53:51 2011 +0200
description:
(mpz_tdiv_qr): Correctly handle dividend value being equal to divisor value.

details:   /var/hg/gmp/rev/5d784d9cb4f2
changeset: 14215:5d784d9cb4f2
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Wed Aug 03 22:02:46 2011 +0200
description:
(mpz_root): Create reasonable starting approximation.
(mpz_sqrt): New function.
(mpz_mul_2exp): Add faster block shifting code, disabled for now.

details:   /var/hg/gmp/rev/271cd7af4f72
changeset: 14216:271cd7af4f72
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Wed Aug 03 22:03:03 2011 +0200
description:
*** empty log message ***

diffstat:

 ChangeLog |   8 ++++++
 dumbmp.c  |  78 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----
 2 files changed, 81 insertions(+), 5 deletions(-)

diffs (125 lines):

diff -r 17adfcccc3db -r 271cd7af4f72 ChangeLog
--- a/ChangeLog	Thu Jul 14 15:25:59 2011 +0200
+++ b/ChangeLog	Wed Aug 03 22:03:03 2011 +0200
@@ -1,3 +1,11 @@
+2011-08-03  Torbjorn Granlund  <tege at gmplib.org>
+
+	* dumbmp.c (mpz_tdiv_qr): Correctly handle dividend value being equal
+	to divisor value.
+	(mpz_root): Create reasonable starting approximation.
+	(mpz_sqrt): New function.
+	(mpz_mul_2exp): Add faster block shifting code, disabled for now.
+
 2011-07-14  Torbjorn Granlund  <tege at gmplib.org>
 
 	* mpn/generic/dcpi1_bdiv_q.c (mpn_dcpi1_bdiv_q): Get mpn_sub_1 size
diff -r 17adfcccc3db -r 271cd7af4f72 dumbmp.c
--- a/dumbmp.c	Thu Jul 14 15:25:59 2011 +0200
+++ b/dumbmp.c	Wed Aug 03 22:03:03 2011 +0200
@@ -517,6 +517,34 @@
 mpz_mul_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
 {
   mpz_set (r, a);
+
+#if 0
+  {
+    unsigned long int lcnt;
+
+    lcnt = bcnt / GMP_LIMB_BITS;
+    if (lcnt > 0)
+      {
+	int rn = ABSIZ (r);
+	mp_limb_t *rp;
+	int i;
+
+	mpz_realloc (r, rn + lcnt);
+	rp = PTR (r);
+
+	for (i = rn - 1; i >= 0; i--)
+	  rp[i + lcnt] = rp[i];
+	for (i = lcnt - 1; i >= 0; i--)
+	  rp[i] = 0;
+
+	rn += lcnt;
+	SIZ (r) = SIZ (r) >= 0 ? rn : -rn;
+
+	bcnt %= GMP_LIMB_BITS;
+      }
+  }
+#endif
+
   while (bcnt)
     {
       mpz_add (r, r, r);
@@ -658,7 +686,7 @@
   mpz_init_set (tmpb, b);
   mpz_set_ui (q, 0L);
 
-  if (mpz_cmp (tmpr, tmpb) > 0)
+  if (mpz_cmp (tmpr, tmpb) >= 0)
     {
       cnt = mpz_sizeinbase (tmpr, 2) - mpz_sizeinbase (tmpb, 2) + 1;
       mpz_mul_2exp (tmpb, tmpb, cnt);
@@ -900,9 +928,12 @@
       mpz_set (x, y);
       return;
     }
-  mpz_init (t);
-  mpz_init_set (u, y);
-  do
+
+  /* One-bit initial approximation */
+  mpz_init_set_ui (u, 1);
+  mpz_mul_2exp (u, u, ((mpz_sizeinbase (y, 2) - 1) / z) + 1);
+
+  for (;;)
     {
       mpz_pow_ui (t, u, z - 1);
       mpz_tdiv_q (t, y, t);
@@ -912,8 +943,45 @@
 	break;
       mpz_set (u, t);
     }
-  while (1);
+
   mpz_set (x, u);
   mpz_clear (t);
   mpz_clear (u);
 }
+
+/* x=floor(y^(1/2)) */
+void
+mpz_sqrt (mpz_t x, mpz_t y)
+{
+  mpz_t t, u;
+
+  if (mpz_sgn (y) < 0)
+    {
+      fprintf (stderr, "mpz_sqrt does not accept negative values\n");
+      abort ();
+    }
+  if (mpz_cmp_ui (y, 1) <= 0)
+    {
+      mpz_set (x, y);
+      return;
+    }
+  mpz_init (t);
+
+  /* One-bit initial approximation */
+  mpz_init_set_ui (u, 1);
+  mpz_mul_2exp (u, u, ((mpz_sizeinbase (y, 2) - 1) / 2) + 1);
+
+  for (;;)
+    {
+      mpz_tdiv_q (t, y, u);
+      mpz_add (t, t, u);
+      mpz_tdiv_q_2exp (t, t, 1);
+      if (mpz_cmp (t, u) >= 0)
+	break;
+      mpz_set (u, t);
+    }
+
+  mpz_set (x, u);
+  mpz_clear (t);
+  mpz_clear (u);
+}