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

   1 /* -*-c-*-
   2  *
   3  * Compute Jacobi symbol
   4  *
   5  * (c) 1999 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include "mp.h"
  31
  32 /*----- Main code ---------------------------------------------------------*/
  33
  34 /* --- @mp_jacobi@ --- *
  35  *
  36  * Arguments:   @mp *a@ = an integer
  37  *              @mp *n@ = another integer
  38  *
  39  * Returns:     @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
  40  *
  41  * Use:         Computes the Kronecker symbol %$\jacobi{a}{n}$%.  If @n@ is
  42  *              prime, this is the Legendre symbol and is equal to 1 if and
  43  *              only if @a@ is a quadratic residue mod @n@.  The result is
  44  *              zero if and only if @a@ and @n@ have a common factor greater
  45  *              than one.
  46  *
  47  *              If @n@ is composite, then this computes the Kronecker symbol
  48  *
  49  *                %$\jacobi{a}{n}=\jacobi{a}{u}\prod_i\jacobi{a}{p_i}^{e_i}$%
  50  *
  51  *              where %$n = u p_0^{e_0} \ldots p_{n-1}^{e_{n-1}}$% is the
  52  *              prime factorization of %$n$%.  The missing bits are:
  53  *
  54  *                * %$\jacobi{a}{1} = 1$%;
  55  *                * %$\jacobi{a}{-1} = 1$% if @a@ is negative, or 1 if
  56  *                  positive;
  57  *                * %$\jacobi{a}{0} = 0$%;
  58  *                * %$\jacobi{a}{2}$ is 0 if @a@ is even, 1 if @a@ is
  59  *                  congruent to 1 or 7 (mod 8), or %$-1$% otherwise.
  60  *
  61  *              If %$n$% is positive and odd, then this is the Jacobi
  62  *              symbol.  (The Kronecker symbol is a consistent domain
  63  *              extension; the Jacobi symbol was implemented first, and the
  64  *              name stuck.)
  65  */
  66
  67 int mp_jacobi(mp *a, mp *n)
  68 {
  69   int s = 1;
  70   size_t p2;
  71
  72   /* --- Handle zero specially --- *
  73    *
  74    * I can't find any specific statement for what to do when %$n = 0$%; PARI
  75    * opts to set %$\jacobi{\pm1}{0} = \pm 1$% and %$\jacobi{a}{0} = 0$% for
  76    * other %$a$%.
  77    */
  78
  79   if (MP_ZEROP(n)) {
  80     if (MP_EQ(a, MP_ONE)) return (+1);
  81     else if (MP_EQ(a, MP_MONE)) return (-1);
  82     else return (0);
  83   }
  84
  85   /* --- Deal with powers of two --- *
  86    *
  87    * This implicitly takes a copy of %$n$%.  Copy %$a$% at the same time to
  88    * make cleanup easier.
  89    */
  90
  91   MP_COPY(a);
  92   n = mp_odd(MP_NEW, n, &p2);
  93   if (p2) {
  94     if (MP_EVENP(a)) {
  95       s = 0;
  96       goto done;
  97     } else if ((p2 & 1) && ((a->v[0] & 7) == 3 || (a->v[0] & 7) == 5))
  98       s = -s;
  99   }
 100
 101   /* --- Deal with negative %$n$% --- */
 102
 103   if (MP_NEGP(n)) {
 104     n = mp_neg(n, n);
 105     if (MP_NEGP(a))
 106       s = -s;
 107   }
 108
 109   /* --- Check for unit %$n$% --- */
 110
 111   if (MP_EQ(n, MP_ONE))
 112     goto done;
 113
 114   /* --- Reduce %$a$% modulo %$n$% --- */
 115
 116   if (MP_NEGP(a) || MP_CMP(a, >=, n))
 117     mp_div(0, &a, a, n);
 118
 119   /* --- Main recursive mess, flattened out into something nice --- */
 120
 121   for (;;) {
 122     mpw nn;
 123     size_t e;
 124
 125     /* --- Some simple special cases --- */
 126
 127     MP_SHRINK(a);
 128     if (MP_ZEROP(a)) {
 129       s = 0;
 130       goto done;
 131     }
 132
 133     /* --- Strip out powers of two from %$a$% --- */
 134
 135     a = mp_odd(a, a, &e);
 136     nn = n->v[0] & 7;
 137     if ((e & 1) && (nn == 3 || nn == 5))
 138       s = -s;
 139     if (MP_LEN(a) == 1 && a->v[0] == 1)
 140       goto done;
 141
 142     /* --- Reduce and swap, applying quadratic residuosity --- */
 143
 144     if ((nn & 3) == 3 && (a->v[0] & 3) == 3)
 145       s = -s;
 146     mp_div(0, &n, n, a);
 147     { mp *t = n; n = a; a = t; }
 148   }
 149
 150   /* --- Wrap everything up --- */
 151
 152 done:
 153   MP_DROP(a);
 154   MP_DROP(n);
 155   return (s);
 156 }
 157
 158 /*----- Test rig ----------------------------------------------------------*/
 159
 160 #ifdef TEST_RIG
 161
 162 #include <mLib/testrig.h>
 163
 164 static int verify(dstr *v)
 165 {
 166   mp *a = *(mp **)v[0].buf;
 167   mp *n = *(mp **)v[1].buf;
 168   int s = *(int *)v[2].buf;
 169   int j = mp_jacobi(a, n);
 170   int ok = 1;
 171
 172   if (s != j) {
 173     fputs("\n*** fail", stderr);
 174     fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 175     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 176     fprintf(stderr, "s = %i\n", s);
 177     fprintf(stderr, "j = %i\n", j);
 178     ok = 0;
 179   }
 180
 181   mp_drop(a);
 182   mp_drop(n);
 183   assert(mparena_count(MPARENA_GLOBAL) == 0);
 184   return (ok);
 185 }
 186
 187 static test_chunk tests[] = {
 188   { "jacobi", verify, { &type_mp, &type_mp, &type_int, 0 } },
 189   { 0, 0, { 0 } }
 190 };
 191
 192 int main(int argc, char *argv[])
 193 {
 194   sub_init();
 195   test_run(argc, argv, tests, SRCDIR "/t/mp");
 196   return (0);
 197 }
 198
 199 #endif
 200
 201 /*----- That's all, folks -------------------------------------------------*/