[Gmp-commit] /var/hg/gmp: mini-gmp: add mpz_rootrem, correct mpz_root and tests.

mercurial at gmplib.org mercurial at gmplib.org
Wed Apr 11 09:05:07 CEST 2012


details:   /var/hg/gmp/rev/c859958b06d1
changeset: 14804:c859958b06d1
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Wed Apr 11 09:04:59 2012 +0200
description:
mini-gmp: add mpz_rootrem, correct mpz_root and tests.

diffstat:

 ChangeLog                |  10 ++++++
 mini-gmp/mini-gmp.c      |  34 ++++++++++++++++---
 mini-gmp/mini-gmp.h      |   3 +-
 mini-gmp/tests/Makefile  |   2 +-
 mini-gmp/tests/t-reuse.c |   4 +-
 mini-gmp/tests/t-root.c  |  80 ++++++++++++++++++++++++++++++++++++++++++++++++
 6 files changed, 123 insertions(+), 10 deletions(-)

diffs (213 lines):

diff -r b04e4e6b6c95 -r c859958b06d1 ChangeLog
--- a/ChangeLog	Tue Apr 10 22:32:14 2012 +0200
+++ b/ChangeLog	Wed Apr 11 09:04:59 2012 +0200
@@ -1,3 +1,13 @@
+2012-04-11 Marco Bodrato <bodrato at mail.dm.unipi.it>
+
+	* mini-gmp/mini-gmp.h (mpz_root, mpz_rootrem): define (correctly).
+	* mini-gmp/mini-gmp.c (mpz_rootrem): Extended code from _root.
+	(mpz_root): Use mpz_rootrem.
+	
+	* mini-gmp/tests/Makefile (CHECK_PROGRAMS): add t-root.
+	* mini-gmp/tests/t-root.c: New file
+	* mini-gmp/tests/t-reuse.c: Enable root{,rem} tests.
+
 2012-04-10 Marco Bodrato <bodrato at mail.dm.unipi.it>
 
 	* gen-fac_ui.c (mpz_root): Remove.
diff -r b04e4e6b6c95 -r c859958b06d1 mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c	Tue Apr 10 22:32:14 2012 +0200
+++ b/mini-gmp/mini-gmp.c	Wed Apr 11 09:04:59 2012 +0200
@@ -3143,19 +3143,23 @@
   mpz_clear (e);
 }
 
-/* x=floor(y^(1/z)) */
+/* x=floor(y^(1/z)), r=y-x^z */
 void
-mpz_root (mpz_t x, const mpz_t y, unsigned long z)
+mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
 {
+  int sgn;
   mpz_t t, u, v, ty;
 
-  if (y->_mp_size < 0 && (z & 1) == 0)
-    gmp_die ("mpz_root: Negative argument, with even root.");
+  sgn = y->_mp_size < 0;
+  if (sgn && (z & 1) == 0)
+    gmp_die ("mpz_rootrem: Negative argument, with even root.");
   if (z == 0)
-    gmp_die ("mpz_root: Zeroth root.");
+    gmp_die ("mpz_rootrem: Zeroth root.");
 
   if (mpz_cmp_ui (y, 1) <= 0) {
     mpz_set (x, y);
+    if (r)
+      r->_mp_size = 0;
     return;
   }
 
@@ -3176,7 +3180,11 @@
     mpz_tdiv_q_ui (t, t, z);
   } while (mpz_cmp (t, u) < 0);
 
-  if (y->_mp_size < 0)
+  if (r) {
+    mpz_pow_ui (t, u, z);
+    mpz_sub (r, y, t);
+  }
+  if (sgn)
     mpz_neg (x, u);
   else
     mpz_set (x, u);
@@ -3185,6 +3193,20 @@
   mpz_clear (t);
 }
 
+int
+mpz_root (mpz_t x, const mpz_t y, unsigned long z)
+{
+  int res;
+  mpz_t r;
+
+  mpz_init (r);
+  mpz_rootrem (x, r, y, z);
+  res = r->_mp_size == 0;
+  mpz_clear (r);
+
+  return res;
+}
+
 
 /* Logical operations and bit manipulation. */
 
diff -r b04e4e6b6c95 -r c859958b06d1 mini-gmp/mini-gmp.h
--- a/mini-gmp/mini-gmp.h	Tue Apr 10 22:32:14 2012 +0200
+++ b/mini-gmp/mini-gmp.h	Wed Apr 11 09:04:59 2012 +0200
@@ -177,7 +177,8 @@
 void mpz_powm (mpz_t, const mpz_t, const mpz_t, const mpz_t);
 void mpz_powm_ui (mpz_t, const mpz_t, unsigned long, const mpz_t);
 
-void mpz_root (mpz_t, const mpz_t, unsigned long);
+void mpz_rootrem (mpz_t, mpz_t, const mpz_t, unsigned long);
+int mpz_root (mpz_t, const mpz_t, unsigned long);
 
 int mpz_tstbit (const mpz_t, mp_bitcnt_t);
 void mpz_setbit (mpz_t, mp_bitcnt_t);
diff -r b04e4e6b6c95 -r c859958b06d1 mini-gmp/tests/Makefile
--- a/mini-gmp/tests/Makefile	Tue Apr 10 22:32:14 2012 +0200
+++ b/mini-gmp/tests/Makefile	Wed Apr 11 09:04:59 2012 +0200
@@ -12,7 +12,7 @@
 
 CHECK_PROGRAMS = t-add t-sub t-mul t-invert t-div t-div_2exp \
 	t-double t-gcd t-lcm \
-	t-sqrt t-powm t-logops t-bitops t-scan t-str \
+	t-sqrt t-root t-powm t-logops t-bitops t-scan t-str \
 	t-reuse
 
 MISC_OBJS = hex-random.o mini-random.o mini-gmp.o
diff -r b04e4e6b6c95 -r c859958b06d1 mini-gmp/tests/t-reuse.c
--- a/mini-gmp/tests/t-reuse.c	Tue Apr 10 22:32:14 2012 +0200
+++ b/mini-gmp/tests/t-reuse.c	Wed Apr 11 09:04:59 2012 +0200
@@ -361,7 +361,7 @@
 	  if (mpz_cmp (ref1, res1) != 0 || mpz_cmp (ref2, res2) != 0)
 	    FAIL2 (mpz_sqrtrem, in1, NULL, NULL);
 	}
-#if 0
+
       if (mpz_sgn (in1) >= 0)
 	{
 	  mpz_root (ref1, in1, in2i % 0x1000 + 1);
@@ -394,7 +394,7 @@
 	  if (mpz_cmp (ref1, res1) != 0 || mpz_cmp (ref2, res2) != 0)
 	    FAIL2 (mpz_rootrem, in1, in2, NULL);
 	}
-#endif
+
       if (pass < reps / 2)	/* run fewer tests since gcdext lots of time */
 	{
 	  mpz_gcdext (ref1, ref2, ref3, in1, in2);
diff -r b04e4e6b6c95 -r c859958b06d1 mini-gmp/tests/t-root.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mini-gmp/tests/t-root.c	Wed Apr 11 09:04:59 2012 +0200
@@ -0,0 +1,80 @@
+#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(root(u,z)), and r = u - s^z */
+static int
+rootrem_valid_p (const mpz_t u, const mpz_t s, const mpz_t r, unsigned long z)
+{
+  mpz_t t;
+
+  mpz_init (t);
+  mpz_pow_ui (t, s, z);
+  mpz_sub (t, u, t);
+  if (mpz_sgn (t) != mpz_sgn(u) || mpz_cmp (t, r) != 0)
+    {
+      mpz_clear (t);
+      return 0;
+    }
+  mpz_add_ui (t, s, 1);
+  mpz_pow_ui (t, t, z);
+  if (mpz_cmpabs (t, u) <= 0)
+    {
+      mpz_clear (t);
+      return 0;
+    }
+
+  mpz_clear (t);
+  return 1;
+}
+
+int
+main (int argc, char **argv)
+{
+  unsigned i;
+  unsigned long e;
+  mpz_t u, s, r, bs;
+
+  hex_random_init ();
+
+  mpz_init (u);
+  mpz_init (s);
+  mpz_init (r);
+  mpz_init (bs);
+
+  for (i = 0; i < COUNT; i++)
+    {
+      mini_rrandomb (u, MAXBITS);
+      mini_rrandomb (bs, 10);
+      e = mpz_getlimbn (bs, 0) % mpz_sizeinbase (u, 2) + 2;
+      mpz_rootrem (s, r, u, e);
+
+      if (!rootrem_valid_p (u, s, r, e))
+	{
+	  fprintf (stderr, "mpz_rootrem(%lu-th) failed:\n", e);
+	  dump ("u", u);
+	  dump ("root", s);
+	  dump ("rem", r);
+	  abort ();
+	}
+    }
+  mpz_clear (bs);
+  mpz_clear (u);
+  mpz_clear (s);
+  mpz_clear (r);
+
+  return 0;
+}


More information about the gmp-commit mailing list