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 2003, 2004, 2005, 2008 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 the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
16 The GNU MP Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 License for more details.
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
28 /* Use binary algorithm to compute G <-- GCD (U, V) for usize, vsize == 2.
29 Both U and V must be odd. */
30 static inline mp_size_t
31 gcd_2 (mp_ptr gp, mp_srcptr up, mp_srcptr vp)
33 mp_limb_t u0, u1, v0, v1;
44 /* Check for u0 != v0 needed to ensure that argument to
45 * count_trailing_zeros is non-zero. */
46 while (u1 != v1 && u0 != v0)
52 u0 = (u0 - v0) & GMP_NUMB_MASK;
53 count_trailing_zeros (r, u0);
54 u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
60 v0 = (v0 - u0) & GMP_NUMB_MASK;
61 count_trailing_zeros (r, v0);
62 v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
67 gp[0] = u0, gp[1] = u1, gn = 1 + (u1 != 0);
69 /* If U == V == GCD, done. Otherwise, compute GCD (V, |U - V|). */
70 if (u1 == v1 && u0 == v0)
73 v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
74 gp[0] = mpn_gcd_1 (gp, gn, v0);
79 /* Temporary storage: n */
81 mpn_gcd_lehmer_n (mp_ptr gp, mp_ptr ap, mp_ptr bp, mp_size_t n, mp_ptr tp)
83 /* Relax this requirement, and normalize at the start? Must disallow
85 ASSERT(ap[n-1] > 0 || bp[n-1] > 0);
89 struct hgcd_matrix1 M;
90 mp_limb_t ah, al, bh, bl;
93 mask = ap[n-1] | bp[n-1];
96 if (mask & GMP_NUMB_HIGHBIT)
98 ah = ap[n-1]; al = ap[n-2];
99 bh = bp[n-1]; bl = bp[n-2];
105 count_leading_zeros (shift, mask);
106 ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]);
107 al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]);
108 bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]);
109 bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]);
112 /* Try an mpn_nhgcd2 step */
113 if (mpn_hgcd2 (ah, al, bh, bl, &M))
115 n = mpn_hgcd_mul_matrix1_inverse_vector (&M, tp, ap, bp, n);
116 MP_PTR_SWAP (ap, tp);
120 /* mpn_hgcd2 has failed. Then either one of a or b is very
121 small, or the difference is very small. Perform one
122 subtraction followed by one division. */
125 /* Temporary storage n */
126 n = mpn_gcd_subdiv_step (gp, &gn, ap, bp, n, tp);
134 *gp = mpn_gcd_1(ap, 1, bp[0]);
138 /* Due to the calling convention for mpn_gcd, at most one can be
142 MP_PTR_SWAP (ap, bp);
148 *gp = mpn_gcd_1 (ap, 2, bp[1]);
151 else if (! (bp[0] & 1))
154 count_trailing_zeros (r, bp[0]);
155 bp[0] = ((bp[1] << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (bp[0] >> r);
159 return gcd_2(gp, ap, bp);