contrib/gmp/mpz/jacobi.c

   1 /* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
   2
   3 Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify it
   8 under the terms of the GNU Lesser General Public License as published by the
   9 Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License
  15 for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License along
  18 with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include <stdio.h>
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23 #include "longlong.h"
  24
  25
  26 /* Change this to "#define TRACE(x) x" for some traces. */
  27 #define TRACE(x)
  28
  29
  30 #define MPN_RSHIFT_OR_COPY(dst,src,size,shift)                  \
  31   do {                                                          \
  32     if ((shift) != 0)                                           \
  33       {                                                         \
  34         ASSERT_NOCARRY (mpn_rshift (dst, src, size, shift));    \
  35         (size) -= ((dst)[(size)-1] == 0);                       \
  36       }                                                         \
  37     else                                                        \
  38       MPN_COPY (dst, src, size);                                \
  39   } while (0)
  40
  41
  42 /* This code does triple duty as mpz_jacobi, mpz_legendre and mpz_kronecker.
  43
  44    mpz_jacobi could assume b is odd, but the improvements from that seem
  45    small compared to other operations, and anything significant should be
  46    checked at run-time since we'd like odd b to go fast in mpz_kronecker
  47    too.
  48
  49    mpz_legendre could assume b is an odd prime, but knowing this doesn't
  50    present any obvious benefits.  Result 0 wouldn't arise (unless "a" is a
  51    multiple of b), but the checking for that takes little time compared to
  52    other operations.
  53
  54    The main loop is just a simple binary GCD with the jacobi symbol result
  55    tracked during the reduction.
  56
  57    The special cases for a or b fitting in one limb let mod_1 or modexact_1
  58    get used, without any copying, and end up just as efficient as the mixed
  59    precision mpz_kronecker_ui etc.
  60
  61    When tdiv_qr is called it's not necessary to make "a" odd or make a
  62    working copy of it, but tdiv_qr is going to be pretty slow so it's not
  63    worth bothering trying to save anything for that case.
  64
  65    Enhancements:
  66
  67    mpn_bdiv_qr should be used instead of mpn_tdiv_qr.
  68
  69    Some sort of multi-step algorithm should be used.  The current subtract
  70    and shift for every bit is very inefficient.  Lehmer (per current gcdext)
  71    would need some low bits included in its calculation to apply the sign
  72    change for reciprocity.  Binary Lehmer keeps low bits to strip twos
  73    anyway, so might be better suited.  Maybe the accelerated GCD style k-ary
  74    reduction would work, if sign changes due to the extra factors it
  75    introduces can be accounted for (or maybe they can be ignored).  */
  76
  77
  78 int
  79 mpz_jacobi (mpz_srcptr a, mpz_srcptr b)
  80 {
  81   mp_srcptr  asrcp, bsrcp;
  82   mp_size_t  asize, bsize;
  83   mp_ptr     ap, bp;
  84   mp_limb_t  alow, blow, ahigh, bhigh, asecond, bsecond;
  85   unsigned   atwos, btwos;
  86   int        result_bit1;
  87   TMP_DECL;
  88
  89   TRACE (printf ("start asize=%d bsize=%d\n", SIZ(a), SIZ(b));
  90          mpz_trace (" a", a);
  91          mpz_trace (" b", b));
  92
  93   asize = SIZ(a);
  94   asrcp = PTR(a);
  95   alow = asrcp[0];
  96
  97   bsize = SIZ(b);
  98   if (bsize == 0)
  99     return JACOBI_LS0 (alow, asize);  /* (a/0) */
 100
 101   bsrcp = PTR(b);
 102   blow = bsrcp[0];
 103
 104   if (asize == 0)
 105     return JACOBI_0LS (blow, bsize);  /* (0/b) */
 106
 107   /* (even/even)=0 */
 108   if (((alow | blow) & 1) == 0)
 109     return 0;
 110
 111   /* account for effect of sign of b, then ignore it */
 112   result_bit1 = JACOBI_BSGN_SS_BIT1 (asize, bsize);
 113   bsize = ABS (bsize);
 114
 115   /* low zero limbs on b can be discarded */
 116   JACOBI_STRIP_LOW_ZEROS (result_bit1, alow, bsrcp, bsize, blow);
 117
 118   count_trailing_zeros (btwos, blow);
 119   TRACE (printf ("b twos %u\n", btwos));
 120
 121   /* establish shifted blow */
 122   blow >>= btwos;
 123   if (bsize > 1)
 124     {
 125       bsecond = bsrcp[1];
 126       if (btwos != 0)
 127         blow |= (bsecond << (GMP_NUMB_BITS - btwos)) & GMP_NUMB_MASK;
 128     }
 129
 130   /* account for effect of sign of a, then ignore it */
 131   result_bit1 ^= JACOBI_ASGN_SU_BIT1 (asize, blow);
 132   asize = ABS (asize);
 133
 134   if (bsize == 1 || (bsize == 2 && (bsecond >> btwos) == 0))
 135     {
 136       /* special case one limb b, use modexact and no copying */
 137
 138       /* (a/2)=(2/a) with a odd, and if b is even then a is odd here */
 139       result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
 140
 141       if (blow == 1)   /* (a/1)=1 always */
 142         return JACOBI_BIT1_TO_PN (result_bit1);
 143
 144       JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow);
 145       TRACE (printf ("base (%lu/%lu) with %d\n",
 146                      alow, blow, JACOBI_BIT1_TO_PN (result_bit1)));
 147       return mpn_jacobi_base (alow, blow, result_bit1);
 148     }
 149
 150   /* Discard low zero limbs of a.  Usually there won't be anything to
 151      strip, hence not bothering with it for the bsize==1 case.  */
 152   JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, asrcp, asize, alow);
 153
 154   count_trailing_zeros (atwos, alow);
 155   TRACE (printf ("a twos %u\n", atwos));
 156   result_bit1 ^= JACOBI_TWOS_U_BIT1 (atwos, blow);
 157
 158   /* establish shifted alow */
 159   alow >>= atwos;
 160   if (asize > 1)
 161     {
 162       asecond = asrcp[1];
 163       if (atwos != 0)
 164         alow |= (asecond << (GMP_NUMB_BITS - atwos)) & GMP_NUMB_MASK;
 165     }
 166
 167   /* (a/2)=(2/a) with a odd */
 168   result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
 169
 170   if (asize == 1 || (asize == 2 && (asecond >> atwos) == 0))
 171     {
 172       /* another special case with modexact and no copying */
 173
 174       if (alow == 1)  /* (1/b)=1 always */
 175         return JACOBI_BIT1_TO_PN (result_bit1);
 176
 177       /* b still has its twos, so cancel out their effect */
 178       result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
 179
 180       result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);  /* now (b/a) */
 181       JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, blow, bsrcp, bsize, alow);
 182       TRACE (printf ("base (%lu/%lu) with %d\n",
 183                      blow, alow, JACOBI_BIT1_TO_PN (result_bit1)));
 184       return mpn_jacobi_base (blow, alow, result_bit1);
 185     }
 186
 187
 188   TMP_MARK;
 189   TMP_ALLOC_LIMBS_2 (ap, asize, bp, bsize);
 190
 191   MPN_RSHIFT_OR_COPY (ap, asrcp, asize, atwos);
 192   ASSERT (alow == ap[0]);
 193   TRACE (mpn_trace ("stripped a", ap, asize));
 194
 195   MPN_RSHIFT_OR_COPY (bp, bsrcp, bsize, btwos);
 196   ASSERT (blow == bp[0]);
 197   TRACE (mpn_trace ("stripped b", bp, bsize));
 198
 199   /* swap if necessary to make a longer than b */
 200   if (asize < bsize)
 201     {
 202       TRACE (printf ("swap\n"));
 203       MPN_PTR_SWAP (ap,asize, bp,bsize);
 204       MP_LIMB_T_SWAP (alow, blow);
 205       result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
 206     }
 207
 208   /* If a is bigger than b then reduce to a mod b.
 209      Division is much faster than chipping away at "a" bit-by-bit. */
 210   if (asize > bsize)
 211     {
 212       mp_ptr  rp, qp;
 213
 214       TRACE (printf ("tdiv_qr asize=%ld bsize=%ld\n", asize, bsize));
 215
 216       TMP_ALLOC_LIMBS_2 (rp, bsize, qp, asize-bsize+1);
 217       mpn_tdiv_qr (qp, rp, (mp_size_t) 0, ap, asize, bp, bsize);
 218       ap = rp;
 219       asize = bsize;
 220       MPN_NORMALIZE (ap, asize);
 221
 222       TRACE (printf ("tdiv_qr asize=%ld bsize=%ld\n", asize, bsize);
 223              mpn_trace (" a", ap, asize);
 224              mpn_trace (" b", bp, bsize));
 225
 226       if (asize == 0)  /* (0/b)=0 for b!=1 */
 227         goto zero;
 228
 229       alow = ap[0];
 230       goto strip_a;
 231     }
 232
 233   for (;;)
 234     {
 235       ASSERT (asize >= 1);         /* a,b non-empty */
 236       ASSERT (bsize >= 1);
 237       ASSERT (ap[asize-1] != 0);   /* a,b normalized (and hence non-zero) */
 238       ASSERT (bp[bsize-1] != 0);
 239       ASSERT (alow == ap[0]);      /* low limb copies should be correct */
 240       ASSERT (blow == bp[0]);
 241       ASSERT (alow & 1);           /* a,b odd */
 242       ASSERT (blow & 1);
 243
 244       TRACE (printf ("top asize=%ld bsize=%ld\n", asize, bsize);
 245              mpn_trace (" a", ap, asize);
 246              mpn_trace (" b", bp, bsize));
 247
 248       /* swap if necessary to make a>=b, applying reciprocity
 249          high limbs are almost always enough to tell which is bigger */
 250       if (asize < bsize
 251           || (asize == bsize
 252               && ((ahigh=ap[asize-1]) < (bhigh=bp[asize-1])
 253                   || (ahigh == bhigh
 254                       && mpn_cmp (ap, bp, asize-1) < 0))))
 255         {
 256           TRACE (printf ("swap\n"));
 257           MPN_PTR_SWAP (ap,asize, bp,bsize);
 258           MP_LIMB_T_SWAP (alow, blow);
 259           result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
 260         }
 261
 262       if (asize == 1)
 263         break;
 264
 265       /* a = a-b */
 266       ASSERT (asize >= bsize);
 267       ASSERT_NOCARRY (mpn_sub (ap, ap, asize, bp, bsize));
 268       MPN_NORMALIZE (ap, asize);
 269       alow = ap[0];
 270
 271       /* (0/b)=0 for b!=1.  b!=1 when a==0 because otherwise would have had
 272          a==1 which is asize==1 and would have exited above.  */
 273       if (asize == 0)
 274         goto zero;
 275
 276     strip_a:
 277       /* low zero limbs on a can be discarded */
 278       JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, ap, asize, alow);
 279
 280       if ((alow & 1) == 0)
 281         {
 282           /* factors of 2 from a */
 283           unsigned  twos;
 284           count_trailing_zeros (twos, alow);
 285           TRACE (printf ("twos %u\n", twos));
 286           result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, blow);
 287           ASSERT_NOCARRY (mpn_rshift (ap, ap, asize, twos));
 288           asize -= (ap[asize-1] == 0);
 289           alow = ap[0];
 290         }
 291     }
 292
 293   ASSERT (asize == 1 && bsize == 1);  /* just alow and blow left */
 294   TMP_FREE;
 295
 296   /* (1/b)=1 always (in this case have b==1 because a>=b) */
 297   if (alow == 1)
 298     return JACOBI_BIT1_TO_PN (result_bit1);
 299
 300   /* swap with reciprocity and do (b/a) */
 301   result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
 302   TRACE (printf ("base (%lu/%lu) with %d\n",
 303                  blow, alow, JACOBI_BIT1_TO_PN (result_bit1)));
 304   return mpn_jacobi_base (blow, alow, result_bit1);
 305
 306  zero:
 307   TMP_FREE;
 308   return 0;
 309 }