1/* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
2
3Copyright 2000, 2001, 2002, 2005 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 it
8under the terms of the GNU Lesser General Public License as published by the
9Free Software Foundation; either version 3 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 or
14FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License
15for more details.
16
17You should have received a copy of the GNU Lesser General Public License along
18with 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
78int
79mpz_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}
310