chiark - git - mdw - catacomb/blob - mp-jacobi.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-jacobi.c,v 1.2 1999/12/10 23:19:02 mdw Exp $
   4  *
   5  * Compute Jacobi symbol
   6  *
   7  * (c) 1999 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mp-jacobi.c,v $
  33  * Revision 1.2  1999/12/10 23:19:02  mdw
  34  * Improve error-checking.
  35  *
  36  * Revision 1.1  1999/11/22 20:50:37  mdw
  37  * Add support for computing Jacobi symbols.
  38  *
  39  */
  40
  41 /*----- Header files ------------------------------------------------------*/
  42
  43 #include "mp.h"
  44
  45 /*----- Main code ---------------------------------------------------------*/
  46
  47 /* --- @mp_jacobi@ --- *
  48  *
  49  * Arguments:   @mp *a@ = an integer less than @n@
  50  *              @mp *n@ = an odd integer
  51  *
  52  * Returns:     @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
  53  *
  54  * Use:         Computes the Jacobi symbol.  If @n@ is prime, this is the
  55  *              Legendre symbol and is equal to 1 if and only if @a@ is a
  56  *              quadratic residue mod @n@.  The result is zero if and only if
  57  *              @a@ and @n@ have a common factor greater than one.
  58  */
  59
  60 int mp_jacobi(mp *a, mp *n)
  61 {
  62   int s = 1;
  63
  64   /* --- Take copies of the arguments --- */
  65
  66   a = MP_COPY(a);
  67   n = MP_COPY(n);
  68
  69   /* --- Main recursive mess, flattened out into something nice --- */
  70
  71   for (;;) {
  72
  73     /* --- Some simple special cases --- */
  74
  75     MP_SHRINK(a);
  76
  77     if (MP_LEN(a) == 0) {
  78       s = 0;
  79       goto done;
  80     }
  81
  82     /* --- Find the power-of-two factor in @a@ --- */
  83
  84     {
  85       mpscan sc;
  86       mpw nn;
  87       unsigned e;
  88
  89       /* --- Scan for a set bit --- */
  90
  91       MP_SCAN(&sc, a);
  92       e = 0;
  93       while (MP_STEP(&sc) && !MP_BIT(&sc))
  94         e++;
  95
  96       /* --- Do the shift --- */
  97
  98       if (e)
  99         a = mp_lsr(a, a, e);
 100
 101       /* --- Maybe adjust the sign of @s@ --- */
 102
 103       nn = n->v[0] & 7;
 104       if ((e & 1) && (nn == 3 || nn == 5))
 105         s = -s;
 106
 107       if (MP_LEN(a) == 1 && a->v[0] == 1)
 108         goto done;
 109
 110       if ((nn & 3) == 3 && (a->v[0] & 3) == 3)
 111         s = -s;
 112     }
 113
 114     /* --- Reduce and swap --- */
 115
 116     mp_div(0, &n, n, a);
 117     { mp *t = n; n = a; a = t; }
 118   }
 119
 120   /* --- Wrap everything up --- */
 121
 122 done:
 123   MP_DROP(a);
 124   MP_DROP(n);
 125   return (s);
 126 }
 127
 128 /*----- Test rig ----------------------------------------------------------*/
 129
 130 #ifdef TEST_RIG
 131
 132 #include <mLib/testrig.h>
 133
 134 static int verify(dstr *v)
 135 {
 136   mp *a = *(mp **)v[0].buf;
 137   mp *n = *(mp **)v[1].buf;
 138   int s = *(int *)v[2].buf;
 139   int j = mp_jacobi(a, n);
 140   int ok = 1;
 141
 142   if (s != j) {
 143     fputs("\n*** fail", stderr);
 144     fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 145     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 146     fprintf(stderr, "s = %i\n", s);
 147     fprintf(stderr, "j = %i\n", j);
 148     ok = 0;
 149   }
 150
 151   mp_drop(a);
 152   mp_drop(n);
 153   assert(mparena_count(MPARENA_GLOBAL) == 0);
 154   return (ok);
 155 }
 156
 157 static test_chunk tests[] = {
 158   { "jacobi", verify, { &type_mp, &type_mp, &type_int, 0 } },
 159   { 0, 0, { 0 } }
 160 };
 161
 162 int main(int argc, char *argv[])
 163 {
 164   sub_init();
 165   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 166   return (0);
 167 }
 168
 169 #endif
 170
 171 /*----- That's all, folks -------------------------------------------------*/