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.
7 Copyright 1996, 1998, 2000-2004, 2008, 2010 Free Software Foundation, Inc.
9 This file is part of the GNU MP Library.
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of either:
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.
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
24 or both in parallel, as here.
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
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/. */
39 #ifndef JACOBI_2_METHOD
40 #define JACOBI_2_METHOD 2
43 /* Computes (a / b) where b is odd, and a and b are otherwise arbitrary
45 #if JACOBI_2_METHOD == 1
47 mpn_jacobi_2 (mp_srcptr ap
, mp_srcptr bp
, unsigned bit
)
49 mp_limb_t ah
, al
, bh
, bl
;
59 bl
= ((bh
<< (GMP_NUMB_BITS
- 1)) & GMP_NUMB_MASK
) | (bl
>> 1);
63 return 1 - 2*(bit
& 1);
72 bit
^= GMP_NUMB_BITS
& (bl
^ (bl
>> 1));
74 count_trailing_zeros (c
, al
);
75 bit
^= c
& (bl
^ (bl
>> 1));
78 if (UNLIKELY (c
== GMP_NUMB_BITS
))
85 al
= ((ah
<< (GMP_NUMB_BITS
- c
)) & GMP_NUMB_MASK
) | (al
>> c
);
88 while ( (ah
| bh
) > 0)
93 sub_ddmmss (th
, tl
, ah
, al
, bh
, bl
);
97 bgta
= LIMB_HIGHBIT_TO_MASK (th
);
99 /* If b > a, invoke reciprocity */
100 bit
^= (bgta
& al
& bl
);
102 /* b <-- min (a, b) */
103 add_ssaaaa (bh
, bl
, bh
, bl
, th
& bgta
, tl
& bgta
);
106 return 1 - 2*(bit
& 1);
109 al
= (bgta
^ tl
) - bgta
;
112 if (UNLIKELY (al
== 0))
114 /* If b > a, al == 0 implies that we have a carry to
118 bit
^= GMP_NUMB_BITS
& (bl
^ (bl
>> 1));
120 count_trailing_zeros (c
, al
);
122 bit
^= c
& (bl
^ (bl
>> 1));
124 if (UNLIKELY (c
== GMP_NUMB_BITS
))
131 al
= ((ah
<< (GMP_NUMB_BITS
- c
)) & GMP_NUMB_MASK
) | (al
>> c
);
138 while ( (al
| bl
) & GMP_LIMB_HIGHBIT
)
140 /* Need an extra comparison to get the mask. */
141 mp_limb_t t
= al
- bl
;
142 mp_limb_t bgta
= - (bl
> al
);
147 /* If b > a, invoke reciprocity */
148 bit
^= (bgta
& al
& bl
);
150 /* b <-- min (a, b) */
154 al
= (t
^ bgta
) - bgta
;
156 /* Number of trailing zeros is the same no matter if we look at
157 * t or a, but using t gives more parallelism. */
158 count_trailing_zeros (c
, t
);
160 /* (2/b) = -1 if b = 3 or 5 mod 8 */
161 bit
^= c
& (bl
^ (bl
>> 1));
163 if (UNLIKELY (c
== GMP_NUMB_BITS
))
164 return 1 - 2*(bit
& 1);
169 /* Here we have a little impedance mismatch. Better to inline it? */
170 return mpn_jacobi_base (2*al
+1, 2*bl
+1, bit
<< 1);
172 #elif JACOBI_2_METHOD == 2
174 mpn_jacobi_2 (mp_srcptr ap
, mp_srcptr bp
, unsigned bit
)
176 mp_limb_t ah
, al
, bh
, bl
;
189 if (bh
== 0 && bl
== 1)
191 return 1 - (bit
& 2);
196 /* (0|b) = 0, b > 1 */
199 count_trailing_zeros (c
, ah
);
200 bit
^= ((GMP_NUMB_BITS
+ c
) << 1) & (bl
^ (bl
>> 1));
207 return 1 - (bit
& 2);
217 count_trailing_zeros (c
, al
);
219 al
= ((ah
<< (GMP_NUMB_BITS
- c
)) & GMP_NUMB_MASK
) | (al
>> c
);
221 bit
^= (c
<< 1) & (bl
^ (bl
>> 1));
228 MP_LIMB_T_SWAP (al
, bl
);
240 sub_ddmmss (ah
, al
, ah
, al
, bh
, bl
);
243 count_trailing_zeros (c
, ah
);
244 bit
^= ((GMP_NUMB_BITS
+ c
) << 1) & (bl
^ (bl
>> 1));
253 count_trailing_zeros (c
, al
);
254 bit
^= (c
<< 1) & (bl
^ (bl
>> 1));
255 al
= ((ah
<< (GMP_NUMB_BITS
- c
)) & GMP_NUMB_MASK
) | (al
>> c
);
264 MP_LIMB_T_SWAP (al
, bl
);
274 sub_ddmmss (bh
, bl
, bh
, bl
, ah
, al
);
277 count_trailing_zeros (c
, bh
);
278 bit
^= ((GMP_NUMB_BITS
+ c
) << 1) & (al
^ (al
>> 1));
284 count_trailing_zeros (c
, bl
);
285 bit
^= (c
<< 1) & (al
^ (al
>> 1));
286 bl
= ((bh
<< (GMP_NUMB_BITS
- c
)) & GMP_NUMB_MASK
) | (bl
>> c
);
297 MP_LIMB_T_SWAP (al
, bl
);
304 count_trailing_zeros (c
, al
);
305 bit
^= (c
<< 1) & (bl
^ (bl
>> 1));
309 return 1 - (bit
& 2);
311 MP_LIMB_T_SWAP (al
, bl
);
318 /* Compute (a|b), with b a single limb. */
323 return 1 - (bit
& 2);
333 count_trailing_zeros (c
, ah
);
334 bit
^= ((GMP_NUMB_BITS
+ c
) << 1) & (bl
^ (bl
>> 1));
338 count_trailing_zeros (c
, al
);
340 al
= ((ah
<< (GMP_NUMB_BITS
- c
)) & GMP_NUMB_MASK
) | (al
>> c
);
342 bit
^= (c
<< 1) & (bl
^ (bl
>> 1));
348 return mpn_jacobi_base (al
, bl
, bit
);
351 #error Unsupported value for JACOBI_2_METHOD