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, 2020 Free Software
8 Foundation, Inc.
9
10 This file is part of the GNU MP Library.
11
12 The GNU MP Library is free software; you can redistribute it and/or modify
13 it under the terms of either:
14
15 * the GNU Lesser General Public License as published by the Free
16 Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
18
19 or
20
21 * the GNU General Public License as published by the Free Software
22 Foundation; either version 2 of the License, or (at your option) any
23 later version.
24
25 or both in parallel, as here.
26
27 The GNU MP Library is distributed in the hope that it will be useful, but
28 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
29 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
30 for more details.
31
32 You should have received copies of the GNU General Public License and the
33 GNU Lesser General Public License along with the GNU MP Library. If not,
34 see https://www.gnu.org/licenses/. */
35
36 #include "gmp-impl.h"
37 #include "longlong.h"
38
39 #include "mpn/generic/hgcd2-div.h"
40
41 #if GMP_NAIL_BITS != 0
42 #error Nails not implemented
43 #endif
44
45 int
46 mpn_hgcd2_jacobi (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl,
47 struct hgcd_matrix1 *M, unsigned *bitsp)
48 {
49 mp_limb_t u00, u01, u10, u11;
50 unsigned bits = *bitsp;
51
52 if (ah < 2 || bh < 2)
53 return 0;
54
55 if (ah > bh || (ah == bh && al > bl))
56 {
57 sub_ddmmss (ah, al, ah, al, bh, bl);
58 if (ah < 2)
59 return 0;
60
61 u00 = u01 = u11 = 1;
62 u10 = 0;
63 bits = mpn_jacobi_update (bits, 1, 1);
64 }
65 else
66 {
67 sub_ddmmss (bh, bl, bh, bl, ah, al);
68 if (bh < 2)
69 return 0;
70
71 u00 = u10 = u11 = 1;
72 u01 = 0;
73 bits = mpn_jacobi_update (bits, 0, 1);
74 }
75
76 if (ah < bh)
77 goto subtract_a;
78
79 for (;;)
80 {
81 ASSERT (ah >= bh);
82 if (ah == bh)
83 goto done;
84
85 if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
86 {
87 ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
88 bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
89
90 break;
91 }
92
93 /* Subtract a -= q b, and multiply M from the right by (1 q ; 0
94 1), affecting the second column of M. */
95 ASSERT (ah > bh);
96 sub_ddmmss (ah, al, ah, al, bh, bl);
97
98 if (ah < 2)
99 goto done;
100
101 if (ah <= bh)
102 {
103 /* Use q = 1 */
104 u01 += u00;
105 u11 += u10;
106 bits = mpn_jacobi_update (bits, 1, 1);
107 }
108 else
109 {
110 mp_limb_t r[2];
111 mp_limb_t q = div2 (r, ah, al, bh, bl);
112 al = r[0]; ah = r[1];
113 if (ah < 2)
114 {
115 /* A is too small, but q is correct. */
116 u01 += q * u00;
117 u11 += q * u10;
118 bits = mpn_jacobi_update (bits, 1, q & 3);
119 goto done;
120 }
121 q++;
122 u01 += q * u00;
123 u11 += q * u10;
124 bits = mpn_jacobi_update (bits, 1, q & 3);
125 }
126 subtract_a:
127 ASSERT (bh >= ah);
128 if (ah == bh)
129 goto done;
130
131 if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
132 {
133 ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
134 bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
135
136 goto subtract_a1;
137 }
138
139 /* Subtract b -= q a, and multiply M from the right by (1 0 ; q
140 1), affecting the first column of M. */
141 sub_ddmmss (bh, bl, bh, bl, ah, al);
142
143 if (bh < 2)
144 goto done;
145
146 if (bh <= ah)
147 {
148 /* Use q = 1 */
149 u00 += u01;
150 u10 += u11;
151 bits = mpn_jacobi_update (bits, 0, 1);
152 }
153 else
154 {
155 mp_limb_t r[2];
156 mp_limb_t q = div2 (r, bh, bl, ah, al);
157 bl = r[0]; bh = r[1];
158 if (bh < 2)
159 {
160 /* B is too small, but q is correct. */
161 u00 += q * u01;
162 u10 += q * u11;
163 bits = mpn_jacobi_update (bits, 0, q & 3);
164 goto done;
165 }
166 q++;
167 u00 += q * u01;
168 u10 += q * u11;
169 bits = mpn_jacobi_update (bits, 0, q & 3);
170 }
171 }
172
173 /* NOTE: Since we discard the least significant half limb, we don't get a
174 truly maximal M (corresponding to |a - b| < 2^{GMP_LIMB_BITS +1}). */
175 /* Single precision loop */
176 for (;;)
177 {
178 ASSERT (ah >= bh);
179
180 ah -= bh;
181 if (ah < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
182 break;
183
184 if (ah <= bh)
185 {
186 /* Use q = 1 */
187 u01 += u00;
188 u11 += u10;
189 bits = mpn_jacobi_update (bits, 1, 1);
190 }
191 else
192 {
193 mp_double_limb_t rq = div1 (ah, bh);
194 mp_limb_t q = rq.d1;
195 ah = rq.d0;
196
197 if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
198 {
199 /* A is too small, but q is correct. */
200 u01 += q * u00;
201 u11 += q * u10;
202 bits = mpn_jacobi_update (bits, 1, q & 3);
203 break;
204 }
205 q++;
206 u01 += q * u00;
207 u11 += q * u10;
208 bits = mpn_jacobi_update (bits, 1, q & 3);
209 }
210 subtract_a1:
211 ASSERT (bh >= ah);
212
213 bh -= ah;
214 if (bh < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
215 break;
216
217 if (bh <= ah)
218 {
219 /* Use q = 1 */
220 u00 += u01;
221 u10 += u11;
222 bits = mpn_jacobi_update (bits, 0, 1);
223 }
224 else
225 {
226 mp_double_limb_t rq = div1 (bh, ah);
227 mp_limb_t q = rq.d1;
228 bh = rq.d0;
229
230 if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
231 {
232 /* B is too small, but q is correct. */
233 u00 += q * u01;
234 u10 += q * u11;
235 bits = mpn_jacobi_update (bits, 0, q & 3);
236 break;
237 }
238 q++;
239 u00 += q * u01;
240 u10 += q * u11;
241 bits = mpn_jacobi_update (bits, 0, q & 3);
242 }
243 }
244
245 done:
246 M->u[0][0] = u00; M->u[0][1] = u01;
247 M->u[1][0] = u10; M->u[1][1] = u11;
248 *bitsp = bits;
249
250 return 1;
251 }