1 /* mpn_gcdext -- Extended Greatest Common Divisor.
2
3 Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008, 2009 Free Software
4 Foundation, Inc.
5
6 This file is part of the GNU MP Library.
7
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MP Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
20
21 #include "gmp.h"
22 #include "gmp-impl.h"
23 #include "longlong.h"
24
25 /* Temporary storage: 3*(n+1) for u. n+1 for the matrix-vector
26 multiplications (if hgcd2 succeeds). If hgcd fails, n+1 limbs are
27 needed for the division, with most n for the quotient, and n+1 for
28 the product q u0. In all, 4n + 3. */
29
30 mp_size_t
mpn_gcdext_lehmer_n(mp_ptr gp,mp_ptr up,mp_size_t * usize,mp_ptr ap,mp_ptr bp,mp_size_t n,mp_ptr tp)31 mpn_gcdext_lehmer_n (mp_ptr gp, mp_ptr up, mp_size_t *usize,
32 mp_ptr ap, mp_ptr bp, mp_size_t n,
33 mp_ptr tp)
34 {
35 mp_size_t ualloc = n + 1;
36
37 /* Keeps track of the second row of the reduction matrix
38 *
39 * M = (v0, v1 ; u0, u1)
40 *
41 * which correspond to the first column of the inverse
42 *
43 * M^{-1} = (u1, -v1; -u0, v0)
44 */
45
46 mp_size_t un;
47 mp_ptr u0;
48 mp_ptr u1;
49 mp_ptr u2;
50
51 MPN_ZERO (tp, 3*ualloc);
52 u0 = tp; tp += ualloc;
53 u1 = tp; tp += ualloc;
54 u2 = tp; tp += ualloc;
55
56 u1[0] = 1; un = 1;
57
58 /* FIXME: Handle n == 2 differently, after the loop? */
59 while (n >= 2)
60 {
61 struct hgcd_matrix1 M;
62 mp_limb_t ah, al, bh, bl;
63 mp_limb_t mask;
64
65 mask = ap[n-1] | bp[n-1];
66 ASSERT (mask > 0);
67
68 if (mask & GMP_NUMB_HIGHBIT)
69 {
70 ah = ap[n-1]; al = ap[n-2];
71 bh = bp[n-1]; bl = bp[n-2];
72 }
73 else if (n == 2)
74 {
75 /* We use the full inputs without truncation, so we can
76 safely shift left. */
77 int shift;
78
79 count_leading_zeros (shift, mask);
80 ah = MPN_EXTRACT_NUMB (shift, ap[1], ap[0]);
81 al = ap[0] << shift;
82 bh = MPN_EXTRACT_NUMB (shift, bp[1], bp[0]);
83 bl = bp[0] << shift;
84 }
85 else
86 {
87 int shift;
88
89 count_leading_zeros (shift, mask);
90 ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]);
91 al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]);
92 bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]);
93 bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]);
94 }
95
96 /* Try an mpn_nhgcd2 step */
97 if (mpn_hgcd2 (ah, al, bh, bl, &M))
98 {
99 n = mpn_hgcd_mul_matrix1_inverse_vector (&M, tp, ap, bp, n);
100 MP_PTR_SWAP (ap, tp);
101 un = mpn_hgcd_mul_matrix1_vector(&M, u2, u0, u1, un);
102 MP_PTR_SWAP (u0, u2);
103 }
104 else
105 {
106 /* mpn_hgcd2 has failed. Then either one of a or b is very
107 small, or the difference is very small. Perform one
108 subtraction followed by one division. */
109 mp_size_t gn;
110 mp_size_t updated_un = un;
111
112 /* Temporary storage n for the quotient and ualloc for the
113 new cofactor. */
114 n = mpn_gcdext_subdiv_step (gp, &gn, up, usize, ap, bp, n,
115 u0, u1, &updated_un, tp, u2);
116 if (n == 0)
117 return gn;
118
119 un = updated_un;
120 }
121 }
122 ASSERT_ALWAYS (ap[0] > 0);
123 ASSERT_ALWAYS (bp[0] > 0);
124
125 if (ap[0] == bp[0])
126 {
127 int c;
128
129 /* Which cofactor to return now? Candidates are +u1 and -u0,
130 depending on which of a and b was most recently reduced,
131 which we don't keep track of. So compare and get the smallest
132 one. */
133
134 gp[0] = ap[0];
135
136 MPN_CMP (c, u0, u1, un);
137 ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1));
138 if (c < 0)
139 {
140 MPN_NORMALIZE (u0, un);
141 MPN_COPY (up, u0, un);
142 *usize = -un;
143 }
144 else
145 {
146 MPN_NORMALIZE_NOT_ZERO (u1, un);
147 MPN_COPY (up, u1, un);
148 *usize = un;
149 }
150 return 1;
151 }
152 else
153 {
154 mp_limb_t uh, vh;
155 mp_limb_signed_t u;
156 mp_limb_signed_t v;
157 int negate;
158
159 gp[0] = mpn_gcdext_1 (&u, &v, ap[0], bp[0]);
160
161 /* Set up = u u1 - v u0. Keep track of size, un grows by one or
162 two limbs. */
163
164 if (u == 0)
165 {
166 ASSERT (v == 1);
167 MPN_NORMALIZE (u0, un);
168 MPN_COPY (up, u0, un);
169 *usize = -un;
170 return 1;
171 }
172 else if (v == 0)
173 {
174 ASSERT (u == 1);
175 MPN_NORMALIZE (u1, un);
176 MPN_COPY (up, u1, un);
177 *usize = un;
178 return 1;
179 }
180 else if (u > 0)
181 {
182 negate = 0;
183 ASSERT (v < 0);
184 v = -v;
185 }
186 else
187 {
188 negate = 1;
189 ASSERT (v > 0);
190 u = -u;
191 }
192
193 uh = mpn_mul_1 (up, u1, un, u);
194 vh = mpn_addmul_1 (up, u0, un, v);
195
196 if ( (uh | vh) > 0)
197 {
198 uh += vh;
199 up[un++] = uh;
200 if (uh < vh)
201 up[un++] = 1;
202 }
203
204 MPN_NORMALIZE_NOT_ZERO (up, un);
205
206 *usize = negate ? -un : un;
207 return 1;
208 }
209 }
210