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

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