git.distorted.org.uk Git - u/mdw/catacomb/blob - mp-jacobi.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-jacobi.c,v 1.1 1999/11/22 20:50:37 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.1  1999/11/22 20:50:37  mdw
  34  * Add support for computing Jacobi symbols.
  35  *
  36  */
  37
  38 /*----- Header files ------------------------------------------------------*/
  39
  40 #include "mp.h"
  41
  42 /*----- Main code ---------------------------------------------------------*/
  43
  44 /* --- @mp_jacobi@ --- *
  45  *
  46  * Arguments:   @mp *a@ = an integer less than @n@
  47  *              @mp *n@ = an odd integer
  48  *
  49  * Returns:     @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
  50  *
  51  * Use:         Computes the Jacobi symbol.  If @n@ is prime, this is the
  52  *              Legendre symbol and is equal to 1 if and only if @a@ is a
  53  *              quadratic residue mod @n@.  The result is zero if and only if
  54  *              @a@ and @n@ have a common factor greater than one.
  55  */
  56
  57 int mp_jacobi(mp *a, mp *n)
  58 {
  59   int s = 1;
  60
  61   /* --- Take copies of the arguments --- */
  62
  63   a = MP_COPY(a);
  64   n = MP_COPY(n);
  65
  66   /* --- Main recursive mess, flattened out into something nice --- */
  67
  68   for (;;) {
  69
  70     /* --- Some simple special cases --- */
  71
  72     MP_SHRINK(a);
  73
  74     if (MP_LEN(a) == 0) {
  75       s = 0;
  76       goto done;
  77     }
  78
  79     /* --- Find the power-of-two factor in @a@ --- */
  80
  81     {
  82       mpscan sc;
  83       mpw nn;
  84       unsigned e;
  85
  86       /* --- Scan for a set bit --- */
  87
  88       MP_SCAN(&sc, a);
  89       e = 0;
  90       while (MP_STEP(&sc) && !MP_BIT(&sc))
  91         e++;
  92
  93       /* --- Do the shift --- */
  94
  95       if (e)
  96         a = mp_lsr(a, a, e);
  97
  98       /* --- Maybe adjust the sign of @s@ --- */
  99
 100       nn = n->v[0] & 7;
 101       if ((e & 1) && (nn == 3 || nn == 5))
 102         s = -s;
 103
 104       if (MP_LEN(a) == 1 && a->v[0] == 1)
 105         goto done;
 106
 107       if ((nn & 3) == 3 && (a->v[0] & 3) == 3)
 108         s = -s;
 109     }
 110
 111     /* --- Reduce and swap --- */
 112
 113     mp_div(0, &n, n, a);
 114     { mp *t = n; n = a; a = t; }
 115   }
 116
 117   /* --- Wrap everything up --- */
 118
 119 done:
 120   MP_DROP(a);
 121   MP_DROP(n);
 122   return (s);
 123 }
 124
 125 /*----- Test rig ----------------------------------------------------------*/
 126
 127 #ifdef TEST_RIG
 128
 129 #include <mLib/testrig.h>
 130
 131 static int verify(dstr *v)
 132 {
 133   mp *a = *(mp **)v[0].buf;
 134   mp *n = *(mp **)v[1].buf;
 135   int s = *(int *)v[2].buf;
 136   int j = mp_jacobi(a, n);
 137   int ok = 1;
 138
 139   if (s != j) {
 140     fputs("\n*** fail", stderr);
 141     fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 142     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 143     fprintf(stderr, "s = %i\n", s);
 144     fprintf(stderr, "j = %i\n", j);
 145     ok = 0;
 146   }
 147
 148   mp_drop(a);
 149   mp_drop(n);
 150   return (ok);
 151 }
 152
 153 static test_chunk tests[] = {
 154   { "jacobi", verify, { &type_mp, &type_mp, &type_int, 0 } },
 155   { 0, 0, { 0 } }
 156 };
 157
 158 int main(int argc, char *argv[])
 159 {
 160   sub_init();
 161   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 162   return (0);
 163 }
 164
 165 #endif
 166
 167 /*----- That's all, folks -------------------------------------------------*/