Add support for computing Jacobi symbols.

diff --git a/mp-jacobi.c b/mp-jacobi.c
new file mode 100644 (file)
index 0000000..69f697e
--- /dev/null
+++ b/mp-jacobi.c
@@ -0,0 +1,167 @@
+/* -*-c-*-
+ *
+ * $Id: mp-jacobi.c,v 1.1 1999/11/22 20:50:37 mdw Exp $
+ *
+ * Compute Jacobi symbol
+ *
+ * (c) 1999 Straylight/Edgeware
+ */
+
+/*----- Licensing notice --------------------------------------------------* 
+ *
+ * This file is part of Catacomb.
+ *
+ * Catacomb is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Library General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ * 
+ * Catacomb 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 Library General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Library General Public
+ * License along with Catacomb; if not, write to the Free
+ * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+ * MA 02111-1307, USA.
+ */
+
+/*----- Revision history --------------------------------------------------* 
+ *
+ * $Log: mp-jacobi.c,v $
+ * Revision 1.1  1999/11/22 20:50:37  mdw
+ * Add support for computing Jacobi symbols.
+ *
+ */
+
+/*----- Header files ------------------------------------------------------*/
+
+#include "mp.h"
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- @mp_jacobi@ --- *
+ *
+ * Arguments:  @mp *a@ = an integer less than @n@
+ *             @mp *n@ = an odd integer
+ *
+ * Returns:    @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
+ *
+ * Use:                Computes the Jacobi symbol.  If @n@ is prime, this is the
+ *             Legendre symbol and is equal to 1 if and only if @a@ is a
+ *             quadratic residue mod @n@.  The result is zero if and only if
+ *             @a@ and @n@ have a common factor greater than one.
+ */
+
+int mp_jacobi(mp *a, mp *n)
+{
+  int s = 1;
+
+  /* --- Take copies of the arguments --- */
+
+  a = MP_COPY(a);
+  n = MP_COPY(n);
+
+  /* --- Main recursive mess, flattened out into something nice --- */
+
+  for (;;) {
+
+    /* --- Some simple special cases --- */
+
+    MP_SHRINK(a);
+
+    if (MP_LEN(a) == 0) {
+      s = 0;
+      goto done;
+    }
+
+    /* --- Find the power-of-two factor in @a@ --- */
+
+    {
+      mpscan sc;
+      mpw nn;
+      unsigned e;
+
+      /* --- Scan for a set bit --- */
+
+      MP_SCAN(&sc, a);
+      e = 0;
+      while (MP_STEP(&sc) && !MP_BIT(&sc))
+       e++;
+
+      /* --- Do the shift --- */
+
+      if (e)
+       a = mp_lsr(a, a, e);
+
+      /* --- Maybe adjust the sign of @s@ --- */
+
+      nn = n->v[0] & 7;
+      if ((e & 1) && (nn == 3 || nn == 5))
+       s = -s;
+
+      if (MP_LEN(a) == 1 && a->v[0] == 1)
+       goto done;
+
+      if ((nn & 3) == 3 && (a->v[0] & 3) == 3)
+       s = -s;
+    }
+
+    /* --- Reduce and swap --- */
+
+    mp_div(0, &n, n, a);
+    { mp *t = n; n = a; a = t; }
+  }
+
+  /* --- Wrap everything up --- */
+
+done:
+  MP_DROP(a);
+  MP_DROP(n);
+  return (s);
+}
+
+/*----- Test rig ----------------------------------------------------------*/
+
+#ifdef TEST_RIG
+
+#include <mLib/testrig.h>
+
+static int verify(dstr *v)
+{
+  mp *a = *(mp **)v[0].buf;
+  mp *n = *(mp **)v[1].buf;
+  int s = *(int *)v[2].buf;
+  int j = mp_jacobi(a, n);
+  int ok = 1;
+
+  if (s != j) {
+    fputs("\n*** fail", stderr);
+    fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
+    fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
+    fprintf(stderr, "s = %i\n", s);
+    fprintf(stderr, "j = %i\n", j);
+    ok = 0;
+  }
+
+  mp_drop(a);
+  mp_drop(n);
+  return (ok);
+}
+
+static test_chunk tests[] = {
+  { "jacobi", verify, { &type_mp, &type_mp, &type_int, 0 } },
+  { 0, 0, { 0 } }
+};
+
+int main(int argc, char *argv[])
+{
+  sub_init();
+  test_run(argc, argv, tests, SRCDIR "/tests/mp");
+  return (0);
+}
+
+#endif
+
+/*----- That's all, folks -------------------------------------------------*/

diff --git a/mp.h b/mp.h
index 13bcd3d..dce8128 100644 (file)
--- a/mp.h
+++ b/mp.h
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mp.h,v 1.4 1999/11/21 22:13:02 mdw Exp $
+ * $Id: mp.h,v 1.5 1999/11/22 20:50:37 mdw Exp $
  *
  * Simple multiprecision arithmetic
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mp.h,v $
+ * Revision 1.5  1999/11/22 20:50:37  mdw
+ * Add support for computing Jacobi symbols.
+ *
  * Revision 1.4  1999/11/21 22:13:02  mdw
  * Add mp version of MPX_BITS.
  *
@@ -639,6 +642,21 @@ extern void mp_div(mp **/*qq*/, mp **/*rr*/,
 extern void mp_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/,
                   mp */*a*/, mp */*b*/);
 
+/* --- @mp_jacobi@ --- *
+ *
+ * Arguments:  @mp *a@ = an integer less than @n@
+ *             @mp *n@ = an odd integer
+ *
+ * Returns:    @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
+ *
+ * Use:                Computes the Jacobi symbol.  If @n@ is prime, this is the
+ *             Legendre symbol and is equal to 1 if and only if @a@ is a
+ *             quadratic residue mod @n@.  The result is zero if and only if
+ *             @a@ and @n@ have a common factor greater than one.
+ */
+
+int mp_jacobi(mp */*a*/, mp */*n*/);
+
 /*----- Test harness support ----------------------------------------------*/
 
 #include <mLib/testrig.h>

diff --git a/tests/mp b/tests/mp
index 27854a2..910b89d 100644 (file)
--- a/tests/mp
+++ b/tests/mp
@@ -1,6 +1,6 @@
 # Test vectors for MP functions
 #
-# $Id: mp,v 1.1 1999/11/17 18:02:17 mdw Exp $
+# $Id: mp,v 1.2 1999/11/22 20:51:09 mdw Exp $
 
 add {
   5 4 9; 5 -4 1; -5 4 -1; -5 -4 -9;
@@ -55,3 +55,10 @@ gcd {
   -4601007896041464028712478963832994007038251361995647370
   514778499400157641662814932021958856708417966520837469125919104431;
 }
+
+jacobi {
+  4 5 1;
+  6 7 -1;
+  15 27 0;
+  2132498039840981 98729378979237498798347932749951 1;
+}