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