Line data Source code
1 : #line 2 "../src/kernel/none/halfgcd.c"
2 : /* Copyright (C) 2019 The PARI group.
3 :
4 : This file is part of the PARI/GP package.
5 :
6 : PARI/GP is free software; you can redistribute it and/or modify it under the
7 : terms of the GNU General Public License as published by the Free Software
8 : Foundation; either version 2 of the License, or (at your option) any later
9 : version. It is distributed in the hope that it will be useful, but WITHOUT
10 : ANY WARRANTY WHATSOEVER.
11 :
12 : Check the License for details. You should have received a copy of it, along
13 : with the package; see the file 'COPYING'. If not, write to the Free Software
14 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
15 :
16 : GEN
17 7132248 : ZM2_mul(GEN A, GEN B)
18 : {
19 7132248 : const long t = ZM2_MUL_LIMIT+2;
20 7132248 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
21 7132248 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
22 7132248 : if (lgefint(A11) < t || lgefint(B11) < t || lgefint(A22) < t || lgefint(B22) < t
23 110656 : || lgefint(A12) < t || lgefint(B12) < t || lgefint(A21) < t || lgefint(B21) < t)
24 : {
25 7024250 : GEN a = mulii(A11, B11), b = mulii(A12, B21);
26 7023731 : GEN c = mulii(A11, B12), d = mulii(A12, B22);
27 7023681 : GEN e = mulii(A21, B11), f = mulii(A22, B21);
28 7023790 : GEN g = mulii(A21, B12), h = mulii(A22, B22);
29 7023710 : retmkmat2(mkcol2(addii(a,b), addii(e,f)), mkcol2(addii(c,d), addii(g,h)));
30 : } else
31 : {
32 107998 : GEN M1 = mulii(addii(A11,A22), addii(B11,B22));
33 107998 : GEN M2 = mulii(addii(A21,A22), B11);
34 107998 : GEN M3 = mulii(A11, subii(B12,B22));
35 107998 : GEN M4 = mulii(A22, subii(B21,B11));
36 107998 : GEN M5 = mulii(addii(A11,A12), B22);
37 107998 : GEN M6 = mulii(subii(A21,A11), addii(B11,B12));
38 107998 : GEN M7 = mulii(subii(A12,A22), addii(B21,B22));
39 107998 : GEN T1 = addii(M1,M4), T2 = subii(M7,M5);
40 107998 : GEN T3 = subii(M1,M2), T4 = addii(M3,M6);
41 107997 : retmkmat2(mkcol2(addii(T1,T2), addii(M2,M4)),
42 : mkcol2(addii(M3,M5), addii(T3,T4)));
43 : }
44 : }
45 :
46 : static GEN
47 666 : matid2(void)
48 : {
49 666 : retmkmat2(mkcol2(gen_1,gen_0),
50 : mkcol2(gen_0,gen_1));
51 : }
52 :
53 : /* Return M*[q,1;1,0] */
54 : static GEN
55 2983043 : mulq(GEN M, GEN q)
56 : {
57 2983043 : GEN u, v, res = cgetg(3, t_MAT);
58 2983014 : u = addii(mulii(gcoeff(M,1,1), q), gcoeff(M,1,2));
59 2983055 : v = addii(mulii(gcoeff(M,2,1), q), gcoeff(M,2,2));
60 2983059 : gel(res,1) = mkcol2(u, v);
61 2983028 : gel(res,2) = gel(M,1);
62 2983028 : return res;
63 : }
64 : static GEN
65 59 : mulqab(GEN M, GEN q, GEN *ap, GEN *bp)
66 : {
67 59 : GEN b = subii(*ap, mulii(*bp, q));
68 59 : *ap = *bp; *bp = b;
69 59 : return mulq(M,q);
70 : }
71 :
72 : /* Return M*[q,1;1,0]^-1 */
73 :
74 : static GEN
75 1169 : mulqi(GEN M, GEN q, GEN *ap, GEN *bp)
76 : {
77 : GEN u, v, res, a;
78 1169 : a = addii(mulii(*ap, q), *bp);
79 1169 : *bp = *ap; *ap = a;
80 1169 : res = cgetg(3, t_MAT);
81 1169 : u = subii(gcoeff(M,1,1),mulii(gcoeff(M,1,2), q));
82 1169 : v = subii(gcoeff(M,2,1),mulii(gcoeff(M,2,2), q));
83 1169 : gel(res,1) = gel(M,2);
84 1169 : gel(res,2) = mkcol2(u,v);
85 1169 : return res;
86 : }
87 :
88 : /* test whether n is a power of 2 */
89 : static long
90 23200698 : isint2n(GEN n)
91 : {
92 : GEN x;
93 23200698 : long lx = lgefint(n), i;
94 23200698 : if (lx == 2) return 0;
95 23200698 : x = int_MSW(n);
96 23200698 : if (*(ulong*)x != 1UL<<expu(*(ulong*)x) ) return 0;
97 519940 : for (i = 3; i < lx; i++)
98 : {
99 512562 : x = int_precW(x); if (*x) return 0;
100 : }
101 7378 : return 1;
102 : }
103 :
104 : static long
105 23200868 : uexpi(GEN a)
106 23200868 : { return expi(a)+!isint2n(a); }
107 :
108 : static GEN
109 107718 : FIXUP0(GEN M, GEN *a, GEN *b, long m)
110 : {
111 107718 : long cnt=0;
112 137130 : while (expi(*b) >= m)
113 : {
114 29412 : GEN r, q = dvmdii(*a, *b, &r);
115 29412 : *a = *b; *b = r;
116 29412 : M = mulq(M, q);
117 29412 : cnt++;
118 : };
119 107718 : if (cnt>6) pari_err_BUG("FIXUP0");
120 107718 : return M;
121 : }
122 :
123 : static long
124 6397137 : signdet(GEN Q)
125 : {
126 6397137 : long a = Mod4(gcoeff(Q,1,1)), b = Mod4(gcoeff(Q,1,2));
127 6397089 : long c = Mod4(gcoeff(Q,2,1)), d = Mod4(gcoeff(Q,2,2));
128 6398325 : return ((a*d-b*c)&3)==1 ? 1 : -1;
129 : }
130 :
131 : static GEN
132 6231786 : ZM_inv2(GEN M)
133 : {
134 6231786 : long e = signdet(M);
135 9867364 : if (e==1) return mkmat22(gcoeff(M,2,2),negi(gcoeff(M,1,2)),
136 3634766 : negi(gcoeff(M,2,1)),gcoeff(M,1,1));
137 2597949 : else return mkmat22(negi(gcoeff(M,2,2)),gcoeff(M,1,2),
138 2597942 : gcoeff(M,2,1),negi(gcoeff(M,1,1)));
139 : }
140 :
141 : static GEN
142 1169 : lastq(GEN Q)
143 : {
144 1169 : GEN p = gcoeff(Q,1,1), q = gcoeff(Q,1,2), s = gcoeff(Q,2,2);
145 1169 : if (signe(q)==0) pari_err_BUG("halfgcd");
146 1169 : if (signe(s)==0) return p;
147 1169 : if (equali1(q)) return subiu(p,1);
148 1169 : return divii(p, q);
149 : }
150 :
151 : static GEN
152 6820 : mulT(GEN Q, GEN *ap, GEN *bp)
153 : {
154 6820 : *ap = addii(*ap, *bp);
155 6820 : *bp = negi(*bp);
156 6820 : return mkmat2(gel(Q,1),
157 6820 : mkcol2(subii(gcoeff(Q,1,1), gcoeff(Q,1,2))
158 6820 : , subii(gcoeff(Q,2,1), gcoeff(Q,2,2))));
159 : }
160 :
161 : static GEN
162 165516 : FIXUP1(GEN M, GEN a, GEN b, long m, long t, GEN *ap, GEN *bp)
163 : {
164 165516 : GEN Q = gel(M,1), a0 = gel(M,2), b0 = gel(M,3);
165 165516 : GEN q, am = remi2n(a, m), bm = remi2n(b, m);
166 165586 : if (signdet(Q)==-1)
167 : {
168 110956 : *ap = subii(mulii(bm, gcoeff(Q,1,2)),mulii(am, gcoeff(Q,2,2)));
169 110970 : *bp = subii(mulii(am, gcoeff(Q,2,1)),mulii(bm, gcoeff(Q,1,1)));
170 111143 : *ap = addii(*ap, shifti(addii(a0, gcoeff(Q,2,2)), m));
171 111088 : *bp = addii(*bp, shifti(subii(b0, gcoeff(Q,2,1)), m));
172 110760 : if (signe(*bp) >= 0)
173 103869 : return Q;
174 6891 : if (expi(addii(*ap,*bp)) >= m+t)
175 6820 : return mulT(Q, ap ,bp);
176 71 : q = lastq(Q);
177 71 : Q = mulqi(Q, q, ap, bp);
178 71 : if (cmpiu(q, 2)>=0)
179 59 : return mulqab(Q, subiu(q,1), ap, bp);
180 : else
181 12 : return mulqi(Q, lastq(Q), ap, bp);
182 : }
183 : else
184 : {
185 54499 : *ap = subii(mulii(am, gcoeff(Q,2,2)),mulii(bm, gcoeff(Q,1,2)));
186 54499 : *bp = subii(mulii(bm, gcoeff(Q,1,1)),mulii(am, gcoeff(Q,2,1)));
187 54499 : *ap = addii(*ap, shifti(subii(a0, gcoeff(Q,2,2)), m));
188 54499 : *bp = addii(*bp, shifti(addii(b0, gcoeff(Q,2,1)), m));
189 54499 : if (expi(*ap) >= m+t)
190 53413 : return FIXUP0(Q, ap, bp, m+t);
191 : else
192 1086 : return signe(gcoeff(Q,1,2))==0? Q: mulqi(Q, lastq(Q), ap, bp);
193 : }
194 : }
195 :
196 : static long
197 23146623 : magic_threshold(GEN a)
198 23146623 : { return (3+uexpi(a))>>1; }
199 :
200 : static GEN
201 15014693 : HGCD_basecase(GEN y, GEN x)
202 : {
203 15014693 : pari_sp av = avma;
204 : GEN d, d1, q, r;
205 : GEN u, u1, v, v1;
206 : ulong xu, xu1, xv, xv1; /* Lehmer stage recurrence matrix */
207 : int lhmres; /* Lehmer stage return value */
208 :
209 15014693 : long m = magic_threshold(y);
210 :
211 : /* There is no special case for single-word numbers since this is
212 : * mainly meant to be used with large moduli. */
213 15014678 : if (cmpii(y,x) <= 0)
214 : {
215 91605 : d = x; d1 = y;
216 91605 : u = gen_1; u1 = gen_0;
217 91605 : v = gen_0; v1 = gen_1;
218 : } else
219 : {
220 14923047 : d = y; d1 = x;
221 14923047 : u = gen_0; u1 = gen_1;
222 14923047 : v = gen_1; v1 = gen_0;
223 : }
224 33336225 : while (lgefint(d) > 3 && expi(d1) >= m + BITS_IN_LONG + 1)
225 : {
226 : /* do a Lehmer-Jebelean round */
227 18654082 : lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, 0);
228 :
229 18659946 : if (lhmres)
230 : {
231 17692620 : if (lhmres == 1 || lhmres == -1)
232 : {
233 447555 : if (xv1 == 1)
234 : {
235 384983 : r = subii(d,d1); d = d1; d1 = r;
236 384960 : r = addii(u,u1); u = u1; u1 = r;
237 384891 : r = addii(v,v1); v = v1; v1 = r;
238 : }
239 : else
240 : {
241 62572 : r = subii(d, mului(xv1,d1)); d = d1; d1 = r;
242 62571 : r = addii(u, mului(xv1,u1)); u = u1; u1 = r;
243 62573 : r = addii(v, mului(xv1,v1)); v = v1; v1 = r;
244 : }
245 : }
246 : else
247 : {
248 17245065 : r = subii(muliu(d,xu), muliu(d1,xv));
249 17241049 : d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
250 17241059 : r = addii(muliu(u,xu), muliu(u1,xv));
251 17239957 : u1 = addii(muliu(u,xu1), muliu(u1,xv1)); u = r;
252 17239910 : r = addii(muliu(v,xu), muliu(v1,xv));
253 17239781 : v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
254 17242287 : if (lhmres&1) togglesign(d); else togglesign(d1);
255 : }
256 : } /* lhmres != 0 */
257 18656930 : if (expi(d1) < m) break;
258 :
259 18321812 : if (lhmres <= 0 && signe(d1))
260 : {
261 1013105 : q = dvmdii(d,d1,&r);
262 1013048 : d = d1; d1 = r;
263 1013048 : r = addii(u, mulii(q,u1)); u = u1; u1 = r;
264 1012775 : r = addii(v, mulii(q,v1)); v = v1; v1 = r;
265 : }
266 18321467 : if (gc_needed(av,1))
267 : {
268 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ratlift");
269 0 : gerepileall(av, 6, &d, &d1, &u, &u1, &v, &v1);
270 : }
271 : }
272 82375112 : while (expi(d1) >= m)
273 : {
274 67360094 : GEN r, q = dvmdii(d,d1, &r);
275 67365682 : d = d1; d1 = r; swap(u,u1); swap(v,v1);
276 67365682 : u1 = addii(mulii(u, q), u1);
277 67357947 : v1 = addii(mulii(v, q), v1);
278 : }
279 15007564 : return gerepilecopy(av, mkvec3(mkmat22(u1,u,v1,v), d, d1));
280 : }
281 :
282 : static GEN HGCD(GEN x, GEN y);
283 :
284 : /*
285 : Based on
286 : Klaus Thull and Chee K. Yap,
287 : A unified approach to HGCD algorithms for polynomials andintegers,
288 : 1990, Manuscript.
289 : URL: http://cs.nyu.edu/cs/faculty/yap/papers.
290 : */
291 :
292 : static GEN
293 111581 : HGCD_split(GEN a, GEN b)
294 : {
295 111581 : pari_sp av = avma;
296 111581 : long m = magic_threshold(a), t, l, k, tp;
297 : GEN a0, b0, ap, bp, c, d, c0, d0, cp, dp, R, S, T, q, r;
298 111519 : if (signe(b) < 0 || cmpii(a,b)<0) pari_err_BUG("HGCD_split");
299 111563 : if (expi(b) < m)
300 444 : return gerepilecopy(av, mkvec3(matid2(), a, b));
301 111077 : a0 = addiu(shifti(a, -m), 1);
302 111114 : if (cmpiu(a0,7) <= 0)
303 : {
304 0 : R = FIXUP0(matid2(), &a, &b, m);
305 0 : return gerepilecopy(av, mkvec3(R, a, b));
306 : }
307 111090 : b0 = shifti(b,-m);
308 111311 : t = magic_threshold(a0);
309 111107 : R = FIXUP1(HGCD(a0,b0),a, b, m, t, &ap, &bp);
310 110953 : if (expi(bp) < m)
311 56571 : return gerepilecopy(av, mkvec3(R, ap, bp));
312 54305 : q = dvmdii(ap, bp, &r);
313 54305 : c = bp; d = r;
314 54305 : if (cmpiu(shifti(c,-m),6) <= 0)
315 : {
316 21 : R = FIXUP0(mulq(R, q), &c, &d, m);
317 21 : return gerepilecopy(av, mkvec3(R, c, d));
318 : }
319 54284 : l = uexpi(c);
320 54284 : k = 2*m-l-1; if (k<0) pari_err_BUG("halfgcd");
321 54284 : c0 = addiu(shifti(c, -k), 1); if (cmpiu(c0,8)<0) pari_err_BUG("halfgcd");
322 54284 : d0 = shifti(d, -k);
323 54284 : tp = magic_threshold(c0);
324 54284 : S = FIXUP1(HGCD(c0,d0), c, d, k, tp, &cp, &dp);
325 54284 : if (!(expi(cp)>=m+1 && m+1 > expi(dp))) pari_err_BUG("halfgcd");
326 54284 : T = FIXUP0(ZM2_mul(mulq(R, q), S), &cp, &dp, m);
327 54284 : return gerepilecopy(av, mkvec3(T, cp, dp));
328 : }
329 :
330 : static GEN
331 15126112 : HGCD(GEN x, GEN y)
332 : {
333 15126112 : if (lgefint(y)-2 < HALFGCD_LIMIT)
334 15014649 : return HGCD_basecase(x, y);
335 : else
336 111463 : return HGCD_split(x, y);
337 : }
338 :
339 : static GEN
340 29551704 : HGCD0(GEN x, GEN y)
341 : {
342 29551704 : if (signe(y) >= 0 && cmpii(x, y) >= 0)
343 14956983 : return HGCD(x, y);
344 14594712 : if (cmpii(x, y) < 0)
345 : {
346 12059327 : GEN M = HGCD0(y, x), Q = gel(M,1);
347 12057976 : return mkvec3(mkmat22(gcoeff(Q,2,1),gcoeff(Q,2,2),gcoeff(Q,1,1),gcoeff(Q,1,2)),
348 12059597 : gel(M,2),gel(M,3));
349 : } /* Now y <= x*/
350 2535560 : if (signe(x) <= 0)
351 : { /* y <= x <=0 */
352 3847 : GEN M = HGCD(negi(y), negi(x)), Q = gel(M,1);
353 3847 : return mkvec3(mkmat22(negi(gcoeff(Q,2,1)),negi(gcoeff(Q,2,2)),
354 3847 : negi(gcoeff(Q,1,1)),negi(gcoeff(Q,1,2))),
355 3847 : gel(M,2),gel(M,3));
356 : }
357 : else /* y <= 0 <=x */
358 : {
359 2531713 : GEN M = HGCD0(x, negi(y)), Q = gel(M,1);
360 2531713 : return mkvec3(mkmat22(gcoeff(Q,1,1),gcoeff(Q,1,2),negi(gcoeff(Q,2,1)),negi(gcoeff(Q,2,2))),
361 2531713 : gel(M,2),gel(M,3));
362 : }
363 : }
364 :
365 : GEN
366 6232403 : halfgcdii(GEN A, GEN B)
367 : {
368 6232403 : pari_sp av = avma;
369 6232403 : GEN M, Q, a, b, m = abscmpii(A, B)>0 ? A: B;
370 6232401 : M = HGCD0(A,B); Q = gel(M,1); a = gel(M,2); b = gel(M,3);
371 9130045 : while (signe(b) && abscmpii(sqri(b), m) >= 0)
372 : {
373 2899148 : GEN r, q = dvmdii(a, b, &r);
374 2899116 : a = b; b = r;
375 2899116 : Q = mulq(Q, q);
376 : }
377 6231495 : return gerepilecopy(av, mkvec2(ZM_inv2(Q),mkcol2(a,b)));
378 : }
|
