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 6454178 : ZM2_mul(GEN A, GEN B)
18 : {
19 6454178 : const long t = 52;
20 6454178 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
21 6454178 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
22 6454178 : if (lgefint(A11) < t || lgefint(B11) < t || lgefint(A22) < t || lgefint(B22) < t
23 33364 : || lgefint(A12) < t || lgefint(B12) < t || lgefint(A21) < t || lgefint(B21) < t)
24 : {
25 6420815 : GEN a = mulii(A11, B11), b = mulii(A12, B21);
26 6419970 : GEN c = mulii(A11, B12), d = mulii(A12, B22);
27 6419990 : GEN e = mulii(A21, B11), f = mulii(A22, B21);
28 6420094 : GEN g = mulii(A21, B12), h = mulii(A22, B22);
29 6419894 : retmkmat2(mkcol2(addii(a,b), addii(e,f)), mkcol2(addii(c,d), addii(g,h)));
30 : } else
31 : {
32 33363 : GEN M1 = mulii(addii(A11,A22), addii(B11,B22));
33 33363 : GEN M2 = mulii(addii(A21,A22), B11);
34 33363 : GEN M3 = mulii(A11, subii(B12,B22));
35 33363 : GEN M4 = mulii(A22, subii(B21,B11));
36 33363 : GEN M5 = mulii(addii(A11,A12), B22);
37 33363 : GEN M6 = mulii(subii(A21,A11), addii(B11,B12));
38 33363 : GEN M7 = mulii(subii(A12,A22), addii(B21,B22));
39 33363 : GEN T1 = addii(M1,M4), T2 = subii(M7,M5);
40 33363 : GEN T3 = subii(M1,M2), T4 = addii(M3,M6);
41 33363 : retmkmat2(mkcol2(addii(T1,T2), addii(M2,M4)),
42 : mkcol2(addii(M3,M5), addii(T3,T4)));
43 : }
44 : }
45 :
46 : static GEN
47 17409 : matid2(void)
48 : {
49 17409 : 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 2906917 : mulq(GEN M, GEN q)
56 : {
57 2906917 : GEN u, v, res = cgetg(3, t_MAT);
58 2906923 : u = addii(mulii(gcoeff(M,1,1), q), gcoeff(M,1,2));
59 2906916 : v = addii(mulii(gcoeff(M,2,1), q), gcoeff(M,2,2));
60 2906917 : gel(res,1) = mkcol2(u, v);
61 2906927 : gel(res,2) = gel(M,1);
62 2906927 : return res;
63 : }
64 : static GEN
65 1072 : mulqab(GEN M, GEN q, GEN *ap, GEN *bp)
66 : {
67 1072 : GEN b = subii(*ap, mulii(*bp, q));
68 1072 : *ap = *bp; *bp = b;
69 1072 : return mulq(M,q);
70 : }
71 :
72 : /* Return M*[q,1;1,0]^-1 */
73 :
74 : static GEN
75 20168 : mulqi(GEN M, GEN q, GEN *ap, GEN *bp)
76 : {
77 : GEN u, v, res, a;
78 20168 : a = addii(mulii(*ap, q), *bp);
79 20168 : *bp = *ap; *ap = a;
80 20168 : res = cgetg(3, t_MAT);
81 20168 : u = subii(gcoeff(M,1,1),mulii(gcoeff(M,1,2), q));
82 20168 : v = subii(gcoeff(M,2,1),mulii(gcoeff(M,2,2), q));
83 20168 : gel(res,1) = gel(M,2);
84 20168 : gel(res,2) = mkcol2(u,v);
85 20168 : return res;
86 : }
87 :
88 : static long
89 5157176 : uexpi(GEN a)
90 5157176 : { return expi(subiu(shifti(a,1),1)); }
91 :
92 : static GEN
93 1574347 : FIXUP0(GEN M, GEN *a, GEN *b, long m)
94 : {
95 1574347 : long cnt=0;
96 2055168 : while (expi(*b) >= m)
97 : {
98 480821 : GEN r, q = dvmdii(*a, *b, &r);
99 480821 : *a = *b; *b = r;
100 480821 : M = mulq(M, q);
101 480821 : cnt++;
102 : };
103 1574347 : if (cnt>6) pari_err_BUG("FIXUP0");
104 1574347 : return M;
105 : }
106 :
107 : static long
108 2711698 : signdet(GEN Q)
109 : {
110 2711698 : long a = Mod4(gcoeff(Q,1,1)), b = Mod4(gcoeff(Q,1,2));
111 2711758 : long c = Mod4(gcoeff(Q,2,1)), d = Mod4(gcoeff(Q,2,2));
112 2711789 : return ((a*d-b*c)&3)==1 ? 1 : -1;
113 : }
114 :
115 : static GEN
116 1115567 : ZM_inv2(GEN M)
117 : {
118 1115567 : long e = signdet(M);
119 1676598 : if (e==1) return mkmat22(gcoeff(M,2,2),negi(gcoeff(M,1,2)),
120 560970 : negi(gcoeff(M,2,1)),gcoeff(M,1,1));
121 554668 : else return mkmat22(negi(gcoeff(M,2,2)),gcoeff(M,1,2),
122 554675 : gcoeff(M,2,1),negi(gcoeff(M,1,1)));
123 : }
124 :
125 : static GEN
126 20168 : lastq(GEN Q)
127 : {
128 20168 : GEN p = gcoeff(Q,1,1), q = gcoeff(Q,1,2), s = gcoeff(Q,2,2);
129 20168 : if (signe(q)==0) pari_err_BUG("halfgcd");
130 20168 : if (signe(s)==0) return p;
131 20168 : if (equali1(q)) return subiu(p,1);
132 20168 : return divii(p, q);
133 : }
134 :
135 : static GEN
136 39549 : mulT(GEN Q, GEN *ap, GEN *bp)
137 : {
138 39549 : *ap = addii(*ap, *bp);
139 39549 : *bp = negi(*bp);
140 39549 : return mkmat2(gel(Q,1),
141 39549 : mkcol2(subii(gcoeff(Q,1,1), gcoeff(Q,1,2))
142 39549 : , subii(gcoeff(Q,2,1), gcoeff(Q,2,2))));
143 : }
144 :
145 : static GEN
146 1596131 : FIXUP1(GEN M, GEN a, GEN b, long m, long t, GEN *ap, GEN *bp)
147 : {
148 1596131 : GEN Q = gel(M,1), a0 = gel(M,2), b0 = gel(M,3);
149 1596131 : GEN q, am = remi2n(a, m), bm = remi2n(b, m);
150 1596131 : if (signdet(Q)==-1)
151 : {
152 800198 : *ap = subii(mulii(bm, gcoeff(Q,1,2)),mulii(am, gcoeff(Q,2,2)));
153 800198 : *bp = subii(mulii(am, gcoeff(Q,2,1)),mulii(bm, gcoeff(Q,1,1)));
154 800198 : *ap = addii(*ap, shifti(addii(a0, gcoeff(Q,2,2)), m));
155 800198 : *bp = addii(*bp, shifti(subii(b0, gcoeff(Q,2,1)), m));
156 800198 : if (signe(*bp) >= 0)
157 759403 : return Q;
158 40795 : if (expi(addii(*ap,*bp)) >= m+t)
159 39549 : return mulT(Q, ap ,bp);
160 1246 : q = lastq(Q);
161 1246 : Q = mulqi(Q, q, ap, bp);
162 1246 : if (cmpiu(q, 2)>=0)
163 1072 : return mulqab(Q, subiu(q,1), ap, bp);
164 : else
165 174 : return mulqi(Q, lastq(Q), ap, bp);
166 : }
167 : else
168 : {
169 795933 : *ap = subii(mulii(am, gcoeff(Q,2,2)),mulii(bm, gcoeff(Q,1,2)));
170 795933 : *bp = subii(mulii(bm, gcoeff(Q,1,1)),mulii(am, gcoeff(Q,2,1)));
171 795933 : *ap = addii(*ap, shifti(subii(a0, gcoeff(Q,2,2)), m));
172 795933 : *bp = addii(*bp, shifti(addii(b0, gcoeff(Q,2,1)), m));
173 795933 : if (expi(*ap) >= m+t)
174 777185 : return FIXUP0(Q, ap, bp, m+t);
175 : else
176 18748 : return signe(gcoeff(Q,1,2))==0? Q: mulqi(Q, lastq(Q), ap, bp);
177 : }
178 : }
179 :
180 : static long
181 4360034 : magic_threshold(GEN a)
182 4360034 : { return (3+uexpi(a))>>1; }
183 :
184 : static GEN
185 1921924 : HGCD_basecase(GEN a, GEN b)
186 : {
187 1921924 : pari_sp av=avma;
188 1921924 : long m = magic_threshold(a);
189 : GEN u,u1,v,v1;
190 1921902 : u1 = v = gen_0;
191 1921902 : u = v1 = gen_1;
192 65324541 : while (expi(b) >= m)
193 : {
194 63403129 : GEN r, q = dvmdii(a,b, &r);
195 63402414 : a = b; b = r; swap(u,u1); swap(v,v1);
196 63402414 : u = addii(mulii(u1, q), u);
197 63402512 : v = addii(mulii(v1, q), v);
198 63402631 : if (gc_needed(av,2))
199 : {
200 0 : if (DEBUGMEM>1) pari_warn(warnmem,"halfgcd (d = %ld)",lgefint(b));
201 0 : gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
202 : }
203 : }
204 1921544 : return gerepilecopy(av, mkvec3(mkmat22(u,u1,v,v1), a, b));
205 : }
206 :
207 : static GEN HGCD(GEN x, GEN y);
208 :
209 : /*
210 : Based on
211 : Klaus Thull and Chee K. Yap,
212 : A unified approach to HGCD algorithms for polynomials andintegers,
213 : 1990, Manuscript.
214 : URL: http://cs.nyu.edu/cs/faculty/yap/papers.
215 : */
216 :
217 : static GEN
218 816396 : HGCD_split(GEN a, GEN b)
219 : {
220 816396 : pari_sp av = avma;
221 816396 : long m = magic_threshold(a), t, l, k, tp;
222 : GEN a0, b0, ap, bp, c, d, c0, d0, cp, dp, R, S, T, q, r;
223 816396 : if (signe(b) < 0 || cmpii(a,b)<0) pari_err_BUG("HGCD_split");
224 816396 : if (expi(b) < m)
225 17409 : return gerepilecopy(av, mkvec3(matid2(), a, b));
226 798987 : a0 = addiu(shifti(a, -m), 1);
227 798987 : if (cmpiu(a0,7) <= 0)
228 : {
229 0 : R = FIXUP0(matid2(), &a, &b, m);
230 0 : return gerepilecopy(av, mkvec3(R, a, b));
231 : }
232 798987 : b0 = shifti(b,-m);
233 798987 : t = magic_threshold(a0);
234 798987 : R = FIXUP1(HGCD(a0,b0),a, b, m, t, &ap, &bp);
235 798987 : if (expi(bp) < m)
236 1825 : return gerepilecopy(av, mkvec3(R, ap, bp));
237 797162 : q = dvmdii(ap, bp, &r);
238 797162 : c = bp; d = r;
239 797162 : if (cmpiu(shifti(c,-m),6) <= 0)
240 : {
241 18 : R = FIXUP0(mulq(R, q), &c, &d, m);
242 18 : return gerepilecopy(av, mkvec3(R, c, d));
243 : }
244 797144 : l = uexpi(c);
245 797144 : k = 2*m-l-1; if (k<0) pari_err_BUG("halfgcd");
246 797144 : c0 = addiu(shifti(c, -k), 1); if (cmpiu(c0,8)<0) pari_err_BUG("halfgcd");
247 797144 : d0 = shifti(d, -k);
248 797144 : tp = magic_threshold(c0);
249 797144 : S = FIXUP1(HGCD(c0,d0), c, d, k, tp, &cp, &dp);
250 797144 : if (!(expi(cp)>=m+1 && m+1 > expi(dp))) pari_err_BUG("halfgcd");
251 797144 : T = FIXUP0(ZM2_mul(mulq(R, q), S), &cp, &dp, m);
252 797144 : return gerepilecopy(av, mkvec3(T, cp, dp));
253 : }
254 :
255 : static GEN
256 2738323 : HGCD(GEN x, GEN y)
257 : {
258 2738323 : if (lgefint(y)-2 < HALFGCD_LIMIT)
259 1921927 : return HGCD_basecase(x, y);
260 : else
261 816396 : return HGCD_split(x, y);
262 : }
263 :
264 : static GEN
265 1175613 : HGCD0(GEN x, GEN y)
266 : {
267 1175613 : if (signe(y) >= 0 && cmpii(x, y) >= 0)
268 1142132 : return HGCD(x, y);
269 33495 : if (cmpii(x, y) < 0)
270 : {
271 22965 : GEN M = HGCD0(y, x), Q = gel(M,1);
272 22965 : return mkvec3(mkmat22(gcoeff(Q,2,1),gcoeff(Q,2,2),gcoeff(Q,1,1),gcoeff(Q,1,2)),
273 22965 : gel(M,2),gel(M,3));
274 : } /* Now y <= x*/
275 10530 : if (signe(x) <= 0)
276 : { /* y <= x <=0 */
277 60 : GEN M = HGCD(negi(y), negi(x)), Q = gel(M,1);
278 60 : return mkvec3(mkmat22(negi(gcoeff(Q,2,1)),negi(gcoeff(Q,2,2)),
279 60 : negi(gcoeff(Q,1,1)),negi(gcoeff(Q,1,2))),
280 60 : gel(M,2),gel(M,3));
281 : }
282 : else /* y <= 0 <=x */
283 : {
284 10470 : GEN M = HGCD0(x, negi(y)), Q = gel(M,1);
285 10470 : return mkvec3(mkmat22(gcoeff(Q,1,1),gcoeff(Q,1,2),negi(gcoeff(Q,2,1)),negi(gcoeff(Q,2,2))),
286 10470 : gel(M,2),gel(M,3));
287 : }
288 : }
289 :
290 : GEN
291 1114890 : halfgcdii(GEN A, GEN B)
292 : {
293 1114890 : pari_sp av = avma;
294 1114890 : GEN M, Q, a, b, m = absi_shallow(abscmpii(A, B)>0 ? A: B);
295 1114940 : M = HGCD0(A,B); Q = gel(M,1); a = gel(M,2); b = gel(M,3);
296 2742305 : while (signe(b) && cmpii(sqri(b), m) >= 0)
297 : {
298 1627353 : GEN r, q = dvmdii(a, b, &r);
299 1627318 : a = b; b = r;
300 1627318 : Q = mulq(Q, q);
301 : }
302 1114871 : return gerepilecopy(av, mkvec2(ZM_inv2(Q),mkcol2(a,b)));
303 : }
|
