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

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