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

   1 /* -*-c-*-
   2  *
   3  * $Id$
   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 /*----- Header files ------------------------------------------------------*/
  31
  32 #include "mp.h"
  33
  34 /*----- Main code ---------------------------------------------------------*/
  35
  36 /* --- @mp_jacobi@ --- *
  37  *
  38  * Arguments:   @mp *a@ = an integer less than @n@
  39  *              @mp *n@ = an odd integer
  40  *
  41  * Returns:     @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
  42  *
  43  * Use:         Computes the Jacobi symbol.  If @n@ is prime, this is the
  44  *              Legendre symbol and is equal to 1 if and only if @a@ is a
  45  *              quadratic residue mod @n@.  The result is zero if and only if
  46  *              @a@ and @n@ have a common factor greater than one.
  47  */
  48
  49 int mp_jacobi(mp *a, mp *n)
  50 {
  51   int s = 1;
  52
  53   assert(MP_ODDP(n));
  54
  55   /* --- Take copies of the arguments --- */
  56
  57   a = MP_COPY(a);
  58   n = MP_COPY(n);
  59
  60   /* --- Main recursive mess, flattened out into something nice --- */
  61
  62   for (;;) {
  63     mpw nn;
  64     size_t e;
  65
  66     /* --- Some simple special cases --- */
  67
  68     MP_SHRINK(a);
  69     if (MP_ZEROP(a)) {
  70       s = 0;
  71       goto done;
  72     }
  73
  74     /* --- Main case with powers of two --- */
  75
  76     a = mp_odd(a, a, &e);
  77     nn = n->v[0] & 7;
  78     if ((e & 1) && (nn == 3 || nn == 5))
  79       s = -s;
  80     if (MP_LEN(a) == 1 && a->v[0] == 1)
  81       goto done;
  82     if ((nn & 3) == 3 && (a->v[0] & 3) == 3)
  83       s = -s;
  84
  85     /* --- Reduce and swap --- */
  86
  87     mp_div(0, &n, n, a);
  88     { mp *t = n; n = a; a = t; }
  89   }
  90
  91   /* --- Wrap everything up --- */
  92
  93 done:
  94   MP_DROP(a);
  95   MP_DROP(n);
  96   return (s);
  97 }
  98
  99 /*----- Test rig ----------------------------------------------------------*/
 100
 101 #ifdef TEST_RIG
 102
 103 #include <mLib/testrig.h>
 104
 105 static int verify(dstr *v)
 106 {
 107   mp *a = *(mp **)v[0].buf;
 108   mp *n = *(mp **)v[1].buf;
 109   int s = *(int *)v[2].buf;
 110   int j = mp_jacobi(a, n);
 111   int ok = 1;
 112
 113   if (s != j) {
 114     fputs("\n*** fail", stderr);
 115     fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 116     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 117     fprintf(stderr, "s = %i\n", s);
 118     fprintf(stderr, "j = %i\n", j);
 119     ok = 0;
 120   }
 121
 122   mp_drop(a);
 123   mp_drop(n);
 124   assert(mparena_count(MPARENA_GLOBAL) == 0);
 125   return (ok);
 126 }
 127
 128 static test_chunk tests[] = {
 129   { "jacobi", verify, { &type_mp, &type_mp, &type_int, 0 } },
 130   { 0, 0, { 0 } }
 131 };
 132
 133 int main(int argc, char *argv[])
 134 {
 135   sub_init();
 136   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 137   return (0);
 138 }
 139
 140 #endif
 141
 142 /*----- That's all, folks -------------------------------------------------*/