1 /* mpn_gcdext -- Extended Greatest Common Divisor.
3 Copyright 1996, 1998, 2000-2005, 2008, 2009 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of either:
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.
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
20 or both in parallel, as here.
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
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/. */
35 #ifndef GCDEXT_1_USE_BINARY
36 #define GCDEXT_1_USE_BINARY 0
39 #ifndef GCDEXT_1_BINARY_METHOD
40 #define GCDEXT_1_BINARY_METHOD 2
47 #if GCDEXT_1_USE_BINARY
50 static unsigned char zerotab
[0x40] = {
51 6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
52 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
53 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
54 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
59 mpn_gcdext_1 (mp_limb_signed_t
*sp
, mp_limb_signed_t
*tp
,
60 mp_limb_t u
, mp_limb_t v
)
67 where U, V are the inputs (without any shared power of two),
68 and the matrix has determinant ± 2^{shift}.
81 #if GCDEXT_1_BINARY_METHOD == 2
88 count_trailing_zeros (zero_bits
, u
| v
);
94 count_trailing_zeros (shift
, u
);
98 else if ((v
& 1) == 0)
100 count_trailing_zeros (shift
, v
);
107 #if GCDEXT_1_BINARY_METHOD == 1
115 count
= zerotab
[u
& 0x3f];
117 if (UNLIKELY (count
== 6))
122 c
= zerotab
[u
& 0x3f];
129 count_trailing_zeros (count
, u
);
132 t0
+= t1
; t1
<<= count
;
133 s0
+= s1
; s1
<<= count
;
139 count
= zerotab
[v
& 0x3f];
141 if (UNLIKELY (count
== 6))
146 c
= zerotab
[v
& 0x3f];
153 count_trailing_zeros (count
, v
);
156 t1
+= t0
; t0
<<= count
;
157 s1
+= s0
; s0
<<= count
;
162 # if GCDEXT_1_BINARY_METHOD == 2
172 mp_limb_t vgtu
= LIMB_HIGHBIT_TO_MASK (d
);
176 /* When v <= u (vgtu == 0), the updates are:
178 (u; v) <-- ( (u - v) >> count; v) (det = +(1<<count) for corr. M factor)
179 (t1, t0) <-- (t1 << count, t0 + t1)
181 and when v > 0, the updates are
183 (u; v) <-- ( (v - u) >> count; u) (det = -(1<<count))
184 (t1, t0) <-- (t0 << count, t0 + t1)
186 and similarly for s1, s0
189 /* v <-- min (u, v) */
193 u
= (d
^ vgtu
) - vgtu
;
195 /* Number of trailing zeros is the same no matter if we look at
196 * d or u, but using d gives more parallelism. */
198 count
= zerotab
[d
& 0x3f];
199 if (UNLIKELY (count
== 6))
205 c
= zerotab
[d
& 0x3f];
211 count_trailing_zeros (count
, d
);
215 tx
= vgtu
& (t0
- t1
);
216 sx
= vgtu
& (s0
- s1
);
229 # else /* GCDEXT_1_BINARY_METHOD == 2 */
230 # error Unknown GCDEXT_1_BINARY_METHOD
234 /* Now u = v = g = gcd (u,v). Compute U/g and V/g */
238 ugh
= ug
/2 + (ug
& 1);
239 vgh
= vg
/2 + (vg
& 1);
241 /* Now ±2^{shift} g = s0 U - t0 V. Get rid of the power of two, using
242 s0 U - t0 V = (s0 + V/g) U - (t0 + U/g) V. */
243 for (i
= 0; i
< shift
; i
++)
245 mp_limb_t mask
= - ( (s0
| t0
) & 1);
252 /* FIXME: Try simplifying this condition. */
253 if ( (s0
> 1 && 2*s0
>= vg
) || (t0
> 1 && 2*t0
>= ug
) )
258 #if GCDEXT_1_BINARY_METHOD == 2
259 /* Conditional negation. */
260 s0
= (s0
^ det_sign
) - det_sign
;
261 t0
= (t0
^ det_sign
) - det_sign
;
266 return u
<< zero_bits
;
269 #else /* !GCDEXT_1_USE_BINARY */
272 /* FIXME: Takes two single-word limbs. It could be extended to a
273 * function that accepts a bignum for the first input, and only
274 * returns the first co-factor. */
277 mpn_gcdext_1 (mp_limb_signed_t
*up
, mp_limb_signed_t
*vp
,
278 mp_limb_t a
, mp_limb_t b
)
285 where A, B are the original inputs.
287 mp_limb_signed_t u0
= 1;
288 mp_limb_signed_t v0
= 0;
289 mp_limb_signed_t u1
= 0;
290 mp_limb_signed_t v1
= 1;
328 #endif /* !GCDEXT_1_USE_BINARY */