[Gmp-commit] /var/hg/gmp-proj/mini-gmp: 2 new changesets

[mercurial at gmplib.org]

details:   /var/hg/gmp-proj/mini-gmp/rev/ee0425fe3e0f
changeset: 46:ee0425fe3e0f
user:      Niels M?ller <nisse at lysator.liu.se>
date:      Wed Jan 04 21:56:54 2012 +0100
description:
Strike mpz_powm from TODO list.

details:   /var/hg/gmp-proj/mini-gmp/rev/d74bc57c9fa5
changeset: 47:d74bc57c9fa5
user:      Niels M?ller <nisse at lysator.liu.se>
date:      Wed Jan 04 23:42:06 2012 +0100
description:
Implemented mpz_sqrtrem.

diffstat:

 .hgignore      |   1 +
 mini-gmp.c     |  59 +++++++++++++++++++++++++++++++++++++++++++--
 mini-gmp.h     |   2 +
 tests/Makefile |   2 +-
 tests/t-sqrt.c |  75 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 135 insertions(+), 4 deletions(-)

diffs (190 lines):

diff -r b753004ece2c -r d74bc57c9fa5 .hgignore
--- a/.hgignore	Wed Jan 04 21:25:37 2012 +0100
+++ b/.hgignore	Wed Jan 04 23:42:06 2012 +0100
@@ -10,6 +10,7 @@
 ^tests/t-div$
 ^tests/t-double$
 ^tests/t-gcd$
+^tests/t-sqrt
 ^tests/t-powm$
 ^tests/t-logops$
 ^tests/t-str$
diff -r b753004ece2c -r d74bc57c9fa5 mini-gmp.c
--- a/mini-gmp.c	Wed Jan 04 21:25:37 2012 +0100
+++ b/mini-gmp.c	Wed Jan 04 23:42:06 2012 +0100
@@ -36,9 +36,6 @@
      mpz_hamdist
      mpz_lcm_ui
      mpz_popcount
-     mpz_powm
-     mpz_sqrtrem
-     mpz_tstbit
  */
 #include <assert.h>
 #include <ctype.h>
@@ -2347,6 +2344,62 @@
   return invertible;
 }
 
+/* Compute s = floor(sqrt(u)) and r = u - s^2 */
+void
+mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
+{
+  mpz_t x, t, dx, q;
+
+  if (u->_mp_size < 0)
+    gmp_die ("mpz_sqrtrem: Negative argument.");
+  if (u->_mp_size == 0)
+    {
+      mpz_set_ui (s, 0);
+      mpz_set_ui (r, 0);
+      return;
+    }
+
+  mpz_init (x);
+  mpz_init (t);
+  mpz_init (dx);
+  mpz_init (q);
+
+  /* Make x > sqrt(a). This will be invariant through the loop. */
+  mpz_setbit (x, (mpz_sizeinbase (u, 2) + 1) / 2);
+
+  for (;;)
+    {
+      /* Compute x^2 - u */
+      mpz_mul (t, x, x);
+      mpz_sub (t, t, u);
+      assert (t->_mp_size > 0);
+
+      mpz_mul_2exp (dx, x, 1);
+      mpz_tdiv_q (q, t, dx);
+      if (q->_mp_size == 0)
+	break;
+      mpz_sub (x, x, q);
+      assert (x->_mp_size > 0);
+    }
+  /* We have 0 < u - x^2 < 2x, which implies that sqrt(a) < x < 1 +
+     sqrt(1+a), and x-1 = floor(sqrt(a)). Then r = a - (x-1)^2 = 2x -
+     1 - t. */
+  mpz_sub_ui (x, x, 1);
+  assert (x->_mp_size > 0);
+
+  mpz_sub_ui (dx, dx, 1);
+  mpz_sub (t, dx, t);
+  assert (t->_mp_size >= 0);
+
+  mpz_swap (s, x);
+  mpz_swap (t, r);
+
+  mpz_clear (x);
+  mpz_clear (t);
+  mpz_clear (dx);
+  mpz_clear (q);
+}
+
 void
 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
 {
diff -r b753004ece2c -r d74bc57c9fa5 mini-gmp.h
--- a/mini-gmp.h	Wed Jan 04 21:25:37 2012 +0100
+++ b/mini-gmp.h	Wed Jan 04 23:42:06 2012 +0100
@@ -136,6 +136,8 @@
 void mpz_gcdext (mpz_t, mpz_t, mpz_t, const mpz_t, const mpz_t);
 int mpz_invert (mpz_t, const mpz_t, const mpz_t);
 
+void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t);
+
 void mpz_powm (mpz_t, const mpz_t, const mpz_t, const mpz_t);
 
 int mpz_tstbit (const mpz_t, mp_bitcnt_t);
diff -r b753004ece2c -r d74bc57c9fa5 tests/Makefile
--- a/tests/Makefile	Wed Jan 04 21:25:37 2012 +0100
+++ b/tests/Makefile	Wed Jan 04 23:42:06 2012 +0100
@@ -7,7 +7,7 @@
 LIBS = -lgmp -lm -lmcheck
 
 CHECK_PROGRAMS = t-add t-sub t-mul t-invert t-div t-double t-gcd \
-	t-powm t-logops t-str 
+	t-sqrt t-powm t-logops t-str 
 
 MISC_OBJS = hex-random.o mini-random.o mini-gmp.o
 
diff -r b753004ece2c -r d74bc57c9fa5 tests/t-sqrt.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/t-sqrt.c	Wed Jan 04 23:42:06 2012 +0100
@@ -0,0 +1,75 @@
+#include <limits.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "mini-random.h"
+
+#define MAXBITS 400
+#define COUNT 10000
+
+static void
+dump (const char *label, const mpz_t x)
+{
+  char *buf = mpz_get_str (NULL, 16, x);
+  fprintf (stderr, "%s: %s\n", label, buf);
+  free (buf);
+}
+
+/* Called when s is supposed to be floor(sqrt(u)), and r = u - s^2 */
+static int
+sqrtrem_valid_p (const mpz_t u, const mpz_t s, const mpz_t r)
+{
+  mpz_t t;
+
+  mpz_init (t);
+  mpz_mul (t, s, s);
+  mpz_sub (t, u, t);
+  if (mpz_sgn (t) < 0 || mpz_cmp (t, r) != 0)
+    {
+      mpz_clear (t);
+      return 0;
+    }
+  mpz_add_ui (t, s, 1);
+  mpz_mul (t, t, t);
+  if (mpz_cmp (t, u) <= 0)
+    {
+      mpz_clear (t);
+      return 0;
+    }
+  
+  mpz_clear (t);
+  return 1;  
+}
+
+int
+main (int argc, char **argv)
+{
+  unsigned i;
+  mpz_t u, s, r;
+  
+  hex_random_init ();
+
+  mpz_init (u);
+  mpz_init (s);
+  mpz_init (r);
+  
+  for (i = 0; i < COUNT; i++)
+    {
+      mini_rrandomb (u, MAXBITS);
+      mpz_sqrtrem (s, r, u);
+
+      if (!sqrtrem_valid_p (u, s, r))
+	{
+	  fprintf (stderr, "mpz_sqrtrem failed:\n");
+	  dump ("u", u);
+	  dump ("sqrt", s);
+	  dump ("rem", r);
+	  abort ();
+	}
+    }
+  mpz_clear (u);
+  mpz_clear (s);
+  mpz_clear (r);
+
+  return 0;
+}