1 /* hgcd2_jacobi.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, 2011 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.h"
36 #include "gmp-impl.h"
37 #include "longlong.h"
38
39 #if GMP_NAIL_BITS > 0
40 #error Nails not supported.
41 #endif
42
43 /* FIXME: Duplicated in hgcd2.c. Should move to gmp-impl.h, and
44 possibly be renamed. */
45 static inline mp_limb_t
div1(mp_ptr rp,mp_limb_t n0,mp_limb_t d0)46 div1 (mp_ptr rp,
47 mp_limb_t n0,
48 mp_limb_t d0)
49 {
50 mp_limb_t q = 0;
51
52 if ((mp_limb_signed_t) n0 < 0)
53 {
54 int cnt;
55 for (cnt = 1; (mp_limb_signed_t) d0 >= 0; cnt++)
56 {
57 d0 = d0 << 1;
58 }
59
60 q = 0;
61 while (cnt)
62 {
63 q <<= 1;
64 if (n0 >= d0)
65 {
66 n0 = n0 - d0;
67 q |= 1;
68 }
69 d0 = d0 >> 1;
70 cnt--;
71 }
72 }
73 else
74 {
75 int cnt;
76 for (cnt = 0; n0 >= d0; cnt++)
77 {
78 d0 = d0 << 1;
79 }
80
81 q = 0;
82 while (cnt)
83 {
84 d0 = d0 >> 1;
85 q <<= 1;
86 if (n0 >= d0)
87 {
88 n0 = n0 - d0;
89 q |= 1;
90 }
91 cnt--;
92 }
93 }
94 *rp = n0;
95 return q;
96 }
97
98 /* Two-limb division optimized for small quotients. */
99 static inline mp_limb_t
div2(mp_ptr rp,mp_limb_t nh,mp_limb_t nl,mp_limb_t dh,mp_limb_t dl)100 div2 (mp_ptr rp,
101 mp_limb_t nh, mp_limb_t nl,
102 mp_limb_t dh, mp_limb_t dl)
103 {
104 mp_limb_t q = 0;
105
106 if ((mp_limb_signed_t) nh < 0)
107 {
108 int cnt;
109 for (cnt = 1; (mp_limb_signed_t) dh >= 0; cnt++)
110 {
111 dh = (dh << 1) | (dl >> (GMP_LIMB_BITS - 1));
112 dl = dl << 1;
113 }
114
115 while (cnt)
116 {
117 q <<= 1;
118 if (nh > dh || (nh == dh && nl >= dl))
119 {
120 sub_ddmmss (nh, nl, nh, nl, dh, dl);
121 q |= 1;
122 }
123 dl = (dh << (GMP_LIMB_BITS - 1)) | (dl >> 1);
124 dh = dh >> 1;
125 cnt--;
126 }
127 }
128 else
129 {
130 int cnt;
131 for (cnt = 0; nh > dh || (nh == dh && nl >= dl); cnt++)
132 {
133 dh = (dh << 1) | (dl >> (GMP_LIMB_BITS - 1));
134 dl = dl << 1;
135 }
136
137 while (cnt)
138 {
139 dl = (dh << (GMP_LIMB_BITS - 1)) | (dl >> 1);
140 dh = dh >> 1;
141 q <<= 1;
142 if (nh > dh || (nh == dh && nl >= dl))
143 {
144 sub_ddmmss (nh, nl, nh, nl, dh, dl);
145 q |= 1;
146 }
147 cnt--;
148 }
149 }
150
151 rp[0] = nl;
152 rp[1] = nh;
153
154 return q;
155 }
156
157 int
mpn_hgcd2_jacobi(mp_limb_t ah,mp_limb_t al,mp_limb_t bh,mp_limb_t bl,struct hgcd_matrix1 * M,unsigned * bitsp)158 mpn_hgcd2_jacobi (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl,
159 struct hgcd_matrix1 *M, unsigned *bitsp)
160 {
161 mp_limb_t u00, u01, u10, u11;
162 unsigned bits = *bitsp;
163
164 if (ah < 2 || bh < 2)
165 return 0;
166
167 if (ah > bh || (ah == bh && al > bl))
168 {
169 sub_ddmmss (ah, al, ah, al, bh, bl);
170 if (ah < 2)
171 return 0;
172
173 u00 = u01 = u11 = 1;
174 u10 = 0;
175 bits = mpn_jacobi_update (bits, 1, 1);
176 }
177 else
178 {
179 sub_ddmmss (bh, bl, bh, bl, ah, al);
180 if (bh < 2)
181 return 0;
182
183 u00 = u10 = u11 = 1;
184 u01 = 0;
185 bits = mpn_jacobi_update (bits, 0, 1);
186 }
187
188 if (ah < bh)
189 goto subtract_a;
190
191 for (;;)
192 {
193 ASSERT (ah >= bh);
194 if (ah == bh)
195 goto done;
196
197 if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
198 {
199 ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
200 bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
201
202 break;
203 }
204
205 /* Subtract a -= q b, and multiply M from the right by (1 q ; 0
206 1), affecting the second column of M. */
207 ASSERT (ah > bh);
208 sub_ddmmss (ah, al, ah, al, bh, bl);
209
210 if (ah < 2)
211 goto done;
212
213 if (ah <= bh)
214 {
215 /* Use q = 1 */
216 u01 += u00;
217 u11 += u10;
218 bits = mpn_jacobi_update (bits, 1, 1);
219 }
220 else
221 {
222 mp_limb_t r[2];
223 mp_limb_t q = div2 (r, ah, al, bh, bl);
224 al = r[0]; ah = r[1];
225 if (ah < 2)
226 {
227 /* A is too small, but q is correct. */
228 u01 += q * u00;
229 u11 += q * u10;
230 bits = mpn_jacobi_update (bits, 1, q & 3);
231 goto done;
232 }
233 q++;
234 u01 += q * u00;
235 u11 += q * u10;
236 bits = mpn_jacobi_update (bits, 1, q & 3);
237 }
238 subtract_a:
239 ASSERT (bh >= ah);
240 if (ah == bh)
241 goto done;
242
243 if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
244 {
245 ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
246 bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
247
248 goto subtract_a1;
249 }
250
251 /* Subtract b -= q a, and multiply M from the right by (1 0 ; q
252 1), affecting the first column of M. */
253 sub_ddmmss (bh, bl, bh, bl, ah, al);
254
255 if (bh < 2)
256 goto done;
257
258 if (bh <= ah)
259 {
260 /* Use q = 1 */
261 u00 += u01;
262 u10 += u11;
263 bits = mpn_jacobi_update (bits, 0, 1);
264 }
265 else
266 {
267 mp_limb_t r[2];
268 mp_limb_t q = div2 (r, bh, bl, ah, al);
269 bl = r[0]; bh = r[1];
270 if (bh < 2)
271 {
272 /* B is too small, but q is correct. */
273 u00 += q * u01;
274 u10 += q * u11;
275 bits = mpn_jacobi_update (bits, 0, q & 3);
276 goto done;
277 }
278 q++;
279 u00 += q * u01;
280 u10 += q * u11;
281 bits = mpn_jacobi_update (bits, 0, q & 3);
282 }
283 }
284
285 /* NOTE: Since we discard the least significant half limb, we don't
286 get a truly maximal M (corresponding to |a - b| <
287 2^{GMP_LIMB_BITS +1}). */
288 /* Single precision loop */
289 for (;;)
290 {
291 ASSERT (ah >= bh);
292 if (ah == bh)
293 break;
294
295 ah -= bh;
296 if (ah < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
297 break;
298
299 if (ah <= bh)
300 {
301 /* Use q = 1 */
302 u01 += u00;
303 u11 += u10;
304 bits = mpn_jacobi_update (bits, 1, 1);
305 }
306 else
307 {
308 mp_limb_t r;
309 mp_limb_t q = div1 (&r, ah, bh);
310 ah = r;
311 if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
312 {
313 /* A is too small, but q is correct. */
314 u01 += q * u00;
315 u11 += q * u10;
316 bits = mpn_jacobi_update (bits, 1, q & 3);
317 break;
318 }
319 q++;
320 u01 += q * u00;
321 u11 += q * u10;
322 bits = mpn_jacobi_update (bits, 1, q & 3);
323 }
324 subtract_a1:
325 ASSERT (bh >= ah);
326 if (ah == bh)
327 break;
328
329 bh -= ah;
330 if (bh < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
331 break;
332
333 if (bh <= ah)
334 {
335 /* Use q = 1 */
336 u00 += u01;
337 u10 += u11;
338 bits = mpn_jacobi_update (bits, 0, 1);
339 }
340 else
341 {
342 mp_limb_t r;
343 mp_limb_t q = div1 (&r, bh, ah);
344 bh = r;
345 if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
346 {
347 /* B is too small, but q is correct. */
348 u00 += q * u01;
349 u10 += q * u11;
350 bits = mpn_jacobi_update (bits, 0, q & 3);
351 break;
352 }
353 q++;
354 u00 += q * u01;
355 u10 += q * u11;
356 bits = mpn_jacobi_update (bits, 0, q & 3);
357 }
358 }
359
360 done:
361 M->u[0][0] = u00; M->u[0][1] = u01;
362 M->u[1][0] = u10; M->u[1][1] = u11;
363 *bitsp = bits;
364
365 return 1;
366 }
367