1 /* jacobi_2.c
2
3 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
4 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
5 GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
6
7 Copyright 1996, 1998, 2000-2004, 2008, 2010 Free Software Foundation, Inc.
8
9 This file is part of the GNU MP Library.
10
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of either:
13
14 * the GNU Lesser General Public License as published by the Free
15 Software Foundation; either version 3 of the License, or (at your
16 option) any later version.
17
18 or
19
20 * the GNU General Public License as published by the Free Software
21 Foundation; either version 2 of the License, or (at your option) any
22 later version.
23
24 or both in parallel, as here.
25
26 The GNU MP Library is distributed in the hope that it will be useful, but
27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
29 for more details.
30
31 You should have received copies of the GNU General Public License and the
32 GNU Lesser General Public License along with the GNU MP Library. If not,
33 see https://www.gnu.org/licenses/. */
34
35 #include "gmp-impl.h"
36 #include "longlong.h"
37
38 #ifndef JACOBI_2_METHOD
39 #define JACOBI_2_METHOD 2
40 #endif
41
42 /* Computes (a / b) where b is odd, and a and b are otherwise arbitrary
43 two-limb numbers. */
44 #if JACOBI_2_METHOD == 1
45 int
46 mpn_jacobi_2 (mp_srcptr ap, mp_srcptr bp, unsigned bit)
47 {
48 mp_limb_t ah, al, bh, bl;
49 int c;
50
51 al = ap[0];
52 ah = ap[1];
53 bl = bp[0];
54 bh = bp[1];
55
56 ASSERT (bl & 1);
57
58 bl = ((bh << (GMP_NUMB_BITS - 1)) & GMP_NUMB_MASK) | (bl >> 1);
59 bh >>= 1;
60
61 if ( (bh | bl) == 0)
62 return 1 - 2*(bit & 1);
63
64 if ( (ah | al) == 0)
65 return 0;
66
67 if (al == 0)
68 {
69 al = ah;
70 ah = 0;
71 bit ^= GMP_NUMB_BITS & (bl ^ (bl >> 1));
72 }
73 count_trailing_zeros (c, al);
74 bit ^= c & (bl ^ (bl >> 1));
75
76 c++;
77 if (UNLIKELY (c == GMP_NUMB_BITS))
78 {
79 al = ah;
80 ah = 0;
81 }
82 else
83 {
84 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
85 ah >>= c;
86 }
87 while ( (ah | bh) > 0)
88 {
89 mp_limb_t th, tl;
90 mp_limb_t bgta;
91
92 sub_ddmmss (th, tl, ah, al, bh, bl);
93 if ( (tl | th) == 0)
94 return 0;
95
96 bgta = LIMB_HIGHBIT_TO_MASK (th);
97
98 /* If b > a, invoke reciprocity */
99 bit ^= (bgta & al & bl);
100
101 /* b <-- min (a, b) */
102 add_ssaaaa (bh, bl, bh, bl, th & bgta, tl & bgta);
103
104 if ( (bh | bl) == 0)
105 return 1 - 2*(bit & 1);
106
107 /* a <-- |a - b| */
108 al = (bgta ^ tl) - bgta;
109 ah = (bgta ^ th);
110
111 if (UNLIKELY (al == 0))
112 {
113 /* If b > a, al == 0 implies that we have a carry to
114 propagate. */
115 al = ah - bgta;
116 ah = 0;
117 bit ^= GMP_NUMB_BITS & (bl ^ (bl >> 1));
118 }
119 count_trailing_zeros (c, al);
120 c++;
121 bit ^= c & (bl ^ (bl >> 1));
122
123 if (UNLIKELY (c == GMP_NUMB_BITS))
124 {
125 al = ah;
126 ah = 0;
127 }
128 else
129 {
130 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
131 ah >>= c;
132 }
133 }
134
135 ASSERT (bl > 0);
136
137 while ( (al | bl) & GMP_LIMB_HIGHBIT)
138 {
139 /* Need an extra comparison to get the mask. */
140 mp_limb_t t = al - bl;
141 mp_limb_t bgta = - (bl > al);
142
143 if (t == 0)
144 return 0;
145
146 /* If b > a, invoke reciprocity */
147 bit ^= (bgta & al & bl);
148
149 /* b <-- min (a, b) */
150 bl += (bgta & t);
151
152 /* a <-- |a - b| */
153 al = (t ^ bgta) - bgta;
154
155 /* Number of trailing zeros is the same no matter if we look at
156 * t or a, but using t gives more parallelism. */
157 count_trailing_zeros (c, t);
158 c ++;
159 /* (2/b) = -1 if b = 3 or 5 mod 8 */
160 bit ^= c & (bl ^ (bl >> 1));
161
162 if (UNLIKELY (c == GMP_NUMB_BITS))
163 return 1 - 2*(bit & 1);
164
165 al >>= c;
166 }
167
168 /* Here we have a little impedance mismatch. Better to inline it? */
169 return mpn_jacobi_base (2*al+1, 2*bl+1, bit << 1);
170 }
171 #elif JACOBI_2_METHOD == 2
172 int
173 mpn_jacobi_2 (mp_srcptr ap, mp_srcptr bp, unsigned bit)
174 {
175 mp_limb_t ah, al, bh, bl;
176 int c;
177
178 al = ap[0];
179 ah = ap[1];
180 bl = bp[0];
181 bh = bp[1];
182
183 ASSERT (bl & 1);
184
185 /* Use bit 1. */
186 bit <<= 1;
187
188 if (bh == 0 && bl == 1)
189 /* (a|1) = 1 */
190 return 1 - (bit & 2);
191
192 if (al == 0)
193 {
194 if (ah == 0)
195 /* (0|b) = 0, b > 1 */
196 return 0;
197
198 count_trailing_zeros (c, ah);
199 bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1));
200
201 al = bl;
202 bl = ah >> c;
203
204 if (bl == 1)
205 /* (1|b) = 1 */
206 return 1 - (bit & 2);
207
208 ah = bh;
209
210 bit ^= al & bl;
211
212 goto b_reduced;
213 }
214 if ( (al & 1) == 0)
215 {
216 count_trailing_zeros (c, al);
217
218 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
219 ah >>= c;
220 bit ^= (c << 1) & (bl ^ (bl >> 1));
221 }
222 if (ah == 0)
223 {
224 if (bh > 0)
225 {
226 bit ^= al & bl;
227 MP_LIMB_T_SWAP (al, bl);
228 ah = bh;
229 goto b_reduced;
230 }
231 goto ab_reduced;
232 }
233
234 while (bh > 0)
235 {
236 /* Compute (a|b) */
237 while (ah > bh)
238 {
239 sub_ddmmss (ah, al, ah, al, bh, bl);
240 if (al == 0)
241 {
242 count_trailing_zeros (c, ah);
243 bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1));
244
245 al = bl;
246 bl = ah >> c;
247 ah = bh;
248
249 bit ^= al & bl;
250 goto b_reduced;
251 }
252 count_trailing_zeros (c, al);
253 bit ^= (c << 1) & (bl ^ (bl >> 1));
254 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
255 ah >>= c;
256 }
257 if (ah == bh)
258 goto cancel_hi;
259
260 if (ah == 0)
261 {
262 bit ^= al & bl;
263 MP_LIMB_T_SWAP (al, bl);
264 ah = bh;
265 break;
266 }
267
268 bit ^= al & bl;
269
270 /* Compute (b|a) */
271 while (bh > ah)
272 {
273 sub_ddmmss (bh, bl, bh, bl, ah, al);
274 if (bl == 0)
275 {
276 count_trailing_zeros (c, bh);
277 bit ^= ((GMP_NUMB_BITS + c) << 1) & (al ^ (al >> 1));
278
279 bl = bh >> c;
280 bit ^= al & bl;
281 goto b_reduced;
282 }
283 count_trailing_zeros (c, bl);
284 bit ^= (c << 1) & (al ^ (al >> 1));
285 bl = ((bh << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (bl >> c);
286 bh >>= c;
287 }
288 bit ^= al & bl;
289
290 /* Compute (a|b) */
291 if (ah == bh)
292 {
293 cancel_hi:
294 if (al < bl)
295 {
296 MP_LIMB_T_SWAP (al, bl);
297 bit ^= al & bl;
298 }
299 al -= bl;
300 if (al == 0)
301 return 0;
302
303 count_trailing_zeros (c, al);
304 bit ^= (c << 1) & (bl ^ (bl >> 1));
305 al >>= c;
306
307 if (al == 1)
308 return 1 - (bit & 2);
309
310 MP_LIMB_T_SWAP (al, bl);
311 bit ^= al & bl;
312 break;
313 }
314 }
315
316 b_reduced:
317 /* Compute (a|b), with b a single limb. */
318 ASSERT (bl & 1);
319
320 if (bl == 1)
321 /* (a|1) = 1 */
322 return 1 - (bit & 2);
323
324 while (ah > 0)
325 {
326 ah -= (al < bl);
327 al -= bl;
328 if (al == 0)
329 {
330 if (ah == 0)
331 return 0;
332 count_trailing_zeros (c, ah);
333 bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1));
334 al = ah >> c;
335 goto ab_reduced;
336 }
337 count_trailing_zeros (c, al);
338
339 al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
340 ah >>= c;
341 bit ^= (c << 1) & (bl ^ (bl >> 1));
342 }
343 ab_reduced:
344 ASSERT (bl & 1);
345 ASSERT (bl > 1);
346
347 return mpn_jacobi_base (al, bl, bit);
348 }
349 #else
350 #error Unsupported value for JACOBI_2_METHOD
351 #endif