Modular square root patch.

Robert Ostling robert at robos.org
Wed May 3 05:31:22 CEST 2006


Hi,

As part of an implementation of the quadratic sieve, I implemented the
Shanks-Tonelli algorithm for modular square roots and now I integrated it into
the GMP sources. A patch (to 4.2.1-rc) and the source file are attached, in
case you think it belongs in the general distribution.

Robert Östling

-------------- next part --------------
diff -Nurp gmp-4.2.1-rc/Makefile.am gmp-4.2.1-rc-sqrtm/Makefile.am
--- gmp-4.2.1-rc/Makefile.am	2006-04-12 18:16:01.000000000 +0200
+++ gmp-4.2.1-rc-sqrtm/Makefile.am	2006-05-03 05:07:42.000000000 +0200
@@ -174,6 +174,7 @@ MPZ_OBJECTS = mpz/abs$U.lo mpz/add$U.lo 
   mpz/root$U.lo mpz/rootrem$U.lo mpz/rrandomb$U.lo mpz/scan0$U.lo	\
   mpz/scan1$U.lo mpz/set$U.lo mpz/set_d$U.lo mpz/set_f$U.lo		\
   mpz/set_q$U.lo mpz/set_si$U.lo mpz/set_str$U.lo mpz/set_ui$U.lo	\
+  mpz/sqrtm$U.lo \
   mpz/setbit$U.lo							\
   mpz/size$U.lo mpz/sizeinbase$U.lo mpz/sqrt$U.lo			\
   mpz/sqrtrem$U.lo mpz/sub$U.lo mpz/sub_ui$U.lo mpz/swap$U.lo		\
diff -Nurp gmp-4.2.1-rc/mpz/Makefile.am gmp-4.2.1-rc-sqrtm/mpz/Makefile.am
--- gmp-4.2.1-rc/mpz/Makefile.am	2006-03-14 16:57:54.000000000 +0100
+++ gmp-4.2.1-rc-sqrtm/mpz/Makefile.am	2006-05-03 05:08:44.000000000 +0200
@@ -50,8 +50,8 @@ libmpz_la_SOURCES = aors.h aors_ui.h fit
   powm_ui.c pprime_p.c random.c random2.c \
   realloc.c realloc2.c remove.c root.c rootrem.c rrandomb.c \
   scan0.c scan1.c set.c set_d.c set_f.c set_q.c set_si.c set_str.c \
-  set_ui.c setbit.c size.c sizeinbase.c sqrt.c sqrtrem.c sub.c sub_ui.c \
-  swap.c tdiv_ui.c tdiv_q.c tdiv_q_2exp.c tdiv_q_ui.c tdiv_qr.c \
+  set_ui.c setbit.c size.c sizeinbase.c sqrt.c sqrtrem.c sqrtm.c sub.c \
+  sub_ui.c swap.c tdiv_ui.c tdiv_q.c tdiv_q_2exp.c tdiv_q_ui.c tdiv_qr.c \
   tdiv_qr_ui.c tdiv_r.c tdiv_r_2exp.c tdiv_r_ui.c tstbit.c ui_pow_ui.c \
   ui_sub.c urandomb.c urandomm.c xor.c
 
diff -Nurp gmp-4.2.1-rc/mpz/sqrtm.c gmp-4.2.1-rc-sqrtm/mpz/sqrtm.c
--- gmp-4.2.1-rc/mpz/sqrtm.c	1970-01-01 01:00:00.000000000 +0100
+++ gmp-4.2.1-rc-sqrtm/mpz/sqrtm.c	2006-05-03 05:24:59.000000000 +0200
@@ -0,0 +1,92 @@
+/* mpz_sqrtm -- modular square roots using Shanks-Tonelli
+
+Copyright 2006 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301,
+USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/* Solve the modular equatioon x^2 = n (mod p) using the Shanks-Tonelli
+ * algorihm. x will be placed in q and 1 returned if the algorithm is
+ * successful. Otherwise 0 is returned (currently in case n is not a quadratic
+ * residue mod p). A check is done if p = 3 (mod 4), in which case the root is
+ * calculated as n ^ ((p+1) / 4) (mod p).
+ *
+ * Note that currently mpz_legendre is called to make sure that n really is a
+ * quadratic residue. The check can be skipped, at the price of going into an
+ * eternal loop if called with a non-residue.
+ */
+int mpz_sqrtm(mpz_t q, const mpz_t n, const mpz_t p) {
+    mpz_t w, n_inv, y;
+    unsigned int i, s;
+    TMP_DECL;
+    TMP_MARK;
+
+    if(mpz_divisible_p(n, p)) {         /* Is n a multiple of p?            */
+        mpz_set_ui(q, 0);               /* Yes, then the square root is 0.  */
+        return 1;                       /* (special case, not caught        */
+    }                                   /* otherwise)                       */
+    if(mpz_legendre(n, p) != 1)         /* Not a quadratic residue?         */
+        return 0;                       /* No, so return error              */
+    if(mpz_tstbit(p, 1) == 1) {         /* p = 3 (mod 4) ?                  */
+        mpz_set(q, p);
+        mpz_add_ui(q, q, 1);
+        mpz_fdiv_q_2exp(q, q, 2);
+        mpz_powm(q, n, q, p);           /* q = n ^ ((p+1) / 4) (mod p)      */
+        return 1;
+    }
+    MPZ_TMP_INIT(y, 2*SIZ(p));
+    MPZ_TMP_INIT(w, 2*SIZ(p));
+    MPZ_TMP_INIT(n_inv, 2*SIZ(p));
+    mpz_set(q, p);
+    mpz_sub_ui(q, q, 1);                /* q = p-1                          */
+    s = 0;                              /* Factor out 2^s from q            */
+    while(mpz_tstbit(q, s) == 0) s++;
+    mpz_fdiv_q_2exp(q, q, s);           /* q = q / 2^s                      */
+    mpz_set_ui(w, 2);                   /* Search for a non-residue mod p   */
+    while(mpz_legendre(w, p) != -1)     /* by picking the first w such that */
+        mpz_add_ui(w, w, 1);            /* (w/p) is -1                      */
+    mpz_powm(w, w, q, p);               /* w = w^q (mod p)                  */
+    mpz_add_ui(q, q, 1);
+    mpz_fdiv_q_2exp(q, q, 1);           /* q = (q+1) / 2                    */
+    mpz_powm(q, n, q, p);               /* q = n^q (mod p)                  */
+    mpz_invert(n_inv, n, p);
+    for(;;) {
+        mpz_powm_ui(y, q, 2, p);        /* y = q^2 (mod p)                  */
+        mpz_mul(y, y, n_inv);
+        mpz_mod(y, y, p);               /* y = y * n^-1 (mod p)             */
+        i = 0;
+        while(mpz_cmp_ui(y, 1) != 0) {
+            i++;
+            mpz_powm_ui(y, y, 2, p);    /*  y = y ^ 2 (mod p)               */
+        }
+        if(i == 0) {                    /* q^2 * n^-1 = 1 (mod p), return   */
+            TMP_FREE;
+            return 1;
+        }
+        if(s-i == 1) {                  /* In case the exponent to w is 1,  */
+            mpz_mul(q, q, w);           /* Don't bother exponentiating      */
+        } else {
+            mpz_powm_ui(y, w, 1 << (s-i-1), p);
+            mpz_mul(q, q, y);
+        }
+        mpz_mod(q, q, p);               /* r = r * w^(2^(s-i-1)) (mod p)    */
+    }
+}
+
-------------- next part --------------
A non-text attachment was scrubbed...
Name: sqrtm.c
Type: text/x-csrc
Size: 4125 bytes
Desc: not available
Url : http://gmplib.org/list-archives/gmp-devel/attachments/20060503/00f40d7e/sqrtm.bin


More information about the gmp-devel mailing list