source: trunk/third/gmp/mpz/jacobi.c @ 18191

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