Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : static int
19 1807294 : check_ZV(GEN x, long l)
20 : {
21 : long i;
22 10566528 : for (i=1; i<l; i++)
23 8759325 : if (typ(gel(x,i)) != t_INT) return 0;
24 1807203 : return 1;
25 : }
26 : void
27 1422200 : RgV_check_ZV(GEN A, const char *s)
28 : {
29 1422200 : if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
30 1422193 : }
31 : void
32 640276 : RgM_check_ZM(GEN A, const char *s)
33 : {
34 640276 : long n = lg(A);
35 640276 : if (n != 1)
36 : {
37 640129 : long j, m = lgcols(A);
38 2447331 : for (j=1; j<n; j++)
39 1807294 : if (!check_ZV(gel(A,j), m))
40 91 : pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
41 : }
42 640184 : }
43 :
44 : /* assume m > 1 */
45 : static long
46 109559742 : ZV_max_lg_i(GEN x, long m)
47 : {
48 109559742 : long i, l = lgefint(gel(x,1));
49 769470523 : for (i = 2; i < m; i++) l = maxss(l, lgefint(gel(x,i)));
50 109560961 : return l;
51 : }
52 : long
53 10640 : ZV_max_lg(GEN x)
54 : {
55 10640 : long m = lg(x);
56 10640 : return m == 1? 2: ZV_max_lg_i(x, m);
57 : }
58 :
59 : /* assume n > 1 and m > 1 */
60 : static long
61 28231315 : ZM_max_lg_i(GEN x, long n, long m)
62 : {
63 28231315 : long j, l = ZV_max_lg_i(gel(x,1), m);
64 109549687 : for (j = 2; j < n; j++) l = maxss(l, ZV_max_lg_i(gel(x,j), m));
65 28231558 : return l;
66 : }
67 : long
68 21647 : ZM_max_lg(GEN x)
69 : {
70 21647 : long n = lg(x), m;
71 21647 : if (n == 1) return 2;
72 21647 : m = lgcols(x); return m == 1? 2: ZM_max_lg_i(x, n, m);
73 : }
74 :
75 : /* assume m > 1 */
76 : static long
77 0 : ZV_max_expi_i(GEN x, long m)
78 : {
79 0 : long i, prec = expi(gel(x,1));
80 0 : for (i = 2; i < m; i++) prec = maxss(prec, expi(gel(x,i)));
81 0 : return prec;
82 : }
83 : long
84 0 : ZV_max_expi(GEN x)
85 : {
86 0 : long m = lg(x);
87 0 : return m == 1? 2: ZV_max_expi_i(x, m);
88 : }
89 :
90 : /* assume n > 1 and m > 1 */
91 : static long
92 0 : ZM_max_expi_i(GEN x, long n, long m)
93 : {
94 0 : long j, prec = ZV_max_expi_i(gel(x,1), m);
95 0 : for (j = 2; j < n; j++) prec = maxss(prec, ZV_max_expi_i(gel(x,j), m));
96 0 : return prec;
97 : }
98 : long
99 0 : ZM_max_expi(GEN x)
100 : {
101 0 : long n = lg(x), m;
102 0 : if (n == 1) return 2;
103 0 : m = lgcols(x); return m == 1? 2: ZM_max_expi_i(x, n, m);
104 : }
105 :
106 : GEN
107 4001 : ZM_supnorm(GEN x)
108 : {
109 4001 : long i, j, h, lx = lg(x);
110 4001 : GEN s = gen_0;
111 4001 : if (lx == 1) return gen_1;
112 4001 : h = lgcols(x);
113 24408 : for (j=1; j<lx; j++)
114 : {
115 20407 : GEN xj = gel(x,j);
116 281690 : for (i=1; i<h; i++)
117 : {
118 261283 : GEN c = gel(xj,i);
119 261283 : if (abscmpii(c, s) > 0) s = c;
120 : }
121 : }
122 4001 : return absi(s);
123 : }
124 :
125 : /********************************************************************/
126 : /** **/
127 : /** MULTIPLICATION **/
128 : /** **/
129 : /********************************************************************/
130 : /* x nonempty ZM, y a compatible nc (dimension > 0). */
131 : static GEN
132 0 : ZM_nc_mul_i(GEN x, GEN y, long c, long l)
133 : {
134 : long i, j;
135 : pari_sp av;
136 0 : GEN z = cgetg(l,t_COL), s;
137 :
138 0 : for (i=1; i<l; i++)
139 : {
140 0 : av = avma; s = muliu(gcoeff(x,i,1),y[1]);
141 0 : for (j=2; j<c; j++)
142 0 : if (y[j]) s = addii(s, muliu(gcoeff(x,i,j),y[j]));
143 0 : gel(z,i) = gerepileuptoint(av,s);
144 : }
145 0 : return z;
146 : }
147 :
148 : /* x ZV, y a compatible zc. */
149 : GEN
150 2229528 : ZV_zc_mul(GEN x, GEN y)
151 : {
152 2229528 : long j, l = lg(x);
153 2229528 : pari_sp av = avma;
154 2229528 : GEN s = mulis(gel(x,1),y[1]);
155 50292998 : for (j=2; j<l; j++)
156 48063470 : if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
157 2229528 : return gerepileuptoint(av,s);
158 : }
159 :
160 : /* x nonempty ZM, y a compatible zc (dimension > 0). */
161 : static GEN
162 20983483 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
163 : {
164 : long i, j;
165 20983483 : GEN z = cgetg(l,t_COL);
166 :
167 126937771 : for (i=1; i<l; i++)
168 : {
169 105959083 : pari_sp av = avma;
170 105959083 : GEN s = mulis(gcoeff(x,i,1),y[1]);
171 1178830286 : for (j=2; j<c; j++)
172 1072872498 : if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
173 105957788 : gel(z,i) = gerepileuptoint(av,s);
174 : }
175 20978688 : return z;
176 : }
177 : GEN
178 18989932 : ZM_zc_mul(GEN x, GEN y) {
179 18989932 : long l = lg(x);
180 18989932 : if (l == 1) return cgetg(1, t_COL);
181 18989932 : return ZM_zc_mul_i(x,y, l, lgcols(x));
182 : }
183 :
184 : /* y nonempty ZM, x a compatible zv (dimension > 0). */
185 : GEN
186 1736 : zv_ZM_mul(GEN x, GEN y) {
187 1736 : long i,j, lx = lg(x), ly = lg(y);
188 : GEN z;
189 1736 : if (lx == 1) return zerovec(ly-1);
190 1736 : z = cgetg(ly,t_VEC);
191 4046 : for (j=1; j<ly; j++)
192 : {
193 2310 : pari_sp av = avma;
194 2310 : GEN s = mulsi(x[1], gcoeff(y,1,j));
195 3990 : for (i=2; i<lx; i++)
196 1680 : if (x[i]) s = addii(s, mulsi(x[i], gcoeff(y,i,j)));
197 2310 : gel(z,j) = gerepileuptoint(av,s);
198 : }
199 1736 : return z;
200 : }
201 :
202 : /* x ZM, y a compatible zm (dimension > 0). */
203 : GEN
204 949041 : ZM_zm_mul(GEN x, GEN y)
205 : {
206 949041 : long j, c, l = lg(x), ly = lg(y);
207 949041 : GEN z = cgetg(ly, t_MAT);
208 949041 : if (l == 1) return z;
209 949034 : c = lgcols(x);
210 2942498 : for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
211 949030 : return z;
212 : }
213 : /* x ZM, y a compatible zn (dimension > 0). */
214 : GEN
215 0 : ZM_nm_mul(GEN x, GEN y)
216 : {
217 0 : long j, c, l = lg(x), ly = lg(y);
218 0 : GEN z = cgetg(ly, t_MAT);
219 0 : if (l == 1) return z;
220 0 : c = lgcols(x);
221 0 : for (j = 1; j < ly; j++) gel(z,j) = ZM_nc_mul_i(x, gel(y,j), l,c);
222 0 : return z;
223 : }
224 :
225 : /* Strassen-Winograd algorithm */
226 :
227 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
228 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
229 : static GEN
230 524320 : add_slices(long m, long n,
231 : GEN A, long ma, long da, long na, long ea,
232 : GEN B, long mb, long db, long nb, long eb)
233 : {
234 524320 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
235 524320 : GEN M = cgetg(n + 1, t_MAT), C;
236 :
237 3318102 : for (j = 1; j <= min_e; j++) {
238 2793782 : gel(M, j) = C = cgetg(m + 1, t_COL);
239 44972485 : for (i = 1; i <= min_d; i++)
240 42178703 : gel(C, i) = addii(gcoeff(A, ma + i, na + j),
241 42178703 : gcoeff(B, mb + i, nb + j));
242 2860151 : for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
243 2793782 : for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
244 2793782 : for (; i <= m; i++) gel(C, i) = gen_0;
245 : }
246 583716 : for (; j <= ea; j++) {
247 59396 : gel(M, j) = C = cgetg(m + 1, t_COL);
248 208387 : for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
249 59396 : for (; i <= m; i++) gel(C, i) = gen_0;
250 : }
251 524320 : for (; j <= eb; j++) {
252 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
253 0 : for (i = 1; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
254 0 : for (; i <= m; i++) gel(C, i) = gen_0;
255 : }
256 524320 : for (; j <= n; j++) gel(M, j) = zerocol(m);
257 524320 : return M;
258 : }
259 :
260 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
261 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
262 : static GEN
263 458780 : subtract_slices(long m, long n,
264 : GEN A, long ma, long da, long na, long ea,
265 : GEN B, long mb, long db, long nb, long eb)
266 : {
267 458780 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
268 458780 : GEN M = cgetg(n + 1, t_MAT), C;
269 :
270 2917513 : for (j = 1; j <= min_e; j++) {
271 2458733 : gel(M, j) = C = cgetg(m + 1, t_COL);
272 40899178 : for (i = 1; i <= min_d; i++)
273 38440445 : gel(C, i) = subii(gcoeff(A, ma + i, na + j),
274 38440445 : gcoeff(B, mb + i, nb + j));
275 2515834 : for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
276 2546786 : for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
277 2458733 : for (; i <= m; i++) gel(C, i) = gen_0;
278 : }
279 458780 : for (; j <= ea; j++) {
280 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
281 0 : for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
282 0 : for (; i <= m; i++) gel(C, i) = gen_0;
283 : }
284 493769 : for (; j <= eb; j++) {
285 34989 : gel(M, j) = C = cgetg(m + 1, t_COL);
286 127034 : for (i = 1; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
287 34989 : for (; i <= m; i++) gel(C, i) = gen_0;
288 : }
289 493769 : for (; j <= n; j++) gel(M, j) = zerocol(m);
290 458780 : return M;
291 : }
292 :
293 : static GEN ZM_mul_i(GEN x, GEN y, long l, long lx, long ly);
294 :
295 : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
296 : static GEN
297 65540 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
298 : {
299 65540 : pari_sp av = avma;
300 65540 : long m1 = (m + 1)/2, m2 = m/2,
301 65540 : n1 = (n + 1)/2, n2 = n/2,
302 65540 : p1 = (p + 1)/2, p2 = p/2;
303 : GEN A11, A12, A22, B11, B21, B22,
304 : S1, S2, S3, S4, T1, T2, T3, T4,
305 : M1, M2, M3, M4, M5, M6, M7,
306 : V1, V2, V3, C11, C12, C21, C22, C;
307 :
308 65540 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
309 65540 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
310 65540 : M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
311 65540 : if (gc_needed(av, 1))
312 0 : gerepileall(av, 2, &T2, &M2); /* destroy S1 */
313 65540 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
314 65540 : if (gc_needed(av, 1))
315 0 : gerepileall(av, 2, &M2, &T3); /* destroy T2 */
316 65540 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
317 65540 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
318 65540 : M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
319 65540 : if (gc_needed(av, 1))
320 0 : gerepileall(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
321 65540 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
322 65540 : if (gc_needed(av, 1))
323 0 : gerepileall(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
324 65540 : A11 = matslice(A, 1, m1, 1, n1);
325 65540 : B11 = matslice(B, 1, n1, 1, p1);
326 65540 : M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
327 65540 : if (gc_needed(av, 1))
328 0 : gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
329 65540 : A12 = matslice(A, 1, m1, n1 + 1, n);
330 65540 : B21 = matslice(B, n1 + 1, n, 1, p1);
331 65540 : M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
332 65540 : if (gc_needed(av, 1))
333 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
334 65540 : C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
335 65540 : if (gc_needed(av, 1))
336 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11); /* destroy M4 */
337 65540 : M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
338 65540 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
339 65540 : if (gc_needed(av, 1))
340 5 : gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4); /* destroy S3 */
341 65540 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
342 65540 : if (gc_needed(av, 1))
343 0 : gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4); /* destroy T3 */
344 65540 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
345 65540 : if (gc_needed(av, 1))
346 1 : gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1); /* destroy M1, M5 */
347 65540 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
348 65540 : M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
349 65540 : if (gc_needed(av, 1))
350 6 : gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6); /* destroy S4, B22 */
351 65540 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
352 65540 : M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
353 65540 : if (gc_needed(av, 1))
354 0 : gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7); /* destroy A22, T4 */
355 65540 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
356 65540 : C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
357 65540 : if (gc_needed(av, 1))
358 6 : gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12); /* destroy V3, M6 */
359 65540 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
360 65540 : if (gc_needed(av, 1))
361 0 : gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2); /* destroy V1, M2 */
362 65540 : C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
363 65540 : if (gc_needed(av, 1))
364 6 : gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21); /* destroy M7 */
365 65540 : C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
366 65540 : if (gc_needed(av, 1))
367 0 : gerepileall(av, 4, &C11, &C12, &C21, &C22); /* destroy V2, M3 */
368 65540 : C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
369 65540 : return gerepilecopy(av, C);
370 : }
371 :
372 : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
373 : static GEN
374 530010023 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
375 : {
376 530010023 : pari_sp av = avma;
377 530010023 : GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
378 : long k;
379 5300893709 : for (k = 2; k < lx; k++)
380 : {
381 4771188858 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
382 4769935258 : if (t != ZERO) c = addii(c, t);
383 : }
384 529704851 : return gerepileuptoint(av, c);
385 : }
386 : GEN
387 134539900 : ZMrow_ZC_mul(GEN x, GEN y, long i)
388 134539900 : { return ZMrow_ZC_mul_i(x, y, i, lg(x)); }
389 :
390 : /* return x * y, 1 < lx = lg(x), l = lgcols(x) */
391 : static GEN
392 74709536 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
393 : {
394 74709536 : GEN z = cgetg(l,t_COL);
395 : long i;
396 470196902 : for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
397 74674051 : return z;
398 : }
399 :
400 : static GEN
401 14050178 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
402 : {
403 : long j;
404 14050178 : GEN z = cgetg(ly, t_MAT);
405 65774522 : for (j = 1; j < ly; j++)
406 51727671 : gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
407 14046851 : return z;
408 : }
409 :
410 : static GEN
411 1119 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
412 : {
413 1119 : pari_sp av = avma;
414 1119 : long i, n = lg(P)-1;
415 : GEN H, T;
416 1119 : if (n == 1)
417 : {
418 0 : ulong p = uel(P,1);
419 0 : GEN a = ZM_to_Flm(A, p);
420 0 : GEN b = ZM_to_Flm(B, p);
421 0 : GEN Hp = gerepileupto(av, Flm_to_ZM(Flm_mul(a, b, p)));
422 0 : *mod = utoi(p); return Hp;
423 : }
424 1119 : T = ZV_producttree(P);
425 1120 : A = ZM_nv_mod_tree(A, P, T);
426 1120 : B = ZM_nv_mod_tree(B, P, T);
427 1120 : H = cgetg(n+1, t_VEC);
428 7938 : for(i=1; i <= n; i++)
429 6818 : gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
430 1120 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
431 1120 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
432 : }
433 :
434 : GEN
435 1120 : ZM_mul_worker(GEN P, GEN A, GEN B)
436 : {
437 1120 : GEN V = cgetg(3, t_VEC);
438 1119 : gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
439 1120 : return V;
440 : }
441 :
442 : static GEN
443 854 : ZM_mul_fast(GEN A, GEN B, long lA, long lB, long sA, long sB)
444 : {
445 854 : pari_sp av = avma;
446 : forprime_t S;
447 : GEN H, worker;
448 : long h;
449 854 : if (sA == 2 || sB == 2) return zeromat(nbrows(A),lB-1);
450 842 : h = 1 + (sA + sB - 4) * BITS_IN_LONG + expu(lA-1);
451 842 : init_modular_big(&S);
452 842 : worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
453 842 : H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
454 : nmV_chinese_center, FpM_center);
455 842 : return gerepileupto(av, H);
456 : }
457 :
458 : /* s = min(log_BIL |x|, log_BIL |y|), use Strassen-Winograd when
459 : * min(dims) > B */
460 : static long
461 14115929 : sw_bound(long s)
462 14115929 : { return s > 60 ? 2: s > 25 ? 4: s > 15 ? 8 : s > 8 ? 16 : 32; }
463 :
464 : /* assume lx > 1 and ly > 1; x is (l-1) x (lx-1), y is (lx-1) x (ly-1).
465 : * Return x * y */
466 : static GEN
467 21286176 : ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
468 : {
469 : long sx, sy, B;
470 : #ifdef LONG_IS_64BIT /* From Flm_mul_i */
471 18072254 : long Flm_sw_bound = 70;
472 : #else
473 3213922 : long Flm_sw_bound = 120;
474 : #endif
475 21286176 : if (l == 1) return zeromat(0, ly-1);
476 21284160 : if (lx==2 && l==2 && ly==2)
477 357290 : { retmkmat(mkcol(mulii(gcoeff(x,1,1), gcoeff(y,1,1)))); }
478 20926870 : if (lx==3 && l==3 && ly==3) return ZM2_mul(x, y);
479 14092551 : sx = ZM_max_lg_i(x, lx, l);
480 14093333 : sy = ZM_max_lg_i(y, ly, lx);
481 : /* Use modular reconstruction if Flm_mul would use Strassen and the input
482 : * sizes look balanced */
483 14093360 : if (lx > Flm_sw_bound && ly > Flm_sw_bound && l > Flm_sw_bound
484 866 : && sx <= 10 * sy && sy <= 10 * sx) return ZM_mul_fast(x,y, lx,ly, sx,sy);
485 :
486 14092506 : B = sw_bound(minss(sx, sy));
487 14092512 : if (l <= B || lx <= B || ly <= B)
488 14027002 : return ZM_mul_classical(x, y, l, lx, ly);
489 : else
490 65510 : return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
491 : }
492 :
493 : GEN
494 20963839 : ZM_mul(GEN x, GEN y)
495 : {
496 20963839 : long lx = lg(x), ly = lg(y);
497 20963839 : if (ly == 1) return cgetg(1,t_MAT);
498 20829047 : if (lx == 1) return zeromat(0, ly-1);
499 20827465 : return ZM_mul_i(x, y, lgcols(x), lx, ly);
500 : }
501 :
502 : GEN
503 548927 : QM_mul(GEN x, GEN y)
504 : {
505 548927 : GEN dx, nx = Q_primitive_part(x, &dx);
506 548927 : GEN dy, ny = Q_primitive_part(y, &dy);
507 548927 : GEN z = ZM_mul(nx, ny);
508 548927 : if (dx || dy)
509 : {
510 465316 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
511 465316 : if (!gequal1(d)) z = ZM_Q_mul(z, d);
512 : }
513 548927 : return z;
514 : }
515 :
516 : GEN
517 700 : QM_sqr(GEN x)
518 : {
519 700 : GEN dx, nx = Q_primitive_part(x, &dx);
520 700 : GEN z = ZM_sqr(nx);
521 700 : if (dx)
522 700 : z = ZM_Q_mul(z, gsqr(dx));
523 700 : return z;
524 : }
525 :
526 : GEN
527 142082 : QM_QC_mul(GEN x, GEN y)
528 : {
529 142082 : GEN dx, nx = Q_primitive_part(x, &dx);
530 142082 : GEN dy, ny = Q_primitive_part(y, &dy);
531 142082 : GEN z = ZM_ZC_mul(nx, ny);
532 142082 : if (dx || dy)
533 : {
534 142054 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
535 142054 : if (!gequal1(d)) z = ZC_Q_mul(z, d);
536 : }
537 142082 : return z;
538 : }
539 :
540 : /* assume result is symmetric */
541 : GEN
542 0 : ZM_multosym(GEN x, GEN y)
543 : {
544 0 : long j, lx, ly = lg(y);
545 : GEN M;
546 0 : if (ly == 1) return cgetg(1,t_MAT);
547 0 : lx = lg(x); /* = lgcols(y) */
548 0 : if (lx == 1) return cgetg(1,t_MAT);
549 : /* ly = lgcols(x) */
550 0 : M = cgetg(ly, t_MAT);
551 0 : for (j=1; j<ly; j++)
552 : {
553 0 : GEN z = cgetg(ly,t_COL), yj = gel(y,j);
554 : long i;
555 0 : for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
556 0 : for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
557 0 : gel(M,j) = z;
558 : }
559 0 : return M;
560 : }
561 :
562 : /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
563 : GEN
564 21 : ZM_mul_diag(GEN m, GEN d)
565 : {
566 : long j, l;
567 21 : GEN y = cgetg_copy(m, &l);
568 77 : for (j=1; j<l; j++)
569 : {
570 56 : GEN c = gel(d,j);
571 56 : gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
572 : }
573 21 : return y;
574 : }
575 : /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
576 : GEN
577 533483 : ZM_diag_mul(GEN d, GEN m)
578 : {
579 533483 : long i, j, l = lg(d), lm = lg(m);
580 533483 : GEN y = cgetg(lm, t_MAT);
581 1520102 : for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
582 1646673 : for (i=1; i<l; i++)
583 : {
584 1113408 : GEN c = gel(d,i);
585 1113408 : if (equali1(c))
586 236732 : for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
587 : else
588 3397013 : for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
589 : }
590 533265 : return y;
591 : }
592 :
593 : /* assume lx > 1 is lg(x) = lg(y) */
594 : static GEN
595 19128195 : ZV_dotproduct_i(GEN x, GEN y, long lx)
596 : {
597 19128195 : pari_sp av = avma;
598 19128195 : GEN c = mulii(gel(x,1), gel(y,1));
599 : long i;
600 143980152 : for (i = 2; i < lx; i++)
601 : {
602 124852856 : GEN t = mulii(gel(x,i), gel(y,i));
603 124852560 : if (t != gen_0) c = addii(c, t);
604 : }
605 19127296 : return gerepileuptoint(av, c);
606 : }
607 :
608 : /* x~ * y, assuming result is symmetric */
609 : GEN
610 527776 : ZM_transmultosym(GEN x, GEN y)
611 : {
612 527776 : long i, j, l, ly = lg(y);
613 : GEN M;
614 527776 : if (ly == 1) return cgetg(1,t_MAT);
615 : /* lg(x) = ly */
616 527776 : l = lgcols(y); /* = lgcols(x) */
617 527776 : M = cgetg(ly, t_MAT);
618 2691202 : for (i=1; i<ly; i++)
619 : {
620 2163426 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
621 2163426 : gel(M,i) = c;
622 7083212 : for (j=1; j<i; j++)
623 4919786 : gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
624 2163426 : gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
625 : }
626 527776 : return M;
627 : }
628 : /* x~ * y */
629 : GEN
630 2289 : ZM_transmul(GEN x, GEN y)
631 : {
632 2289 : long i, j, l, lx, ly = lg(y);
633 : GEN M;
634 2289 : if (ly == 1) return cgetg(1,t_MAT);
635 2289 : lx = lg(x);
636 2289 : l = lgcols(y);
637 2289 : if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
638 2289 : M = cgetg(ly, t_MAT);
639 6993 : for (i=1; i<ly; i++)
640 : {
641 4704 : GEN yi = gel(y,i), c = cgetg(lx,t_COL);
642 4704 : gel(M,i) = c;
643 12229 : for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
644 : }
645 2289 : return M;
646 : }
647 :
648 : /* assume l > 1; x is (l-1) x (l-1), return x^2.
649 : * FIXME: we ultimately rely on Strassen-Winograd which uses 7M + 15A.
650 : * Should use Bodrato's variant of Winograd, using 3M + 4S + 11A */
651 : static GEN
652 24192 : ZM_sqr_i(GEN x, long l)
653 : {
654 : long s;
655 24192 : if (l == 2) { retmkmat(mkcol(sqri(gcoeff(x,1,1)))); }
656 24192 : if (l == 3) return ZM2_sqr(x);
657 23422 : s = ZM_max_lg_i(x, l, l);
658 23422 : if (l > 70) return ZM_mul_fast(x, x, l, l, s, s);
659 23422 : if (l <= sw_bound(s))
660 23392 : return ZM_mul_classical(x, x, l, l, l);
661 : else
662 30 : return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
663 : }
664 :
665 : GEN
666 24192 : ZM_sqr(GEN x)
667 : {
668 24192 : long lx=lg(x);
669 24192 : if (lx==1) return cgetg(1,t_MAT);
670 24192 : return ZM_sqr_i(x, lx);
671 : }
672 : GEN
673 23070992 : ZM_ZC_mul(GEN x, GEN y)
674 : {
675 23070992 : long lx = lg(x);
676 23070992 : return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
677 : }
678 :
679 : GEN
680 4493177 : ZC_Z_div(GEN x, GEN c)
681 19443794 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
682 :
683 : GEN
684 15782 : ZM_Z_div(GEN x, GEN c)
685 177995 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
686 :
687 : GEN
688 2697405 : ZC_Q_mul(GEN A, GEN z)
689 : {
690 2697405 : pari_sp av = avma;
691 2697405 : long i, l = lg(A);
692 : GEN d, n, Ad, B, u;
693 2697405 : if (typ(z)==t_INT) return ZC_Z_mul(A,z);
694 2691147 : n = gel(z, 1); d = gel(z, 2);
695 2691147 : Ad = FpC_red(A, d);
696 2691098 : u = gcdii(d, FpV_factorback(Ad, NULL, d));
697 2691094 : B = cgetg(l, t_COL);
698 2691093 : if (equali1(u))
699 : {
700 420634 : for(i=1; i<l; i++)
701 354570 : gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
702 : } else
703 : {
704 18180929 : for(i=1; i<l; i++)
705 : {
706 15555974 : GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
707 15555919 : if (equalii(d, di))
708 10626603 : gel(B, i) = ni;
709 : else
710 4929281 : gel(B, i) = mkfrac(ni, diviiexact(d, di));
711 : }
712 : }
713 2691019 : return gerepilecopy(av, B);
714 : }
715 :
716 : GEN
717 1092252 : ZM_Q_mul(GEN x, GEN z)
718 : {
719 1092252 : if (typ(z)==t_INT) return ZM_Z_mul(x,z);
720 3215026 : pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
721 : }
722 :
723 : long
724 196399119 : zv_dotproduct(GEN x, GEN y)
725 : {
726 196399119 : long i, lx = lg(x);
727 : ulong c;
728 196399119 : if (lx == 1) return 0;
729 196399119 : c = uel(x,1)*uel(y,1);
730 3047359154 : for (i = 2; i < lx; i++)
731 2850960035 : c += uel(x,i)*uel(y,i);
732 196399119 : return c;
733 : }
734 :
735 : GEN
736 230629 : ZV_ZM_mul(GEN x, GEN y)
737 : {
738 230629 : long i, lx = lg(x), ly = lg(y);
739 : GEN z;
740 230629 : if (lx == 1) return zerovec(ly-1);
741 230510 : z = cgetg(ly, t_VEC);
742 882063 : for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
743 230510 : return z;
744 : }
745 :
746 : GEN
747 0 : ZC_ZV_mul(GEN x, GEN y)
748 : {
749 0 : long i,j, lx=lg(x), ly=lg(y);
750 : GEN z;
751 0 : if (ly==1) return cgetg(1,t_MAT);
752 0 : z = cgetg(ly,t_MAT);
753 0 : for (j=1; j < ly; j++)
754 : {
755 0 : gel(z,j) = cgetg(lx,t_COL);
756 0 : for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
757 : }
758 0 : return z;
759 : }
760 :
761 : GEN
762 6636290 : ZV_dotsquare(GEN x)
763 : {
764 : long i, lx;
765 : pari_sp av;
766 : GEN z;
767 6636290 : lx = lg(x);
768 6636290 : if (lx == 1) return gen_0;
769 6636290 : av = avma; z = sqri(gel(x,1));
770 26297705 : for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
771 6634988 : return gerepileuptoint(av,z);
772 : }
773 :
774 : GEN
775 16154257 : ZV_dotproduct(GEN x,GEN y)
776 : {
777 : long lx;
778 16154257 : if (x == y) return ZV_dotsquare(x);
779 11385176 : lx = lg(x);
780 11385176 : if (lx == 1) return gen_0;
781 11385176 : return ZV_dotproduct_i(x, y, lx);
782 : }
783 :
784 : static GEN
785 280 : _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
786 280 : { (void)data; return ZM_mul(x,y); }
787 : static GEN
788 23184 : _ZM_sqr(void *data /*ignored*/, GEN x)
789 23184 : { (void)data; return ZM_sqr(x); }
790 : /* FIXME: Using Bodrato's squaring, precomputations attached to fixed
791 : * multiplicand should be reused. And some postcomputations can be fused */
792 : GEN
793 0 : ZM_pow(GEN x, GEN n)
794 : {
795 0 : pari_sp av = avma;
796 0 : if (!signe(n)) return matid(lg(x)-1);
797 0 : return gerepilecopy(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
798 : }
799 : GEN
800 22638 : ZM_powu(GEN x, ulong n)
801 : {
802 22638 : pari_sp av = avma;
803 22638 : if (!n) return matid(lg(x)-1);
804 22638 : return gerepilecopy(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
805 : }
806 : /********************************************************************/
807 : /** **/
808 : /** ADD, SUB **/
809 : /** **/
810 : /********************************************************************/
811 : static GEN
812 34468943 : ZC_add_i(GEN x, GEN y, long lx)
813 : {
814 34468943 : GEN A = cgetg(lx, t_COL);
815 : long i;
816 442121183 : for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
817 34461910 : return A;
818 : }
819 : GEN
820 26387477 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
821 : GEN
822 363179 : ZC_Z_add(GEN x, GEN y)
823 : {
824 363179 : long k, lx = lg(x);
825 363179 : GEN z = cgetg(lx, t_COL);
826 363179 : if (lx == 1) pari_err_TYPE2("+",x,y);
827 363179 : gel(z,1) = addii(y,gel(x,1));
828 2500521 : for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
829 363182 : return z;
830 : }
831 :
832 : static GEN
833 8951504 : ZC_sub_i(GEN x, GEN y, long lx)
834 : {
835 : long i;
836 8951504 : GEN A = cgetg(lx, t_COL);
837 64020079 : for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
838 8951017 : return A;
839 : }
840 : GEN
841 8722259 : ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
842 : GEN
843 0 : ZC_Z_sub(GEN x, GEN y)
844 : {
845 0 : long k, lx = lg(x);
846 0 : GEN z = cgetg(lx, t_COL);
847 0 : if (lx == 1) pari_err_TYPE2("+",x,y);
848 0 : gel(z,1) = subii(gel(x,1), y);
849 0 : for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
850 0 : return z;
851 : }
852 : GEN
853 544242 : Z_ZC_sub(GEN a, GEN x)
854 : {
855 544242 : long k, lx = lg(x);
856 544242 : GEN z = cgetg(lx, t_COL);
857 544240 : if (lx == 1) pari_err_TYPE2("-",a,x);
858 544240 : gel(z,1) = subii(a, gel(x,1));
859 1508464 : for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
860 544228 : return z;
861 : }
862 :
863 : GEN
864 747540 : ZM_add(GEN x, GEN y)
865 : {
866 747540 : long lx = lg(x), l, j;
867 : GEN z;
868 747540 : if (lx == 1) return cgetg(1, t_MAT);
869 668202 : z = cgetg(lx, t_MAT); l = lgcols(x);
870 8749558 : for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
871 668208 : return z;
872 : }
873 : GEN
874 65151 : ZM_sub(GEN x, GEN y)
875 : {
876 65151 : long lx = lg(x), l, j;
877 : GEN z;
878 65151 : if (lx == 1) return cgetg(1, t_MAT);
879 65151 : z = cgetg(lx, t_MAT); l = lgcols(x);
880 294038 : for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
881 65151 : return z;
882 : }
883 : /********************************************************************/
884 : /** **/
885 : /** LINEAR COMBINATION **/
886 : /** **/
887 : /********************************************************************/
888 : /* return X/c assuming division is exact */
889 : GEN
890 4432178 : ZC_Z_divexact(GEN x, GEN c)
891 45854529 : { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
892 : GEN
893 2261 : ZC_divexactu(GEN x, ulong c)
894 11375 : { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
895 :
896 : GEN
897 429817 : ZM_Z_divexact(GEN x, GEN c)
898 2798570 : { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
899 :
900 : GEN
901 441 : ZM_divexactu(GEN x, ulong c)
902 2702 : { pari_APPLY_same(ZC_divexactu(gel(x,i), c)) }
903 :
904 : GEN
905 35109534 : ZC_Z_mul(GEN x, GEN c)
906 : {
907 35109534 : if (!signe(c)) return zerocol(lg(x)-1);
908 33753314 : if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
909 254340526 : pari_APPLY_type(t_COL, mulii(gel(x,i), c))
910 : }
911 :
912 : GEN
913 61732 : ZC_z_mul(GEN x, long c)
914 : {
915 61732 : if (!c) return zerocol(lg(x)-1);
916 53245 : if (c == 1) return ZC_copy(x);
917 48391 : if (c ==-1) return ZC_neg(x);
918 479003 : pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
919 : }
920 :
921 : GEN
922 28352 : zv_z_mul(GEN x, long n)
923 430576 : { pari_APPLY_long(x[i]*n) }
924 :
925 : /* return a ZM */
926 : GEN
927 0 : nm_Z_mul(GEN X, GEN c)
928 : {
929 0 : long i, j, h, l = lg(X), s = signe(c);
930 : GEN A;
931 0 : if (l == 1) return cgetg(1, t_MAT);
932 0 : h = lgcols(X);
933 0 : if (!s) return zeromat(h-1, l-1);
934 0 : if (is_pm1(c)) {
935 0 : if (s > 0) return Flm_to_ZM(X);
936 0 : X = Flm_to_ZM(X); ZM_togglesign(X); return X;
937 : }
938 0 : A = cgetg(l, t_MAT);
939 0 : for (j = 1; j < l; j++)
940 : {
941 0 : GEN a = cgetg(h, t_COL), x = gel(X, j);
942 0 : for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
943 0 : gel(A,j) = a;
944 : }
945 0 : return A;
946 : }
947 : GEN
948 2984697 : ZM_Z_mul(GEN X, GEN c)
949 : {
950 2984697 : long i, j, h, l = lg(X);
951 : GEN A;
952 2984697 : if (l == 1) return cgetg(1, t_MAT);
953 2984697 : h = lgcols(X);
954 2984687 : if (!signe(c)) return zeromat(h-1, l-1);
955 2939421 : if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
956 2654048 : A = cgetg(l, t_MAT);
957 13338530 : for (j = 1; j < l; j++)
958 : {
959 10684712 : GEN a = cgetg(h, t_COL), x = gel(X, j);
960 168378740 : for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
961 10684481 : gel(A,j) = a;
962 : }
963 2653818 : return A;
964 : }
965 : void
966 71952648 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
967 : {
968 : long i;
969 1564328809 : for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
970 71898557 : }
971 : /* X <- X + v Y (elementary col operation) */
972 : void
973 61215603 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
974 : {
975 61215603 : if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
976 : }
977 : void
978 29461390 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
979 : {
980 : long i;
981 29461390 : if (!v) return; /* v = 0 */
982 741665341 : for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
983 : }
984 :
985 : /* X + v Y, wasteful if (v = 0) */
986 : static GEN
987 15322893 : ZC_lincomb1(GEN v, GEN x, GEN y)
988 121655814 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
989 :
990 : /* -X + vY */
991 : static GEN
992 729303 : ZC_lincomb_1(GEN v, GEN x, GEN y)
993 4529026 : { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
994 :
995 : /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
996 : GEN
997 32858045 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
998 : {
999 : long su, sv;
1000 : GEN A;
1001 :
1002 32858045 : su = signe(u); if (!su) return ZC_Z_mul(Y, v);
1003 32857086 : sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
1004 32427890 : if (is_pm1(v))
1005 : {
1006 11023242 : if (is_pm1(u))
1007 : {
1008 9717842 : if (su != sv) A = ZC_sub(X, Y);
1009 2700693 : else A = ZC_add(X, Y);
1010 9717486 : if (su < 0) ZV_togglesign(A); /* in place but was created above */
1011 : }
1012 : else
1013 : {
1014 1305319 : if (sv > 0) A = ZC_lincomb1 (u, Y, X);
1015 600126 : else A = ZC_lincomb_1(u, Y, X);
1016 : }
1017 : }
1018 21405281 : else if (is_pm1(u))
1019 : {
1020 14746258 : if (su > 0) A = ZC_lincomb1 (v, X, Y);
1021 128503 : else A = ZC_lincomb_1(v, X, Y);
1022 : }
1023 : else
1024 : { /* not cgetg_copy: x may be a t_VEC */
1025 6662931 : long i, lx = lg(X);
1026 6662931 : A = cgetg(lx,t_COL);
1027 43590222 : for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
1028 : }
1029 32422151 : return A;
1030 : }
1031 :
1032 : /********************************************************************/
1033 : /** **/
1034 : /** CONVERSIONS **/
1035 : /** **/
1036 : /********************************************************************/
1037 : GEN
1038 400846 : ZV_to_nv(GEN x)
1039 749163 : { pari_APPLY_ulong(itou(gel(x,i))) }
1040 :
1041 : GEN
1042 130047 : zm_to_ZM(GEN x)
1043 815164 : { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
1044 :
1045 : GEN
1046 126 : zmV_to_ZMV(GEN x)
1047 791 : { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
1048 :
1049 : /* same as Flm_to_ZM but do not assume positivity */
1050 : GEN
1051 1022 : ZM_to_zm(GEN x)
1052 17472 : { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
1053 :
1054 : GEN
1055 366646 : zv_to_Flv(GEN x, ulong p)
1056 5418812 : { pari_APPLY_ulong(umodsu(x[i], p)) }
1057 :
1058 : GEN
1059 22694 : zm_to_Flm(GEN x, ulong p)
1060 351750 : { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
1061 :
1062 : GEN
1063 49 : ZMV_to_zmV(GEN x)
1064 399 : { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
1065 :
1066 : /********************************************************************/
1067 : /** **/
1068 : /** COPY, NEGATION **/
1069 : /** **/
1070 : /********************************************************************/
1071 : GEN
1072 13268627 : ZC_copy(GEN x)
1073 : {
1074 13268627 : long i, lx = lg(x);
1075 13268627 : GEN y = cgetg(lx, t_COL);
1076 83864562 : for (i=1; i<lx; i++)
1077 : {
1078 70595993 : GEN c = gel(x,i);
1079 70595993 : gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
1080 : }
1081 13268569 : return y;
1082 : }
1083 :
1084 : GEN
1085 509457 : ZM_copy(GEN x)
1086 4284537 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
1087 :
1088 : void
1089 342109 : ZV_neg_inplace(GEN M)
1090 : {
1091 342109 : long l = lg(M);
1092 1287764 : while (--l > 0) gel(M,l) = negi(gel(M,l));
1093 342146 : }
1094 : GEN
1095 6264502 : ZC_neg(GEN x)
1096 36497843 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
1097 :
1098 : GEN
1099 51070 : zv_neg(GEN x)
1100 655336 : { pari_APPLY_long(-x[i]) }
1101 : GEN
1102 126 : zv_neg_inplace(GEN M)
1103 : {
1104 126 : long l = lg(M);
1105 420 : while (--l > 0) M[l] = -M[l];
1106 126 : return M;
1107 : }
1108 : GEN
1109 77 : zv_abs(GEN x)
1110 5446 : { pari_APPLY_ulong(labs(x[i])) }
1111 : GEN
1112 1655047 : ZM_neg(GEN x)
1113 5039228 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
1114 :
1115 : void
1116 4919019 : ZV_togglesign(GEN M)
1117 : {
1118 4919019 : long l = lg(M);
1119 73284490 : while (--l > 0) togglesign_safe(&gel(M,l));
1120 4919108 : }
1121 : void
1122 0 : ZM_togglesign(GEN M)
1123 : {
1124 0 : long l = lg(M);
1125 0 : while (--l > 0) ZV_togglesign(gel(M,l));
1126 0 : }
1127 :
1128 : /********************************************************************/
1129 : /** **/
1130 : /** "DIVISION" mod HNF **/
1131 : /** **/
1132 : /********************************************************************/
1133 : /* Reduce ZC x modulo ZM y in HNF, may return x itself (not a copy) */
1134 : GEN
1135 10690741 : ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
1136 : {
1137 10690741 : long i, l = lg(x);
1138 : GEN q;
1139 :
1140 10690741 : if (Q) *Q = cgetg(l,t_COL);
1141 10690741 : if (l == 1) return cgetg(1,t_COL);
1142 61545811 : for (i = l-1; i>0; i--)
1143 : {
1144 50857721 : q = diviiround(gel(x,i), gcoeff(y,i,i));
1145 50856776 : if (signe(q)) {
1146 22936413 : togglesign(q);
1147 22936488 : x = ZC_lincomb(gen_1, q, x, gel(y,i));
1148 : }
1149 50855070 : if (Q) gel(*Q, i) = q;
1150 : }
1151 10688090 : return x;
1152 : }
1153 :
1154 : /* x = y Q + R, may return some columns of x (not copies) */
1155 : GEN
1156 453865 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
1157 : {
1158 453865 : long lx = lg(x), i;
1159 453865 : GEN R = cgetg(lx, t_MAT);
1160 453880 : if (Q)
1161 : {
1162 127353 : GEN q = cgetg(lx, t_MAT); *Q = q;
1163 188176 : for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
1164 : }
1165 : else
1166 823205 : for (i=1; i<lx; i++)
1167 : {
1168 496676 : pari_sp av = avma;
1169 496676 : GEN z = ZC_hnfrem(gel(x,i),y);
1170 496638 : gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
1171 : }
1172 453869 : return R;
1173 : }
1174 :
1175 : /********************************************************************/
1176 : /** **/
1177 : /** TESTS **/
1178 : /** **/
1179 : /********************************************************************/
1180 : int
1181 23103506 : zv_equal0(GEN V)
1182 : {
1183 23103506 : long l = lg(V);
1184 37579914 : while (--l > 0)
1185 30805753 : if (V[l]) return 0;
1186 6774161 : return 1;
1187 : }
1188 :
1189 : int
1190 14368470 : ZV_equal0(GEN V)
1191 : {
1192 14368470 : long l = lg(V);
1193 25401654 : while (--l > 0)
1194 24975475 : if (signe(gel(V,l))) return 0;
1195 426179 : return 1;
1196 : }
1197 : int
1198 16409 : ZMrow_equal0(GEN V, long i)
1199 : {
1200 16409 : long l = lg(V);
1201 25507 : while (--l > 0)
1202 21857 : if (signe(gcoeff(V,i,l))) return 0;
1203 3650 : return 1;
1204 : }
1205 :
1206 : static int
1207 6258198 : ZV_equal_lg(GEN V, GEN W, long l)
1208 : {
1209 25927897 : while (--l > 0)
1210 20151675 : if (!equalii(gel(V,l), gel(W,l))) return 0;
1211 5776222 : return 1;
1212 : }
1213 : int
1214 296304 : ZV_equal(GEN V, GEN W)
1215 : {
1216 296304 : long l = lg(V);
1217 296304 : if (lg(W) != l) return 0;
1218 296297 : return ZV_equal_lg(V, W, l);
1219 : }
1220 : int
1221 3347248 : ZM_equal(GEN A, GEN B)
1222 : {
1223 3347248 : long i, m, l = lg(A);
1224 3347248 : if (lg(B) != l) return 0;
1225 3347194 : if (l == 1) return 1;
1226 3347194 : m = lgcols(A);
1227 3347189 : if (lgcols(B) != m) return 0;
1228 9047969 : for (i = 1; i < l; i++)
1229 5961903 : if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
1230 3086066 : return 1;
1231 : }
1232 : int
1233 72400 : ZM_equal0(GEN A)
1234 : {
1235 72400 : long i, j, m, l = lg(A);
1236 72400 : if (l == 1) return 1;
1237 72400 : m = lgcols(A);
1238 123463 : for (j = 1; j < l; j++)
1239 2715086 : for (i = 1; i < m; i++)
1240 2664023 : if (signe(gcoeff(A,i,j))) return 0;
1241 15585 : return 1;
1242 : }
1243 : int
1244 30837951 : zv_equal(GEN V, GEN W)
1245 : {
1246 30837951 : long l = lg(V);
1247 30837951 : if (lg(W) != l) return 0;
1248 262401709 : while (--l > 0)
1249 232681637 : if (V[l] != W[l]) return 0;
1250 29720072 : return 1;
1251 : }
1252 :
1253 : int
1254 1638 : zvV_equal(GEN V, GEN W)
1255 : {
1256 1638 : long l = lg(V);
1257 1638 : if (lg(W) != l) return 0;
1258 80388 : while (--l > 0)
1259 79912 : if (!zv_equal(gel(V,l),gel(W,l))) return 0;
1260 476 : return 1;
1261 : }
1262 :
1263 : int
1264 759289 : ZM_ishnf(GEN x)
1265 : {
1266 759289 : long i,j, lx = lg(x);
1267 2282337 : for (i=1; i<lx; i++)
1268 : {
1269 1635784 : GEN xii = gcoeff(x,i,i);
1270 1635784 : if (signe(xii) <= 0) return 0;
1271 3176101 : for (j=1; j<i; j++)
1272 1577718 : if (signe(gcoeff(x,i,j))) return 0;
1273 3234607 : for (j=i+1; j<lx; j++)
1274 : {
1275 1711559 : GEN xij = gcoeff(x,i,j);
1276 1711559 : if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
1277 : }
1278 : }
1279 646553 : return 1;
1280 : }
1281 : int
1282 638057 : ZM_isidentity(GEN x)
1283 : {
1284 638057 : long i,j, lx = lg(x);
1285 :
1286 638057 : if (lx == 1) return 1;
1287 638050 : if (lx != lgcols(x)) return 0;
1288 3117070 : for (j=1; j<lx; j++)
1289 : {
1290 2481078 : GEN c = gel(x,j);
1291 7694945 : for (i=1; i<j; )
1292 5213903 : if (signe(gel(c,i++))) return 0;
1293 : /* i = j */
1294 2481042 : if (!equali1(gel(c,i++))) return 0;
1295 7701015 : for ( ; i<lx; )
1296 5221981 : if (signe(gel(c,i++))) return 0;
1297 : }
1298 635992 : return 1;
1299 : }
1300 : int
1301 556010 : ZM_isdiagonal(GEN x)
1302 : {
1303 556010 : long i,j, lx = lg(x);
1304 556010 : if (lx == 1) return 1;
1305 556010 : if (lx != lgcols(x)) return 0;
1306 :
1307 1437341 : for (j=1; j<lx; j++)
1308 : {
1309 1207881 : GEN c = gel(x,j);
1310 1695679 : for (i=1; i<j; i++)
1311 814307 : if (signe(gel(c,i))) return 0;
1312 2018906 : for (i++; i<lx; i++)
1313 1137569 : if (signe(gel(c,i))) return 0;
1314 : }
1315 229460 : return 1;
1316 : }
1317 : int
1318 103962 : ZM_isscalar(GEN x, GEN s)
1319 : {
1320 103962 : long i, j, lx = lg(x);
1321 :
1322 103962 : if (lx == 1) return 1;
1323 103962 : if (!s) s = gcoeff(x,1,1);
1324 103962 : if (equali1(s)) return ZM_isidentity(x);
1325 102751 : if (lx != lgcols(x)) return 0;
1326 547663 : for (j=1; j<lx; j++)
1327 : {
1328 456895 : GEN c = gel(x,j);
1329 2251257 : for (i=1; i<j; )
1330 1804356 : if (signe(gel(c,i++))) return 0;
1331 : /* i = j */
1332 446901 : if (!equalii(gel(c,i++), s)) return 0;
1333 2278861 : for ( ; i<lx; )
1334 1833949 : if (signe(gel(c,i++))) return 0;
1335 : }
1336 90768 : return 1;
1337 : }
1338 :
1339 : long
1340 106942 : ZC_is_ei(GEN x)
1341 : {
1342 106942 : long i, j = 0, l = lg(x);
1343 990877 : for (i = 1; i < l; i++)
1344 : {
1345 883936 : GEN c = gel(x,i);
1346 883936 : long s = signe(c);
1347 883936 : if (!s) continue;
1348 106936 : if (s < 0 || !is_pm1(c) || j) return 0;
1349 106935 : j = i;
1350 : }
1351 106941 : return j;
1352 : }
1353 :
1354 : /********************************************************************/
1355 : /** **/
1356 : /** MISCELLANEOUS **/
1357 : /** **/
1358 : /********************************************************************/
1359 : /* assume lg(x) = lg(y), x,y in Z^n */
1360 : int
1361 3132950 : ZV_cmp(GEN x, GEN y)
1362 : {
1363 3132950 : long fl,i, lx = lg(x);
1364 6344512 : for (i=1; i<lx; i++)
1365 5057434 : if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
1366 1287078 : return 0;
1367 : }
1368 : /* assume lg(x) = lg(y), x,y in Z^n */
1369 : int
1370 19740 : ZV_abscmp(GEN x, GEN y)
1371 : {
1372 19740 : long fl,i, lx = lg(x);
1373 54049 : for (i=1; i<lx; i++)
1374 53922 : if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
1375 127 : return 0;
1376 : }
1377 :
1378 : long
1379 20721470 : zv_content(GEN x)
1380 : {
1381 20721470 : long i, s, l = lg(x);
1382 20721470 : if (l == 1) return 0;
1383 20721463 : s = labs(x[1]);
1384 46522111 : for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
1385 20722422 : return s;
1386 : }
1387 : GEN
1388 297322 : ZV_content(GEN x)
1389 : {
1390 297322 : long i, l = lg(x);
1391 297322 : pari_sp av = avma;
1392 : GEN c;
1393 297322 : if (l == 1) return gen_0;
1394 297322 : if (l == 2) return absi(gel(x,1));
1395 204796 : c = gel(x,1);
1396 557875 : for (i = 2; i < l; i++)
1397 : {
1398 403575 : c = gcdii(c, gel(x,i));
1399 403575 : if (is_pm1(c)) { set_avma(av); return gen_1; }
1400 : }
1401 154300 : return gerepileuptoint(av, c);
1402 : }
1403 :
1404 : GEN
1405 3868301 : ZM_det_triangular(GEN mat)
1406 : {
1407 : pari_sp av;
1408 3868301 : long i,l = lg(mat);
1409 : GEN s;
1410 :
1411 3868301 : if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
1412 3465265 : av = avma; s = gcoeff(mat,1,1);
1413 9398645 : for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
1414 3464960 : return gerepileuptoint(av,s);
1415 : }
1416 :
1417 : /* assumes no overflow */
1418 : long
1419 949851 : zv_prod(GEN v)
1420 : {
1421 949851 : long n, i, l = lg(v);
1422 949851 : if (l == 1) return 1;
1423 959935 : n = v[1]; for (i = 2; i < l; i++) n *= v[i];
1424 771519 : return n;
1425 : }
1426 :
1427 : static GEN
1428 318095654 : _mulii(void *E, GEN a, GEN b)
1429 318095654 : { (void) E; return mulii(a, b); }
1430 :
1431 : /* product of ulongs */
1432 : GEN
1433 1864327 : zv_prod_Z(GEN v)
1434 : {
1435 : pari_sp av;
1436 1864327 : long k, m, n = lg(v)-1;
1437 1864327 : int stop = 0;
1438 : GEN V;
1439 1864327 : switch(n) {
1440 20902 : case 0: return gen_1;
1441 129878 : case 1: return utoi(v[1]);
1442 976692 : case 2: return muluu(v[1], v[2]);
1443 : }
1444 736855 : av = avma; m = n >> 1;
1445 736855 : V = cgetg(m + (odd(n)? 2: 1), t_VEC);
1446 151571307 : for (k = n; k; k--) /* start from the end: v is usually sorted */
1447 150835494 : if (v[k] & HIGHMASK) { stop = 1; break; }
1448 2422591 : while (!stop)
1449 : { /* HACK: handle V as a t_VECSMALL; gain a few iterations */
1450 81551667 : for (k = 1; k <= m; k++)
1451 : {
1452 79256797 : V[k] = uel(v,k<<1) * uel(v,(k<<1)-1);
1453 79256797 : if (V[k] & HIGHMASK) stop = 1; /* last "free" iteration */
1454 : }
1455 2294870 : if (odd(n))
1456 : {
1457 1352299 : if (n == 1) { set_avma(av); return utoi(v[1]); }
1458 743165 : V[++m] = v[n];
1459 : }
1460 1685736 : v = V; n = m; m = n >> 1;
1461 : }
1462 : /* n > 1; m > 0 */
1463 127721 : if (n == 2) { set_avma(av); return muluu(v[1], v[2]); }
1464 46319412 : for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
1465 87485 : if (odd(n)) gel(V, ++m) = utoipos(v[n]);
1466 87484 : setlg(V, m+1); /* HACK: now V is a bona fide t_VEC */
1467 87484 : return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
1468 : }
1469 : GEN
1470 14694393 : vecsmall_prod(GEN v)
1471 : {
1472 14694393 : pari_sp av = avma;
1473 14694393 : long k, m, n = lg(v)-1;
1474 : GEN V;
1475 14694393 : switch (n) {
1476 0 : case 0: return gen_1;
1477 0 : case 1: return stoi(v[1]);
1478 21 : case 2: return mulss(v[1], v[2]);
1479 : }
1480 14694372 : m = n >> 1;
1481 14694372 : V = cgetg(m + (odd(n)? 2: 1), t_VEC);
1482 161556906 : for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
1483 14694372 : if (odd(n)) gel(V,k) = stoi(v[n]);
1484 14694372 : return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
1485 : }
1486 :
1487 : GEN
1488 8819073 : ZV_prod(GEN v)
1489 : {
1490 8819073 : pari_sp av = avma;
1491 8819073 : long i, l = lg(v);
1492 : GEN n;
1493 8819073 : if (l == 1) return gen_1;
1494 8635125 : if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
1495 1292346 : n = gel(v,1);
1496 1292346 : if (l == 2) return icopy(n);
1497 2091108 : for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
1498 840374 : return gerepileuptoint(av, n);
1499 : }
1500 : /* assumes no overflow */
1501 : long
1502 16445 : zv_sum(GEN v)
1503 : {
1504 16445 : long n, i, l = lg(v);
1505 16445 : if (l == 1) return 0;
1506 95742 : n = v[1]; for (i = 2; i < l; i++) n += v[i];
1507 16424 : return n;
1508 : }
1509 : /* assumes no overflow and 0 <= n <= #v */
1510 : long
1511 0 : zv_sumpart(GEN v, long n)
1512 : {
1513 : long i, p;
1514 0 : if (!n) return 0;
1515 0 : p = v[1]; for (i = 2; i <= n; i++) p += v[i];
1516 0 : return p;
1517 : }
1518 : GEN
1519 77 : ZV_sum(GEN v)
1520 : {
1521 77 : pari_sp av = avma;
1522 77 : long i, l = lg(v);
1523 : GEN n;
1524 77 : if (l == 1) return gen_0;
1525 77 : n = gel(v,1);
1526 77 : if (l == 2) return icopy(n);
1527 581 : for (i = 2; i < l; i++) n = addii(n, gel(v,i));
1528 77 : return gerepileuptoint(av, n);
1529 : }
1530 :
1531 : /********************************************************************/
1532 : /** **/
1533 : /** GRAM SCHMIDT REDUCTION (integer matrices) **/
1534 : /** **/
1535 : /********************************************************************/
1536 :
1537 : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
1538 : static void
1539 362272 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
1540 : {
1541 362272 : long i, qq = itos_or_0(q);
1542 362274 : if (!qq)
1543 : {
1544 33377 : for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
1545 7069 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
1546 7069 : return;
1547 : }
1548 355205 : if (qq == 1) {
1549 148803 : for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
1550 102093 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
1551 253112 : } else if (qq == -1) {
1552 151483 : for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
1553 89199 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
1554 : } else {
1555 289616 : for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
1556 163934 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
1557 : }
1558 : }
1559 :
1560 : /* update L[k,] */
1561 : static void
1562 1033394 : ZRED(long k, long l, GEN x, GEN L, GEN B)
1563 : {
1564 1033394 : GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
1565 1033488 : if (!signe(q)) return;
1566 362259 : q = negi(q);
1567 362266 : Zupdate_row(k,l,q,L,B);
1568 362258 : gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
1569 : }
1570 :
1571 : /* Gram-Schmidt reduction, x a ZM */
1572 : static void
1573 1183946 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
1574 : {
1575 : long i, j;
1576 3780730 : for (j=1; j<=k; j++)
1577 : {
1578 2597475 : pari_sp av = avma;
1579 2597475 : GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
1580 5613363 : for (i=1; i<j; i++)
1581 : {
1582 3017180 : u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
1583 3016713 : u = diviiexact(u, gel(B,i));
1584 : }
1585 2596183 : gcoeff(L,k,j) = gerepileuptoint(av, u);
1586 : }
1587 1183255 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
1588 1183255 : }
1589 :
1590 : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
1591 : * Very inefficient if y is not LLL-reduced of maximal rank */
1592 : static GEN
1593 110487 : ZC_reducemodmatrix_i(GEN v, GEN y)
1594 : {
1595 110487 : GEN B, L, x = shallowconcat(y, v);
1596 110487 : long k, lx = lg(x), nx = lx-1;
1597 :
1598 110487 : B = scalarcol_shallow(gen_1, lx);
1599 110485 : L = zeromatcopy(nx, nx);
1600 454507 : for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
1601 344023 : for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1602 110485 : return gel(x,nx);
1603 : }
1604 : GEN
1605 110487 : ZC_reducemodmatrix(GEN v, GEN y) {
1606 110487 : pari_sp av = avma;
1607 110487 : return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
1608 : }
1609 :
1610 : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
1611 : * Very inefficient if y is not LLL-reduced of maximal rank */
1612 : static GEN
1613 226907 : ZM_reducemodmatrix_i(GEN v, GEN y)
1614 : {
1615 : GEN B, L, V;
1616 226907 : long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
1617 :
1618 226907 : V = cgetg(lv, t_MAT);
1619 226919 : B = scalarcol_shallow(gen_1, lx);
1620 226921 : L = zeromatcopy(nx, nx);
1621 600996 : for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
1622 692734 : for (j = 1; j < lg(v); j++)
1623 : {
1624 465852 : GEN x = shallowconcat(y, gel(v,j));
1625 465853 : ZincrementalGS(x, L, B, nx); /* overwrite last */
1626 1265730 : for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1627 465820 : gel(V,j) = gel(x,nx);
1628 : }
1629 226882 : return V;
1630 : }
1631 : GEN
1632 226907 : ZM_reducemodmatrix(GEN v, GEN y) {
1633 226907 : pari_sp av = avma;
1634 226907 : return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
1635 : }
1636 :
1637 : GEN
1638 98541 : ZC_reducemodlll(GEN x,GEN y)
1639 : {
1640 98541 : pari_sp av = avma;
1641 98541 : GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
1642 98545 : return gerepilecopy(av, z);
1643 : }
1644 : GEN
1645 0 : ZM_reducemodlll(GEN x,GEN y)
1646 : {
1647 0 : pari_sp av = avma;
1648 0 : GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
1649 0 : return gerepilecopy(av, z);
1650 : }
|
