Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpz/jacobi.c
Line
Count
Source (jump to first uncovered line)
1
/* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
2
3
Copyright 2000-2002, 2005, 2010-2012 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
8
it under the terms of either:
9
10
  * the GNU Lesser General Public License as published by the Free
11
    Software Foundation; either version 3 of the License, or (at your
12
    option) any later version.
13
14
or
15
16
  * the GNU General Public License as published by the Free Software
17
    Foundation; either version 2 of the License, or (at your option) any
18
    later version.
19
20
or both in parallel, as here.
21
22
The GNU MP Library is distributed in the hope that it will be useful, but
23
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
25
for more details.
26
27
You should have received copies of the GNU General Public License and the
28
GNU Lesser General Public License along with the GNU MP Library.  If not,
29
see https://www.gnu.org/licenses/.  */
30
31
#include <stdio.h>
32
#include "gmp-impl.h"
33
#include "longlong.h"
34
35
36
/* This code does triple duty as mpz_jacobi, mpz_legendre and
37
   mpz_kronecker. For ABI compatibility, the link symbol is
38
   __gmpz_jacobi, not __gmpz_kronecker, even though the latter would
39
   be more logical.
40
41
   mpz_jacobi could assume b is odd, but the improvements from that seem
42
   small compared to other operations, and anything significant should be
43
   checked at run-time since we'd like odd b to go fast in mpz_kronecker
44
   too.
45
46
   mpz_legendre could assume b is an odd prime, but knowing this doesn't
47
   present any obvious benefits.  Result 0 wouldn't arise (unless "a" is a
48
   multiple of b), but the checking for that takes little time compared to
49
   other operations.
50
51
   Enhancements:
52
53
   mpn_bdiv_qr should be used instead of mpn_tdiv_qr.
54
55
*/
56
57
int
58
mpz_jacobi (mpz_srcptr a, mpz_srcptr b)
59
720
{
60
720
  mp_srcptr  asrcp, bsrcp;
61
720
  mp_size_t  asize, bsize;
62
720
  mp_limb_t  alow, blow;
63
720
  mp_ptr     ap, bp;
64
720
  unsigned   btwos;
65
720
  int        result_bit1;
66
720
  int        res;
67
720
  TMP_DECL;
68
69
720
  asize = SIZ(a);
70
720
  asrcp = PTR(a);
71
720
  alow = asrcp[0];
72
73
720
  bsize = SIZ(b);
74
720
  bsrcp = PTR(b);
75
720
  blow = bsrcp[0];
76
77
  /* The MPN jacobi functions require positive a and b, and b odd. So
78
     we must to handle the cases of a or b zero, then signs, and then
79
     the case of even b.
80
  */
81
82
720
  if (bsize == 0)
83
    /* (a/0) = [ a = 1 or a = -1 ] */
84
9
    return JACOBI_LS0 (alow, asize);
85
86
711
  if (asize == 0)
87
    /* (0/b) = [ b = 1 or b = - 1 ] */
88
5
    return JACOBI_0LS (blow, bsize);
89
90
706
  if ( (((alow | blow) & 1) == 0))
91
    /* Common factor of 2 ==> (a/b) = 0 */
92
8
    return 0;
93
94
698
  if (bsize < 0)
95
0
    {
96
      /* (a/-1) = -1 if a < 0, +1 if a >= 0 */
97
0
      result_bit1 = (asize < 0) << 1;
98
0
      bsize = -bsize;
99
0
    }
100
698
  else
101
698
    result_bit1 = 0;
102
103
698
  JACOBI_STRIP_LOW_ZEROS (result_bit1, alow, bsrcp, bsize, blow);
104
105
698
  count_trailing_zeros (btwos, blow);
106
698
  blow >>= btwos;
107
108
698
  if (bsize > 1 && btwos > 0)
109
111
    {
110
111
      mp_limb_t b1 = bsrcp[1];
111
111
      blow |= b1 << (GMP_NUMB_BITS - btwos);
112
111
      if (bsize == 2 && (b1 >> btwos) == 0)
113
0
  bsize = 1;
114
111
    }
115
116
698
  if (asize < 0)
117
0
    {
118
      /* (-1/b) = -1 iff b = 3 (mod 4) */
119
0
      result_bit1 ^= JACOBI_N1B_BIT1(blow);
120
0
      asize = -asize;
121
0
    }
122
123
698
  JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, asrcp, asize, alow);
124
125
  /* Ensure asize >= bsize. Take advantage of the generalized
126
     reciprocity law (a/b*2^n) = (b*2^n / a) * RECIP(a,b) */
127
128
698
  if (asize < bsize)
129
237
    {
130
237
      MPN_SRCPTR_SWAP (asrcp, asize, bsrcp, bsize);
131
237
      MP_LIMB_T_SWAP (alow, blow);
132
133
      /* NOTE: The value of alow (old blow) is a bit subtle. For this code
134
   path, we get alow as the low, always odd, limb of shifted A. Which is
135
   what we need for the reciprocity update below.
136
137
   However, all other uses of alow assumes that it is *not*
138
   shifted. Luckily, alow matters only when either
139
140
   + btwos > 0, in which case A is always odd
141
142
   + asize == bsize == 1, in which case this code path is never
143
     taken. */
144
145
237
      count_trailing_zeros (btwos, blow);
146
237
      blow >>= btwos;
147
148
237
      if (bsize > 1 && btwos > 0)
149
124
  {
150
124
    mp_limb_t b1 = bsrcp[1];
151
124
    blow |= b1 << (GMP_NUMB_BITS - btwos);
152
124
    if (bsize == 2 && (b1 >> btwos) == 0)
153
2
      bsize = 1;
154
124
  }
155
156
237
      result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
157
237
    }
158
159
698
  if (bsize == 1)
160
56
    {
161
56
      result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow);
162
163
56
      if (blow == 1)
164
4
  return JACOBI_BIT1_TO_PN (result_bit1);
165
166
52
      if (asize > 1)
167
15
  JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow);
168
169
52
      return mpn_jacobi_base (alow, blow, result_bit1);
170
52
    }
171
172
  /* Allocation strategy: For A, we allocate a working copy only for A % B, but
173
     when A is much larger than B, we have to allocate space for the large
174
     quotient. We use the same area, pointed to by bp, for both the quotient
175
     A/B and the working copy of B. */
176
177
642
  TMP_MARK;
178
179
642
  if (asize >= 2*bsize)
180
66
    TMP_ALLOC_LIMBS_2 (ap, bsize, bp, asize - bsize + 1);
181
576
  else
182
576
    TMP_ALLOC_LIMBS_2 (ap, bsize, bp, bsize);
183
184
  /* In the case of even B, we conceptually shift out the powers of two first,
185
     and then divide A mod B. Hence, when taking those powers of two into
186
     account, we must use alow *before* the division. Doing the actual division
187
     first is ok, because the point is to remove multiples of B from A, and
188
     multiples of 2^k B are good enough. */
189
642
  if (asize > bsize)
190
455
    mpn_tdiv_qr (bp, ap, 0, asrcp, asize, bsrcp, bsize);
191
187
  else
192
187
    MPN_COPY (ap, asrcp, bsize);
193
194
642
  if (btwos > 0)
195
213
    {
196
213
      result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow);
197
198
213
      ASSERT_NOCARRY (mpn_rshift (bp, bsrcp, bsize, btwos));
199
213
      bsize -= (ap[bsize-1] | bp[bsize-1]) == 0;
200
213
    }
201
429
  else
202
429
    MPN_COPY (bp, bsrcp, bsize);
203
204
642
  ASSERT (blow == bp[0]);
205
642
  res = mpn_jacobi_n (ap, bp, bsize,
206
642
          mpn_jacobi_init (ap[0], blow, (result_bit1>>1) & 1));
207
208
642
  TMP_FREE;
209
642
  return res;
210
642
}