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 1939155 : check_ZV(GEN x, long l)
20 : {
21 : long i;
22 13517710 : for (i=1; i<l; i++)
23 11578611 : if (typ(gel(x,i)) != t_INT) return 0;
24 1939099 : return 1;
25 : }
26 : void
27 1497225 : RgV_check_ZV(GEN A, const char *s)
28 : {
29 1497225 : if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
30 1497218 : }
31 : void
32 598635 : RgM_check_ZM(GEN A, const char *s)
33 : {
34 598635 : long n = lg(A);
35 598635 : if (n != 1)
36 : {
37 598502 : long j, m = lgcols(A);
38 2537601 : for (j=1; j<n; j++)
39 1939155 : if (!check_ZV(gel(A,j), m))
40 56 : pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
41 : }
42 598579 : }
43 :
44 : /* assume m > 1 */
45 : static long
46 114773290 : ZV_max_lg_i(GEN x, long m)
47 : {
48 114773290 : long i, l = lgefint(gel(x,1));
49 954965175 : for (i = 2; i < m; i++) l = maxss(l, lgefint(gel(x,i)));
50 114773290 : 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 27950360 : ZM_max_lg_i(GEN x, long n, long m)
62 : {
63 27950360 : long j, l = ZV_max_lg_i(gel(x,1), m);
64 114762650 : for (j = 2; j < n; j++) l = maxss(l, ZV_max_lg_i(gel(x,j), m));
65 27950360 : return l;
66 : }
67 : long
68 23338 : ZM_max_lg(GEN x)
69 : {
70 23338 : long n = lg(x), m;
71 23338 : if (n == 1) return 2;
72 23338 : 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 4260 : ZM_supnorm(GEN x)
108 : {
109 4260 : long i, j, h, lx = lg(x);
110 4260 : GEN s = gen_0;
111 4260 : if (lx == 1) return gen_1;
112 4260 : h = lgcols(x);
113 26221 : for (j=1; j<lx; j++)
114 : {
115 21961 : GEN xj = gel(x,j);
116 295088 : for (i=1; i<h; i++)
117 : {
118 273127 : GEN c = gel(xj,i);
119 273127 : if (abscmpii(c, s) > 0) s = c;
120 : }
121 : }
122 4260 : 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) = gc_INT(av,s);
144 : }
145 0 : return z;
146 : }
147 :
148 : /* x ZV, y a compatible zc. */
149 : GEN
150 2229556 : ZV_zc_mul(GEN x, GEN y)
151 : {
152 2229556 : long j, l = lg(x);
153 2229556 : pari_sp av = avma;
154 2229556 : GEN s = mulis(gel(x,1),y[1]);
155 50293138 : for (j=2; j<l; j++)
156 48063582 : if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
157 2229556 : return gc_INT(av,s);
158 : }
159 :
160 : /* x nonempty ZM, y a compatible zc (dimension > 0). */
161 : static GEN
162 20113387 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
163 : {
164 : long i, j;
165 20113387 : GEN z = cgetg(l,t_COL);
166 :
167 124642224 : for (i=1; i<l; i++)
168 : {
169 104528837 : pari_sp av = avma;
170 104528837 : GEN s = mulis(gcoeff(x,i,1),y[1]);
171 1182222756 : for (j=2; j<c; j++)
172 1077693919 : if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
173 104528837 : gel(z,i) = gc_INT(av,s);
174 : }
175 20113387 : return z;
176 : }
177 : GEN
178 17962095 : ZM_zc_mul(GEN x, GEN y) {
179 17962095 : long l = lg(x);
180 17962095 : if (l == 1) return cgetg(1, t_COL);
181 17962095 : 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 2408 : zv_ZM_mul(GEN x, GEN y) {
187 2408 : long i,j, lx = lg(x), ly = lg(y);
188 : GEN z;
189 2408 : if (lx == 1) return zerovec(ly-1);
190 2408 : z = cgetg(ly,t_VEC);
191 6734 : for (j=1; j<ly; j++)
192 : {
193 4326 : pari_sp av = avma;
194 4326 : GEN s = mulsi(x[1], gcoeff(y,1,j));
195 10038 : for (i=2; i<lx; i++)
196 5712 : if (x[i]) s = addii(s, mulsi(x[i], gcoeff(y,i,j)));
197 4326 : gel(z,j) = gc_INT(av,s);
198 : }
199 2408 : return z;
200 : }
201 :
202 : /* x ZM, y a compatible zm (dimension > 0). */
203 : GEN
204 1001399 : ZM_zm_mul(GEN x, GEN y)
205 : {
206 1001399 : long j, c, l = lg(x), ly = lg(y);
207 1001399 : GEN z = cgetg(ly, t_MAT);
208 1001399 : if (l == 1) return z;
209 1001392 : c = lgcols(x);
210 3152684 : for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
211 1001392 : 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 617072 : 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 617072 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
235 617072 : GEN M = cgetg(n + 1, t_MAT), C;
236 :
237 5195225 : for (j = 1; j <= min_e; j++) {
238 4578153 : gel(M, j) = C = cgetg(m + 1, t_COL);
239 89600305 : for (i = 1; i <= min_d; i++)
240 85022152 : gel(C, i) = addii(gcoeff(A, ma + i, na + j),
241 85022152 : gcoeff(B, mb + i, nb + j));
242 4645376 : for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
243 4578153 : for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
244 4578153 : for (; i <= m; i++) gel(C, i) = gen_0;
245 : }
246 678467 : for (; j <= ea; j++) {
247 61395 : gel(M, j) = C = cgetg(m + 1, t_COL);
248 227639 : for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
249 61395 : for (; i <= m; i++) gel(C, i) = gen_0;
250 : }
251 617072 : 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 617072 : for (; j <= n; j++) gel(M, j) = zerocol(m);
257 617072 : 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 539938 : 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 539938 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
268 539938 : GEN M = cgetg(n + 1, t_MAT), C;
269 :
270 4563530 : for (j = 1; j <= min_e; j++) {
271 4023592 : gel(M, j) = C = cgetg(m + 1, t_COL);
272 80069022 : for (i = 1; i <= min_d; i++)
273 76045430 : gel(C, i) = subii(gcoeff(A, ma + i, na + j),
274 76045430 : gcoeff(B, mb + i, nb + j));
275 4082765 : for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
276 4112702 : for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
277 4023592 : for (; i <= m; i++) gel(C, i) = gen_0;
278 : }
279 539938 : 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 575739 : for (; j <= eb; j++) {
285 35801 : gel(M, j) = C = cgetg(m + 1, t_COL);
286 130300 : for (i = 1; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
287 35801 : for (; i <= m; i++) gel(C, i) = gen_0;
288 : }
289 575739 : for (; j <= n; j++) gel(M, j) = zerocol(m);
290 539938 : 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 77134 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
298 : {
299 77134 : pari_sp av = avma;
300 77134 : long m1 = (m + 1)/2, m2 = m/2,
301 77134 : n1 = (n + 1)/2, n2 = n/2,
302 77134 : 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 77134 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
309 77134 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
310 77134 : M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
311 77134 : if (gc_needed(av, 1))
312 0 : (void)gc_all(av, 2, &T2, &M2); /* destroy S1 */
313 77134 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
314 77134 : if (gc_needed(av, 1))
315 0 : (void)gc_all(av, 2, &M2, &T3); /* destroy T2 */
316 77134 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
317 77134 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
318 77134 : M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
319 77134 : if (gc_needed(av, 1))
320 0 : (void)gc_all(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
321 77134 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
322 77134 : if (gc_needed(av, 1))
323 0 : (void)gc_all(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
324 77134 : A11 = matslice(A, 1, m1, 1, n1);
325 77134 : B11 = matslice(B, 1, n1, 1, p1);
326 77134 : M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
327 77134 : if (gc_needed(av, 1))
328 0 : (void)gc_all(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
329 77134 : A12 = matslice(A, 1, m1, n1 + 1, n);
330 77134 : B21 = matslice(B, n1 + 1, n, 1, p1);
331 77134 : M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
332 77134 : if (gc_needed(av, 1))
333 0 : (void)gc_all(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
334 77134 : C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
335 77134 : if (gc_needed(av, 1))
336 0 : (void)gc_all(av, 6, &M2, &T3, &M3, &S3, &M1, &C11); /* destroy M4 */
337 77134 : M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
338 77134 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
339 77134 : if (gc_needed(av, 1))
340 5 : (void)gc_all(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4); /* destroy S3 */
341 77134 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
342 77134 : if (gc_needed(av, 1))
343 0 : (void)gc_all(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4); /* destroy T3 */
344 77134 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
345 77134 : if (gc_needed(av, 1))
346 1 : (void)gc_all(av, 6, &M2, &M3, &S4, &T4, &C11, &V1); /* destroy M1, M5 */
347 77134 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
348 77134 : M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
349 77134 : if (gc_needed(av, 1))
350 6 : (void)gc_all(av, 6, &M2, &M3, &T4, &C11, &V1, &M6); /* destroy S4, B22 */
351 77134 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
352 77134 : M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
353 77134 : if (gc_needed(av, 1))
354 0 : (void)gc_all(av, 6, &M2, &M3, &C11, &V1, &M6, &M7); /* destroy A22, T4 */
355 77134 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
356 77134 : C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
357 77134 : if (gc_needed(av, 1))
358 6 : (void)gc_all(av, 6, &M2, &M3, &C11, &V1, &M7, &C12); /* destroy V3, M6 */
359 77134 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
360 77134 : if (gc_needed(av, 1))
361 0 : (void)gc_all(av, 5, &M3, &C11, &M7, &C12, &V2); /* destroy V1, M2 */
362 77134 : C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
363 77134 : if (gc_needed(av, 1))
364 6 : (void)gc_all(av, 5, &M3, &C11, &C12, &V2, &C21); /* destroy M7 */
365 77134 : C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
366 77134 : if (gc_needed(av, 1))
367 0 : (void)gc_all(av, 4, &C11, &C12, &C21, &C22); /* destroy V2, M3 */
368 77134 : C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
369 77134 : return gc_GEN(av, C);
370 : }
371 :
372 : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
373 : static GEN
374 609105018 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
375 : {
376 609105018 : pari_sp av = avma;
377 609105018 : GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
378 : long k;
379 7076828504 : for (k = 2; k < lx; k++)
380 : {
381 6467723486 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
382 6467723486 : if (t != ZERO) c = addii(c, t);
383 : }
384 609105018 : return gc_INT(av, c);
385 : }
386 : GEN
387 137967296 : ZMrow_ZC_mul(GEN x, GEN y, long i)
388 137967296 : { 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 79223357 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
393 : {
394 79223357 : GEN z = cgetg(l,t_COL);
395 : long i;
396 550361079 : for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
397 79223357 : return z;
398 : }
399 :
400 : static GEN
401 13897760 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
402 : {
403 : long j;
404 13897760 : GEN z = cgetg(ly, t_MAT);
405 67351597 : for (j = 1; j < ly; j++)
406 53453837 : gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
407 13897760 : return z;
408 : }
409 :
410 : static GEN
411 1105 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
412 : {
413 1105 : pari_sp av = avma;
414 1105 : long i, n = lg(P)-1;
415 : GEN H, T;
416 1105 : 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 = gc_upto(av, Flm_to_ZM(Flm_mul(a, b, p)));
422 0 : *mod = utoi(p); return Hp;
423 : }
424 1105 : T = ZV_producttree(P);
425 1105 : A = ZM_nv_mod_tree(A, P, T);
426 1105 : B = ZM_nv_mod_tree(B, P, T);
427 1105 : H = cgetg(n+1, t_VEC);
428 7878 : for(i=1; i <= n; i++)
429 6773 : gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
430 1105 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
431 1105 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
432 : }
433 :
434 : GEN
435 1105 : ZM_mul_worker(GEN P, GEN A, GEN B)
436 : {
437 1105 : GEN V = cgetg(3, t_VEC);
438 1105 : gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
439 1105 : return V;
440 : }
441 :
442 : static GEN
443 839 : ZM_mul_fast(GEN A, GEN B, long lA, long lB, long sA, long sB)
444 : {
445 839 : pari_sp av = avma;
446 : forprime_t S;
447 : GEN H, worker;
448 : long h;
449 839 : if (sA == 2 || sB == 2) return zeromat(nbrows(A),lB-1);
450 827 : h = 1 + (sA + sB - 4) * BITS_IN_LONG + expu(lA-1);
451 827 : init_modular_big(&S);
452 827 : worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
453 827 : H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
454 : nmV_chinese_center, FpM_center);
455 827 : return gc_upto(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 13974894 : sw_bound(long s)
462 13974894 : { 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 21531066 : 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 18283294 : long Flm_sw_bound = 70;
472 : #else
473 3247772 : long Flm_sw_bound = 120;
474 : #endif
475 21531066 : if (l == 1) return zeromat(0, ly-1);
476 21529166 : if (lx==2 && l==2 && ly==2)
477 374650 : { retmkmat(mkcol(mulii(gcoeff(x,1,1), gcoeff(y,1,1)))); }
478 21154516 : if (lx==3 && l==3 && ly==3) return ZM2_mul(x, y);
479 13951289 : sx = ZM_max_lg_i(x, lx, l);
480 13951289 : 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 13951289 : if (lx > Flm_sw_bound && ly > Flm_sw_bound && l > Flm_sw_bound
484 851 : && sx <= 10 * sy && sy <= 10 * sx) return ZM_mul_fast(x,y, lx,ly, sx,sy);
485 :
486 13950450 : B = sw_bound(minss(sx, sy));
487 13950450 : if (l <= B || lx <= B || ly <= B)
488 13873346 : return ZM_mul_classical(x, y, l, lx, ly);
489 : else
490 77104 : return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
491 : }
492 :
493 : GEN
494 21128776 : ZM_mul(GEN x, GEN y)
495 : {
496 21128776 : long lx = lg(x), ly = lg(y);
497 21128776 : if (ly == 1) return cgetg(1,t_MAT);
498 20992731 : if (lx == 1) return zeromat(0, ly-1);
499 20991128 : return ZM_mul_i(x, y, lgcols(x), lx, ly);
500 : }
501 :
502 : static GEN
503 0 : ZM_sqr_slice(GEN A, GEN P, GEN *mod)
504 : {
505 0 : pari_sp av = avma;
506 0 : long i, n = lg(P)-1;
507 : GEN H, T;
508 0 : if (n == 1)
509 : {
510 0 : ulong p = uel(P,1);
511 0 : GEN a = ZM_to_Flm(A, p);
512 0 : GEN Hp = gc_upto(av, Flm_to_ZM(Flm_sqr(a, p)));
513 0 : *mod = utoi(p); return Hp;
514 : }
515 0 : T = ZV_producttree(P);
516 0 : A = ZM_nv_mod_tree(A, P, T);
517 0 : H = cgetg(n+1, t_VEC);
518 0 : for(i=1; i <= n; i++)
519 0 : gel(H,i) = Flm_sqr(gel(A,i), P[i]);
520 0 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
521 0 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
522 : }
523 :
524 : GEN
525 0 : ZM_sqr_worker(GEN P, GEN A)
526 : {
527 0 : GEN V = cgetg(3, t_VEC);
528 0 : gel(V,1) = ZM_sqr_slice(A, P, &gel(V,2));
529 0 : return V;
530 : }
531 :
532 : static GEN
533 0 : ZM_sqr_fast(GEN A, long l, long s)
534 : {
535 0 : pari_sp av = avma;
536 : forprime_t S;
537 : GEN H, worker;
538 : long h;
539 0 : if (s == 2) return zeromat(l-1,l-1);
540 0 : h = 1 + (2*s - 4) * BITS_IN_LONG + expu(l-1);
541 0 : init_modular_big(&S);
542 0 : worker = snm_closure(is_entry("_ZM_sqr_worker"), mkvec(A));
543 0 : H = gen_crt("ZM_sqr", worker, &S, NULL, h, 0, NULL,
544 : nmV_chinese_center, FpM_center);
545 0 : return gc_upto(av, H);
546 : }
547 :
548 : GEN
549 559203 : QM_mul(GEN x, GEN y)
550 : {
551 559203 : GEN dx, nx = Q_primitive_part(x, &dx);
552 559203 : GEN dy, ny = Q_primitive_part(y, &dy);
553 559203 : GEN z = ZM_mul(nx, ny);
554 559203 : if (dx || dy)
555 : {
556 475088 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
557 475088 : if (!gequal1(d)) z = ZM_Q_mul(z, d);
558 : }
559 559203 : return z;
560 : }
561 :
562 : GEN
563 700 : QM_sqr(GEN x)
564 : {
565 700 : GEN dx, nx = Q_primitive_part(x, &dx);
566 700 : GEN z = ZM_sqr(nx);
567 700 : if (dx)
568 700 : z = ZM_Q_mul(z, gsqr(dx));
569 700 : return z;
570 : }
571 :
572 : GEN
573 127093 : QM_QC_mul(GEN x, GEN y)
574 : {
575 127093 : GEN dx, nx = Q_primitive_part(x, &dx);
576 127093 : GEN dy, ny = Q_primitive_part(y, &dy);
577 127093 : GEN z = ZM_ZC_mul(nx, ny);
578 127093 : if (dx || dy)
579 : {
580 127072 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
581 127072 : if (!gequal1(d)) z = ZC_Q_mul(z, d);
582 : }
583 127093 : return z;
584 : }
585 :
586 : /* assume result is symmetric */
587 : GEN
588 0 : ZM_multosym(GEN x, GEN y)
589 : {
590 0 : long j, lx, ly = lg(y);
591 : GEN M;
592 0 : if (ly == 1) return cgetg(1,t_MAT);
593 0 : lx = lg(x); /* = lgcols(y) */
594 0 : if (lx == 1) return cgetg(1,t_MAT);
595 : /* ly = lgcols(x) */
596 0 : M = cgetg(ly, t_MAT);
597 0 : for (j=1; j<ly; j++)
598 : {
599 0 : GEN z = cgetg(ly,t_COL), yj = gel(y,j);
600 : long i;
601 0 : for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
602 0 : for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
603 0 : gel(M,j) = z;
604 : }
605 0 : return M;
606 : }
607 :
608 : /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
609 : GEN
610 21 : ZM_mul_diag(GEN m, GEN d)
611 : {
612 : long j, l;
613 21 : GEN y = cgetg_copy(m, &l);
614 77 : for (j=1; j<l; j++)
615 : {
616 56 : GEN c = gel(d,j);
617 56 : gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
618 : }
619 21 : return y;
620 : }
621 : /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
622 : GEN
623 593153 : ZM_diag_mul(GEN d, GEN m)
624 : {
625 593153 : long i, j, l = lg(d), lm = lg(m);
626 593153 : GEN y = cgetg(lm, t_MAT);
627 1704469 : for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
628 1850415 : for (i=1; i<l; i++)
629 : {
630 1257262 : GEN c = gel(d,i);
631 1257262 : if (equali1(c))
632 274688 : for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
633 : else
634 3964063 : for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
635 : }
636 593153 : return y;
637 : }
638 :
639 : /* assume lx > 1 is lg(x) = lg(y) */
640 : static GEN
641 19537963 : ZV_dotproduct_i(GEN x, GEN y, long lx)
642 : {
643 19537963 : pari_sp av = avma;
644 19537963 : GEN c = mulii(gel(x,1), gel(y,1));
645 : long i;
646 148538791 : for (i = 2; i < lx; i++)
647 : {
648 129000828 : GEN t = mulii(gel(x,i), gel(y,i));
649 129000828 : if (t != gen_0) c = addii(c, t);
650 : }
651 19537963 : return gc_INT(av, c);
652 : }
653 :
654 : /* x~ * y, assuming result is symmetric */
655 : GEN
656 532275 : ZM_transmultosym(GEN x, GEN y)
657 : {
658 532275 : long i, j, l, ly = lg(y);
659 : GEN M;
660 532275 : if (ly == 1) return cgetg(1,t_MAT);
661 : /* lg(x) = ly */
662 532275 : l = lgcols(y); /* = lgcols(x) */
663 532275 : M = cgetg(ly, t_MAT);
664 2708943 : for (i=1; i<ly; i++)
665 : {
666 2176668 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
667 2176668 : gel(M,i) = c;
668 7109315 : for (j=1; j<i; j++)
669 4932647 : gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
670 2176668 : gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
671 : }
672 532275 : return M;
673 : }
674 : /* x~ * y */
675 : GEN
676 2289 : ZM_transmul(GEN x, GEN y)
677 : {
678 2289 : long i, j, l, lx, ly = lg(y);
679 : GEN M;
680 2289 : if (ly == 1) return cgetg(1,t_MAT);
681 2289 : lx = lg(x);
682 2289 : l = lgcols(y);
683 2289 : if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
684 2289 : M = cgetg(ly, t_MAT);
685 6993 : for (i=1; i<ly; i++)
686 : {
687 4704 : GEN yi = gel(y,i), c = cgetg(lx,t_COL);
688 4704 : gel(M,i) = c;
689 12229 : for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
690 : }
691 2289 : return M;
692 : }
693 :
694 : /* assume l > 1; x is (l-1) x (l-1), return x^2.
695 : * FIXME: we ultimately rely on Strassen-Winograd which uses 7M + 15A.
696 : * Should use Bodrato's variant of Winograd, using 3M + 4S + 11A */
697 : static GEN
698 25221 : ZM_sqr_i(GEN x, long l)
699 : {
700 : long s;
701 25221 : if (l == 2) { retmkmat(mkcol(sqri(gcoeff(x,1,1)))); }
702 25221 : if (l == 3) return ZM2_sqr(x);
703 24444 : s = ZM_max_lg_i(x, l, l);
704 24444 : if (l > 70) return ZM_sqr_fast(x, l, s);
705 24444 : if (l <= sw_bound(s))
706 24414 : return ZM_mul_classical(x, x, l, l, l);
707 : else
708 30 : return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
709 : }
710 :
711 : GEN
712 25221 : ZM_sqr(GEN x)
713 : {
714 25221 : long lx=lg(x);
715 25221 : if (lx==1) return cgetg(1,t_MAT);
716 25221 : return ZM_sqr_i(x, lx);
717 : }
718 : GEN
719 25859610 : ZM_ZC_mul(GEN x, GEN y)
720 : {
721 25859610 : long lx = lg(x);
722 25859610 : return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
723 : }
724 :
725 : GEN
726 3676481 : ZC_Z_div(GEN x, GEN c)
727 17415866 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
728 :
729 : GEN
730 19308 : ZM_Z_div(GEN x, GEN c)
731 221027 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
732 :
733 : GEN
734 2792946 : ZC_Q_mul(GEN A, GEN z)
735 : {
736 2792946 : pari_sp av = avma;
737 2792946 : long i, l = lg(A);
738 : GEN d, n, Ad, B, u;
739 2792946 : if (typ(z)==t_INT) return ZC_Z_mul(A,z);
740 2788557 : n = gel(z, 1); d = gel(z, 2);
741 2788557 : Ad = FpC_red(A, d);
742 2788557 : u = gcdii(d, FpV_factorback(Ad, NULL, d));
743 2788557 : B = cgetg(l, t_COL);
744 2788557 : if (equali1(u))
745 : {
746 414857 : for(i=1; i<l; i++)
747 350172 : gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
748 : } else
749 : {
750 18942502 : for(i=1; i<l; i++)
751 : {
752 16218630 : GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
753 16218630 : if (equalii(d, di))
754 11206206 : gel(B, i) = ni;
755 : else
756 5012424 : gel(B, i) = mkfrac(ni, diviiexact(d, di));
757 : }
758 : }
759 2788557 : return gc_GEN(av, B);
760 : }
761 :
762 : GEN
763 1166404 : ZM_Q_mul(GEN x, GEN z)
764 : {
765 1166404 : if (typ(z)==t_INT) return ZM_Z_mul(x,z);
766 3359715 : pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
767 : }
768 :
769 : long
770 196399119 : zv_dotproduct(GEN x, GEN y)
771 : {
772 196399119 : long i, lx = lg(x);
773 : ulong c;
774 196399119 : if (lx == 1) return 0;
775 196399119 : c = uel(x,1)*uel(y,1);
776 3047359154 : for (i = 2; i < lx; i++)
777 2850960035 : c += uel(x,i)*uel(y,i);
778 196399119 : return c;
779 : }
780 :
781 : GEN
782 231351 : ZV_ZM_mul(GEN x, GEN y)
783 : {
784 231351 : long i, lx = lg(x), ly = lg(y);
785 : GEN z;
786 231351 : if (lx == 1) return zerovec(ly-1);
787 231232 : z = cgetg(ly, t_VEC);
788 883969 : for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
789 231232 : return z;
790 : }
791 :
792 : GEN
793 0 : ZC_ZV_mul(GEN x, GEN y)
794 : {
795 0 : long i,j, lx=lg(x), ly=lg(y);
796 : GEN z;
797 0 : if (ly==1) return cgetg(1,t_MAT);
798 0 : z = cgetg(ly,t_MAT);
799 0 : for (j=1; j < ly; j++)
800 : {
801 0 : gel(z,j) = cgetg(lx,t_COL);
802 0 : for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
803 : }
804 0 : return z;
805 : }
806 :
807 : GEN
808 6869396 : ZV_dotsquare(GEN x)
809 : {
810 : long i, lx;
811 : pari_sp av;
812 : GEN z;
813 6869396 : lx = lg(x);
814 6869396 : if (lx == 1) return gen_0;
815 6869396 : av = avma; z = sqri(gel(x,1));
816 27055201 : for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
817 6869396 : return gc_INT(av,z);
818 : }
819 :
820 : GEN
821 16715602 : ZV_dotproduct(GEN x,GEN y)
822 : {
823 : long lx;
824 16715602 : if (x == y) return ZV_dotsquare(x);
825 11768386 : lx = lg(x);
826 11768386 : if (lx == 1) return gen_0;
827 11768386 : return ZV_dotproduct_i(x, y, lx);
828 : }
829 :
830 : static GEN
831 280 : _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
832 280 : { (void)data; return ZM_mul(x,y); }
833 : static GEN
834 24143 : _ZM_sqr(void *data /*ignored*/, GEN x)
835 24143 : { (void)data; return ZM_sqr(x); }
836 : /* FIXME: Using Bodrato's squaring, precomputations attached to fixed
837 : * multiplicand should be reused. And some postcomputations can be fused */
838 : GEN
839 0 : ZM_pow(GEN x, GEN n)
840 : {
841 0 : pari_sp av = avma;
842 0 : if (!signe(n)) return matid(lg(x)-1);
843 0 : return gc_GEN(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
844 : }
845 : GEN
846 23618 : ZM_powu(GEN x, ulong n)
847 : {
848 23618 : pari_sp av = avma;
849 23618 : if (!n) return matid(lg(x)-1);
850 23618 : return gc_GEN(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
851 : }
852 :
853 : GEN
854 0 : ZMV_prod(GEN v)
855 0 : { return gen_product(v, NULL, _ZM_mul); }
856 :
857 : /********************************************************************/
858 : /** **/
859 : /** ADD, SUB **/
860 : /** **/
861 : /********************************************************************/
862 : static GEN
863 37856433 : ZC_add_i(GEN x, GEN y, long lx)
864 : {
865 37856433 : GEN A = cgetg(lx, t_COL);
866 : long i;
867 534453867 : for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
868 37856433 : return A;
869 : }
870 : GEN
871 27825886 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
872 : GEN
873 370297 : ZC_Z_add(GEN x, GEN y)
874 : {
875 370297 : long k, lx = lg(x);
876 370297 : GEN z = cgetg(lx, t_COL);
877 370297 : if (lx == 1) pari_err_TYPE2("+",x,y);
878 370297 : gel(z,1) = addii(y,gel(x,1));
879 2533084 : for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
880 370297 : return z;
881 : }
882 :
883 : static GEN
884 9577973 : ZC_sub_i(GEN x, GEN y, long lx)
885 : {
886 : long i;
887 9577973 : GEN A = cgetg(lx, t_COL);
888 66421389 : for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
889 9577973 : return A;
890 : }
891 : GEN
892 9162751 : ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
893 : GEN
894 0 : ZC_Z_sub(GEN x, GEN y)
895 : {
896 0 : long k, lx = lg(x);
897 0 : GEN z = cgetg(lx, t_COL);
898 0 : if (lx == 1) pari_err_TYPE2("+",x,y);
899 0 : gel(z,1) = subii(gel(x,1), y);
900 0 : for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
901 0 : return z;
902 : }
903 : GEN
904 584242 : Z_ZC_sub(GEN a, GEN x)
905 : {
906 584242 : long k, lx = lg(x);
907 584242 : GEN z = cgetg(lx, t_COL);
908 584242 : if (lx == 1) pari_err_TYPE2("-",a,x);
909 584242 : gel(z,1) = subii(a, gel(x,1));
910 1596039 : for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
911 584242 : return z;
912 : }
913 :
914 : GEN
915 844758 : ZM_add(GEN x, GEN y)
916 : {
917 844758 : long lx = lg(x), l, j;
918 : GEN z;
919 844758 : if (lx == 1) return cgetg(1, t_MAT);
920 804676 : z = cgetg(lx, t_MAT); l = lgcols(x);
921 10835223 : for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
922 804676 : return z;
923 : }
924 : GEN
925 181433 : ZM_sub(GEN x, GEN y)
926 : {
927 181433 : long lx = lg(x), l, j;
928 : GEN z;
929 181433 : if (lx == 1) return cgetg(1, t_MAT);
930 141351 : z = cgetg(lx, t_MAT); l = lgcols(x);
931 556573 : for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
932 141351 : return z;
933 : }
934 : /********************************************************************/
935 : /** **/
936 : /** LINEAR COMBINATION **/
937 : /** **/
938 : /********************************************************************/
939 : /* return X/c assuming division is exact */
940 : GEN
941 5892984 : ZC_Z_divexact(GEN x, GEN c)
942 80670331 : { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
943 : GEN
944 2261 : ZC_divexactu(GEN x, ulong c)
945 11375 : { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
946 :
947 : GEN
948 530810 : ZM_Z_divexact(GEN x, GEN c)
949 4043963 : { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
950 :
951 : GEN
952 441 : ZM_divexactu(GEN x, ulong c)
953 2702 : { pari_APPLY_same(ZC_divexactu(gel(x,i), c)) }
954 :
955 : GEN
956 36989779 : ZC_Z_mul(GEN x, GEN c)
957 : {
958 36989779 : if (!signe(c)) return zerocol(lg(x)-1);
959 35455455 : if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
960 300412409 : pari_APPLY_type(t_COL, mulii(gel(x,i), c))
961 : }
962 :
963 : GEN
964 64224 : ZC_z_mul(GEN x, long c)
965 : {
966 64224 : if (!c) return zerocol(lg(x)-1);
967 54296 : if (c == 1) return ZC_copy(x);
968 49367 : if (c ==-1) return ZC_neg(x);
969 482526 : pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
970 : }
971 :
972 : GEN
973 28722 : zv_z_mul(GEN x, long n)
974 436904 : { pari_APPLY_long(x[i]*n) }
975 :
976 : /* return a ZM */
977 : GEN
978 0 : nm_Z_mul(GEN X, GEN c)
979 : {
980 0 : long i, j, h, l = lg(X), s = signe(c);
981 : GEN A;
982 0 : if (l == 1) return cgetg(1, t_MAT);
983 0 : h = lgcols(X);
984 0 : if (!s) return zeromat(h-1, l-1);
985 0 : if (is_pm1(c)) {
986 0 : if (s > 0) return Flm_to_ZM(X);
987 0 : X = Flm_to_ZM(X); ZM_togglesign(X); return X;
988 : }
989 0 : A = cgetg(l, t_MAT);
990 0 : for (j = 1; j < l; j++)
991 : {
992 0 : GEN a = cgetg(h, t_COL), x = gel(X, j);
993 0 : for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
994 0 : gel(A,j) = a;
995 : }
996 0 : return A;
997 : }
998 : GEN
999 3392510 : ZM_Z_mul(GEN X, GEN c)
1000 : {
1001 3392510 : long i, j, h, l = lg(X);
1002 : GEN A;
1003 3392510 : if (l == 1) return cgetg(1, t_MAT);
1004 3392510 : h = lgcols(X);
1005 3392510 : if (!signe(c)) return zeromat(h-1, l-1);
1006 3347139 : if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
1007 2930166 : A = cgetg(l, t_MAT);
1008 16004081 : for (j = 1; j < l; j++)
1009 : {
1010 13073915 : GEN a = cgetg(h, t_COL), x = gel(X, j);
1011 220357295 : for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
1012 13073915 : gel(A,j) = a;
1013 : }
1014 2930166 : return A;
1015 : }
1016 : void
1017 77355682 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
1018 : {
1019 : long i;
1020 1687880942 : for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
1021 77355682 : }
1022 : /* X <- X + v Y (elementary col operation) */
1023 : void
1024 66147664 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
1025 : {
1026 66147664 : if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
1027 : }
1028 : void
1029 31167276 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
1030 : {
1031 : long i;
1032 31167276 : if (!v) return; /* v = 0 */
1033 803153868 : for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
1034 : }
1035 :
1036 : /* X + v Y, wasteful if (v = 0) */
1037 : static GEN
1038 16046480 : ZC_lincomb1(GEN v, GEN x, GEN y)
1039 124241198 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
1040 :
1041 : /* -X + vY */
1042 : static GEN
1043 759404 : ZC_lincomb_1(GEN v, GEN x, GEN y)
1044 4657692 : { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
1045 :
1046 : /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
1047 : GEN
1048 34920453 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
1049 : {
1050 : long su, sv;
1051 : GEN A;
1052 :
1053 34920453 : su = signe(u); if (!su) return ZC_Z_mul(Y, v);
1054 34918185 : sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
1055 34201797 : if (is_pm1(v))
1056 : {
1057 11561260 : if (is_pm1(u))
1058 : {
1059 10215253 : if (su != sv) A = ZC_sub(X, Y);
1060 2832298 : else A = ZC_add(X, Y);
1061 10215253 : if (su < 0) ZV_togglesign(A); /* in place but was created above */
1062 : }
1063 : else
1064 : {
1065 1346007 : if (sv > 0) A = ZC_lincomb1 (u, Y, X);
1066 610588 : else A = ZC_lincomb_1(u, Y, X);
1067 : }
1068 : }
1069 22640537 : else if (is_pm1(u))
1070 : {
1071 15459877 : if (su > 0) A = ZC_lincomb1 (v, X, Y);
1072 148816 : else A = ZC_lincomb_1(v, X, Y);
1073 : }
1074 : else
1075 : { /* not cgetg_copy: x may be a t_VEC */
1076 7180660 : long i, lx = lg(X);
1077 7180660 : A = cgetg(lx,t_COL);
1078 45536212 : for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
1079 : }
1080 34201797 : return A;
1081 : }
1082 :
1083 : /********************************************************************/
1084 : /** **/
1085 : /** CONVERSIONS **/
1086 : /** **/
1087 : /********************************************************************/
1088 : GEN
1089 431621 : ZV_to_nv(GEN x)
1090 796094 : { pari_APPLY_ulong(itou(gel(x,i))) }
1091 :
1092 : GEN
1093 183511 : zm_to_ZM(GEN x)
1094 1030260 : { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
1095 :
1096 : GEN
1097 126 : zmV_to_ZMV(GEN x)
1098 791 : { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
1099 :
1100 : /* same as Flm_to_ZM but do not assume positivity */
1101 : GEN
1102 1022 : ZM_to_zm(GEN x)
1103 17472 : { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
1104 :
1105 : GEN
1106 366646 : zv_to_Flv(GEN x, ulong p)
1107 5418812 : { pari_APPLY_ulong(umodsu(x[i], p)) }
1108 :
1109 : GEN
1110 22694 : zm_to_Flm(GEN x, ulong p)
1111 351750 : { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
1112 :
1113 : GEN
1114 49 : ZMV_to_zmV(GEN x)
1115 399 : { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
1116 :
1117 : /********************************************************************/
1118 : /** **/
1119 : /** COPY, NEGATION **/
1120 : /** **/
1121 : /********************************************************************/
1122 : GEN
1123 18827188 : ZC_copy(GEN x)
1124 : {
1125 18827188 : long i, lx = lg(x);
1126 18827188 : GEN y = cgetg(lx, t_COL);
1127 141759283 : for (i=1; i<lx; i++)
1128 : {
1129 122932095 : GEN c = gel(x,i);
1130 122932095 : gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
1131 : }
1132 18827188 : return y;
1133 : }
1134 :
1135 : GEN
1136 750004 : ZM_copy(GEN x)
1137 6754780 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
1138 :
1139 : void
1140 391786 : ZV_neg_inplace(GEN M)
1141 : {
1142 391786 : long l = lg(M);
1143 1614810 : while (--l > 0) gel(M,l) = negi(gel(M,l));
1144 391786 : }
1145 : GEN
1146 6378023 : ZC_neg(GEN x)
1147 37040838 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
1148 :
1149 : GEN
1150 51880 : zv_neg(GEN x)
1151 662825 : { pari_APPLY_long(-x[i]) }
1152 : GEN
1153 126 : zv_neg_inplace(GEN M)
1154 : {
1155 126 : long l = lg(M);
1156 427 : while (--l > 0) M[l] = -M[l];
1157 126 : return M;
1158 : }
1159 : GEN
1160 77 : zv_abs(GEN x)
1161 5446 : { pari_APPLY_ulong(labs(x[i])) }
1162 : GEN
1163 1734732 : ZM_neg(GEN x)
1164 5197619 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
1165 : int
1166 22921521 : zv_canon_inplace(GEN x)
1167 : {
1168 22921521 : long l = lg(x), j, k;
1169 41338311 : for (j = 1; j < l && x[j] == 0; ++j);
1170 22921521 : if (j < l && x[j] < 0)
1171 : {
1172 52277022 : for (k = j; k < l; ++k) x[k] = -x[k];
1173 10255819 : return -1;
1174 : }
1175 12665702 : return 1;
1176 : }
1177 :
1178 : void
1179 7284631 : ZV_togglesign(GEN M)
1180 : {
1181 7284631 : long l = lg(M);
1182 97146799 : while (--l > 0) togglesign_safe(&gel(M,l));
1183 7284631 : }
1184 : void
1185 0 : ZM_togglesign(GEN M)
1186 : {
1187 0 : long l = lg(M);
1188 0 : while (--l > 0) ZV_togglesign(gel(M,l));
1189 0 : }
1190 :
1191 : /********************************************************************/
1192 : /** **/
1193 : /** "DIVISION" mod HNF **/
1194 : /** **/
1195 : /********************************************************************/
1196 : /* Reduce ZC x modulo ZM y in HNF */
1197 : static GEN
1198 11406113 : ZC_hnfdivrem_i(GEN x, GEN y, GEN *Q, GEN (*div)(GEN,GEN))
1199 : {
1200 11406113 : long i, l = lg(x);
1201 11406113 : pari_sp av = avma;
1202 :
1203 11406113 : if (Q) *Q = cgetg(l,t_COL);
1204 11406113 : if (l == 1) return cgetg(1,t_COL);
1205 64214680 : for (i = l-1; i>0; i--)
1206 : {
1207 52808567 : GEN q = div(gel(x,i), gcoeff(y,i,i));
1208 52808567 : if (signe(q)) x = ZC_lincomb(gen_1, negi(q), x, gel(y,i));
1209 52808567 : if (Q) gel(*Q, i) = q;
1210 : }
1211 11406113 : if (avma == av) return ZC_copy(x);
1212 8558877 : if (!Q) return gc_upto(av, x);
1213 60894 : (void)gc_all(av, 2, &x, Q); return x;
1214 : }
1215 : GEN
1216 10699892 : ZC_hnfdivrem(GEN x, GEN y, GEN *Q)
1217 10699892 : { return ZC_hnfdivrem_i(x, y, Q, diviiround); }
1218 : GEN
1219 280 : ZC_modhnf(GEN x, GEN y, GEN *Q)
1220 280 : { return ZC_hnfdivrem_i(x, y, Q, truedivii); }
1221 :
1222 : /* Return R such that x = y Q + R, y integral HNF */
1223 : static GEN
1224 529693 : ZM_hnfdivrem_i(GEN x, GEN y, GEN *Q, GEN (*div)(GEN,GEN))
1225 : {
1226 : long l, i;
1227 529693 : GEN R = cgetg_copy(x, &l);
1228 529693 : if (Q)
1229 : {
1230 128135 : GEN q = cgetg(l, t_MAT); *Q = q;
1231 189015 : for (i = 1; i < l; i++)
1232 60880 : gel(R,i) = ZC_hnfdivrem_i(gel(x,i),y,&gel(q,i),div);
1233 : }
1234 : else
1235 1046619 : for (i = 1; i < l; i++)
1236 645061 : gel(R,i) = ZC_hnfdivrem_i(gel(x,i),y,NULL,div);
1237 529693 : return R;
1238 : }
1239 : GEN
1240 529679 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
1241 529679 : { return ZM_hnfdivrem_i(x, y, Q, diviiround); }
1242 : GEN
1243 14 : ZM_modhnf(GEN x, GEN y, GEN *Q)
1244 14 : { return ZM_hnfdivrem_i(x, y, Q, truedivii); }
1245 :
1246 : static GEN
1247 7 : ZV_ZV_divrem(GEN x, GEN y, GEN *pQ)
1248 : {
1249 7 : long i, l = lg(x), tx = typ(x);
1250 : GEN Q, R;
1251 :
1252 7 : if (!pQ) return ZV_ZV_mod(x, y);
1253 0 : Q = cgetg(l,tx);
1254 0 : R = cgetg(l,tx);
1255 0 : for (i = 1; i < l; i++) gel(Q,i) = truedvmdii(gel(x,i), gel(y,i), &gel(R,i));
1256 0 : *pQ = Q; return R;
1257 : }
1258 : static GEN
1259 0 : ZM_ZV_divrem(GEN x, GEN y, GEN *Q)
1260 0 : { if (!Q) return ZM_ZV_mod(x, y);
1261 0 : pari_APPLY_same(ZV_ZV_divrem(gel(x,i), y, Q)); }
1262 :
1263 : static int
1264 42 : RgM_issquare(GEN x) { long l = lg(x); return l == 1 || lg(gel(x,1)) == l; }
1265 : static void
1266 98 : matmodhnf_check(GEN x)
1267 : {
1268 98 : switch(typ(x))
1269 : {
1270 42 : case t_VEC: case t_COL:
1271 42 : if (!RgV_is_ZV(x)) pari_err_TYPE("matmodhnf", x);
1272 42 : break;
1273 56 : case t_MAT:
1274 56 : if (!RgM_is_ZM(x)) pari_err_TYPE("matmodhnf", x);
1275 56 : break;
1276 0 : default: pari_err_TYPE("matmodhnf", x);
1277 : }
1278 98 : }
1279 : GEN
1280 49 : matmodhnf(GEN x, GEN y, GEN *Q)
1281 : {
1282 49 : long tx = typ(x), ty = typ(y), ly, lx;
1283 49 : matmodhnf_check(x); lx = lg(x);
1284 49 : matmodhnf_check(y); ly = lg(y);
1285 49 : if (ty == t_MAT && !RgM_issquare(y)) pari_err_TYPE("matmodhnf", y);
1286 49 : if (tx == t_MAT && lx == 1)
1287 : {
1288 0 : if (ly != 1) pari_err_DIM("matmodhnf");
1289 0 : if (!Q) *Q = cgetg(1, t_MAT);
1290 0 : return cgetg(1, t_MAT);
1291 : }
1292 49 : if (is_vec_t(ty))
1293 7 : return tx == t_MAT? ZM_ZV_divrem(x, y, Q): ZV_ZV_divrem(x, y, Q);
1294 : /* ty = t_MAT */
1295 42 : if (tx == t_MAT) return ZM_modhnf(x, y, Q);
1296 28 : x = ZC_modhnf(x, y, Q);
1297 28 : if (tx == t_VEC) { settyp(x, tx); if (Q) settyp(*Q, tx); }
1298 28 : return x;
1299 : }
1300 :
1301 : /********************************************************************/
1302 : /** **/
1303 : /** TESTS **/
1304 : /** **/
1305 : /********************************************************************/
1306 : int
1307 23113620 : zv_equal0(GEN V)
1308 : {
1309 23113620 : long l = lg(V);
1310 37599001 : while (--l > 0)
1311 30822295 : if (V[l]) return 0;
1312 6776706 : return 1;
1313 : }
1314 :
1315 : int
1316 14226697 : ZV_equal0(GEN V)
1317 : {
1318 14226697 : long l = lg(V);
1319 25441262 : while (--l > 0)
1320 24991465 : if (signe(gel(V,l))) return 0;
1321 449797 : return 1;
1322 : }
1323 : int
1324 16231 : ZMrow_equal0(GEN V, long i)
1325 : {
1326 16231 : long l = lg(V);
1327 25183 : while (--l > 0)
1328 21679 : if (signe(gcoeff(V,i,l))) return 0;
1329 3504 : return 1;
1330 : }
1331 :
1332 : static int
1333 6528520 : ZV_equal_lg(GEN V, GEN W, long l)
1334 : {
1335 26972809 : while (--l > 0)
1336 20949265 : if (!equalii(gel(V,l), gel(W,l))) return 0;
1337 6023544 : return 1;
1338 : }
1339 : int
1340 293378 : ZV_equal(GEN V, GEN W)
1341 : {
1342 293378 : long l = lg(V);
1343 293378 : if (lg(W) != l) return 0;
1344 293357 : return ZV_equal_lg(V, W, l);
1345 : }
1346 : int
1347 3503816 : ZM_equal(GEN A, GEN B)
1348 : {
1349 3503816 : long i, m, l = lg(A);
1350 3503816 : if (lg(B) != l) return 0;
1351 3503762 : if (l == 1) return 1;
1352 3503762 : m = lgcols(A);
1353 3503762 : if (lgcols(B) != m) return 0;
1354 9452019 : for (i = 1; i < l; i++)
1355 6235163 : if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
1356 3216856 : return 1;
1357 : }
1358 : int
1359 32306 : ZM_equal0(GEN A)
1360 : {
1361 32306 : long i, j, m, l = lg(A);
1362 32306 : if (l == 1) return 1;
1363 32306 : m = lgcols(A);
1364 88094 : for (j = 1; j < l; j++)
1365 2730075 : for (i = 1; i < m; i++)
1366 2674287 : if (signe(gcoeff(A,i,j))) return 0;
1367 17195 : return 1;
1368 : }
1369 : int
1370 68047584 : zv_equal(GEN V, GEN W)
1371 : {
1372 68047584 : long l = lg(V);
1373 68047584 : if (lg(W) != l) return 0;
1374 490781610 : while (--l > 0)
1375 424036269 : if (V[l] != W[l]) return 0;
1376 66745341 : return 1;
1377 : }
1378 :
1379 : int
1380 1638 : zvV_equal(GEN V, GEN W)
1381 : {
1382 1638 : long l = lg(V);
1383 1638 : if (lg(W) != l) return 0;
1384 80388 : while (--l > 0)
1385 79912 : if (!zv_equal(gel(V,l),gel(W,l))) return 0;
1386 476 : return 1;
1387 : }
1388 :
1389 : int
1390 806949 : ZM_ishnf(GEN x)
1391 : {
1392 806949 : long i,j, lx = lg(x);
1393 2529713 : for (i=1; i<lx; i++)
1394 : {
1395 1792361 : GEN xii = gcoeff(x,i,i);
1396 1792361 : if (signe(xii) <= 0) return 0;
1397 3500497 : for (j=1; j<i; j++)
1398 1715528 : if (signe(gcoeff(x,i,j))) return 0;
1399 3524087 : for (j=i+1; j<lx; j++)
1400 : {
1401 1801323 : GEN xij = gcoeff(x,i,j);
1402 1801323 : if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
1403 : }
1404 : }
1405 737352 : return 1;
1406 : }
1407 : int
1408 48524 : QM_ishnf(GEN x)
1409 : {
1410 48524 : long i,j, lx = lg(x);
1411 94795 : for (i=1; i<lx; i++)
1412 : {
1413 89538 : GEN xii = gcoeff(x,i,i);
1414 89538 : if (gsigne(xii) <= 0) return 0;
1415 80283 : for (j=1; j<i; j++)
1416 20867 : if (gsigne(gcoeff(x,i,j))) return 0;
1417 115298 : for (j=i+1; j<lx; j++)
1418 : {
1419 69027 : GEN xij = gcoeff(x,i,j);
1420 69027 : if (gsigne(xij)<0 || gcmp(xij,xii)>=0) return 0;
1421 : }
1422 : }
1423 5257 : return 1;
1424 : }
1425 : int
1426 658577 : ZM_isidentity(GEN x)
1427 : {
1428 658577 : long i,j, lx = lg(x);
1429 :
1430 658577 : if (lx == 1) return 1;
1431 658570 : if (lx != lgcols(x)) return 0;
1432 3223320 : for (j=1; j<lx; j++)
1433 : {
1434 2564828 : GEN c = gel(x,j);
1435 8110569 : for (i=1; i<j; )
1436 5545763 : if (signe(gel(c,i++))) return 0;
1437 : /* i = j */
1438 2564806 : if (!equali1(gel(c,i++))) return 0;
1439 8110584 : for ( ; i<lx; )
1440 5545799 : if (signe(gel(c,i++))) return 0;
1441 : }
1442 658492 : return 1;
1443 : }
1444 : int
1445 677498 : ZM_isdiagonal(GEN x)
1446 : {
1447 677498 : long i,j, lx = lg(x);
1448 677498 : if (lx == 1) return 1;
1449 677498 : if (lx != lgcols(x)) return 0;
1450 :
1451 1738179 : for (j=1; j<lx; j++)
1452 : {
1453 1462260 : GEN c = gel(x,j);
1454 2029210 : for (i=1; i<j; i++)
1455 968494 : if (signe(gel(c,i))) return 0;
1456 2451478 : for (i++; i<lx; i++)
1457 1390797 : if (signe(gel(c,i))) return 0;
1458 : }
1459 275919 : return 1;
1460 : }
1461 : int
1462 162818 : ZM_isscalar(GEN x, GEN s)
1463 : {
1464 162818 : long i, j, lx = lg(x);
1465 :
1466 162818 : if (lx == 1) return 1;
1467 162818 : if (!s) s = gcoeff(x,1,1);
1468 162818 : if (equali1(s)) return ZM_isidentity(x);
1469 161362 : if (lx != lgcols(x)) return 0;
1470 726209 : for (j=1; j<lx; j++)
1471 : {
1472 627250 : GEN c = gel(x,j);
1473 2911276 : for (i=1; i<j; )
1474 2344331 : if (signe(gel(c,i++))) return 0;
1475 : /* i = j */
1476 566945 : if (!equalii(gel(c,i++), s)) return 0;
1477 2939644 : for ( ; i<lx; )
1478 2374797 : if (signe(gel(c,i++))) return 0;
1479 : }
1480 98959 : return 1;
1481 : }
1482 :
1483 : long
1484 174510 : ZC_is_ei(GEN x)
1485 : {
1486 174510 : long i, j = 0, l = lg(x);
1487 1790445 : for (i = 1; i < l; i++)
1488 : {
1489 1615936 : GEN c = gel(x,i);
1490 1615936 : long s = signe(c);
1491 1615936 : if (!s) continue;
1492 174504 : if (s < 0 || !is_pm1(c) || j) return 0;
1493 174503 : j = i;
1494 : }
1495 174509 : return j;
1496 : }
1497 :
1498 : /********************************************************************/
1499 : /** **/
1500 : /** MISCELLANEOUS **/
1501 : /** **/
1502 : /********************************************************************/
1503 : /* assume lg(x) = lg(y), x,y in Z^n */
1504 : int
1505 3256131 : ZV_cmp(GEN x, GEN y)
1506 : {
1507 3256131 : long fl,i, lx = lg(x);
1508 6522733 : for (i=1; i<lx; i++)
1509 5209242 : if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
1510 1313491 : return 0;
1511 : }
1512 : /* assume lg(x) = lg(y), x,y in Z^n */
1513 : int
1514 19747 : ZV_abscmp(GEN x, GEN y)
1515 : {
1516 19747 : long fl,i, lx = lg(x);
1517 54113 : for (i=1; i<lx; i++)
1518 53984 : if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
1519 129 : return 0;
1520 : }
1521 :
1522 : long
1523 19541151 : zv_content(GEN x)
1524 : {
1525 19541151 : long i, s, l = lg(x);
1526 19541151 : if (l == 1) return 0;
1527 19541144 : s = labs(x[1]);
1528 44329780 : for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
1529 19541144 : return s;
1530 : }
1531 : GEN
1532 300484 : ZV_content(GEN x)
1533 : {
1534 300484 : long i, l = lg(x);
1535 300484 : pari_sp av = avma;
1536 : GEN c;
1537 300484 : if (l == 1) return gen_0;
1538 300484 : if (l == 2) return absi(gel(x,1));
1539 207517 : c = gel(x,1);
1540 562637 : for (i = 2; i < l; i++)
1541 : {
1542 407313 : c = gcdii(c, gel(x,i));
1543 407313 : if (is_pm1(c)) { set_avma(av); return gen_1; }
1544 : }
1545 155324 : return gc_INT(av, c);
1546 : }
1547 :
1548 : GEN
1549 4270765 : ZM_det_triangular(GEN mat)
1550 : {
1551 : pari_sp av;
1552 4270765 : long i,l = lg(mat);
1553 : GEN s;
1554 :
1555 4270765 : if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
1556 3838717 : av = avma; s = gcoeff(mat,1,1);
1557 10370296 : for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
1558 3838717 : return gc_INT(av,s);
1559 : }
1560 :
1561 : /* assumes no overflow */
1562 : long
1563 958465 : zv_prod(GEN v)
1564 : {
1565 958465 : long n, i, l = lg(v);
1566 958465 : if (l == 1) return 1;
1567 969354 : n = v[1]; for (i = 2; i < l; i++) n *= v[i];
1568 779916 : return n;
1569 : }
1570 :
1571 : static GEN
1572 319253540 : _mulii(void *E, GEN a, GEN b)
1573 319253540 : { (void) E; return mulii(a, b); }
1574 :
1575 : /* product of ulongs */
1576 : GEN
1577 1976461 : zv_prod_Z(GEN v)
1578 : {
1579 : pari_sp av;
1580 1976461 : long k, m, n = lg(v)-1;
1581 1976461 : int stop = 0;
1582 : GEN V;
1583 1976461 : switch(n) {
1584 21056 : case 0: return gen_1;
1585 129234 : case 1: return utoi(v[1]);
1586 1038124 : case 2: return muluu(v[1], v[2]);
1587 : }
1588 788047 : av = avma; m = n >> 1;
1589 788047 : V = cgetg(m + (odd(n)? 2: 1), t_VEC);
1590 155539366 : for (k = n; k; k--) /* start from the end: v is usually sorted */
1591 154752401 : if (v[k] & HIGHMASK) { stop = 1; break; }
1592 2631134 : while (!stop)
1593 : { /* HACK: handle V as a t_VECSMALL; gain a few iterations */
1594 84039262 : for (k = 1; k <= m; k++)
1595 : {
1596 81585376 : V[k] = uel(v,k<<1) * uel(v,(k<<1)-1);
1597 81585376 : if (V[k] & HIGHMASK) stop = 1; /* last "free" iteration */
1598 : }
1599 2453886 : if (odd(n))
1600 : {
1601 1432336 : if (n == 1) { set_avma(av); return utoi(v[1]); }
1602 821537 : V[++m] = v[n];
1603 : }
1604 1843087 : v = V; n = m; m = n >> 1;
1605 : }
1606 : /* n > 1; m > 0 */
1607 177248 : if (n == 2) { set_avma(av); return muluu(v[1], v[2]); }
1608 47497390 : for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
1609 123601 : if (odd(n)) gel(V, ++m) = utoipos(v[n]);
1610 123601 : setlg(V, m+1); /* HACK: now V is a bona fide t_VEC */
1611 123601 : return gc_INT(av, gen_product(V, NULL, &_mulii));
1612 : }
1613 : GEN
1614 14694393 : vecsmall_prod(GEN v)
1615 : {
1616 14694393 : pari_sp av = avma;
1617 14694393 : long k, m, n = lg(v)-1;
1618 : GEN V;
1619 14694393 : switch (n) {
1620 0 : case 0: return gen_1;
1621 0 : case 1: return stoi(v[1]);
1622 21 : case 2: return mulss(v[1], v[2]);
1623 : }
1624 14694372 : m = n >> 1;
1625 14694372 : V = cgetg(m + (odd(n)? 2: 1), t_VEC);
1626 161556906 : for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
1627 14694372 : if (odd(n)) gel(V,k) = stoi(v[n]);
1628 14694372 : return gc_INT(av, gen_product(V, NULL, &_mulii));
1629 : }
1630 :
1631 : GEN
1632 8856317 : ZV_prod(GEN v)
1633 : {
1634 8856317 : pari_sp av = avma;
1635 8856317 : long i, l = lg(v);
1636 : GEN n;
1637 8856317 : if (l == 1) return gen_1;
1638 8671013 : if (l > 7) return gc_INT(av, gen_product(v, NULL, _mulii));
1639 1325460 : n = gel(v,1);
1640 1325460 : if (l == 2) return icopy(n);
1641 2184642 : for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
1642 867558 : return gc_INT(av, n);
1643 : }
1644 : /* assumes no overflow */
1645 : long
1646 16525 : zv_sum(GEN v)
1647 : {
1648 16525 : long n, i, l = lg(v);
1649 16525 : if (l == 1) return 0;
1650 99703 : n = v[1]; for (i = 2; i < l; i++) n += v[i];
1651 16504 : return n;
1652 : }
1653 : /* assumes no overflow and 0 <= n <= #v */
1654 : long
1655 0 : zv_sumpart(GEN v, long n)
1656 : {
1657 : long i, p;
1658 0 : if (!n) return 0;
1659 0 : p = v[1]; for (i = 2; i <= n; i++) p += v[i];
1660 0 : return p;
1661 : }
1662 : GEN
1663 77 : ZV_sum(GEN v)
1664 : {
1665 77 : pari_sp av = avma;
1666 77 : long i, l = lg(v);
1667 : GEN n;
1668 77 : if (l == 1) return gen_0;
1669 77 : n = gel(v,1);
1670 77 : if (l == 2) return icopy(n);
1671 581 : for (i = 2; i < l; i++) n = addii(n, gel(v,i));
1672 77 : return gc_INT(av, n);
1673 : }
1674 :
1675 : /********************************************************************/
1676 : /** **/
1677 : /** GRAM SCHMIDT REDUCTION (integer matrices) **/
1678 : /** **/
1679 : /********************************************************************/
1680 :
1681 : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
1682 : static void
1683 402314 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
1684 : {
1685 402314 : long i, qq = itos_or_0(q);
1686 402314 : if (!qq)
1687 : {
1688 36730 : for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
1689 7587 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
1690 7587 : return;
1691 : }
1692 394727 : if (qq == 1) {
1693 182194 : for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
1694 123623 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
1695 271104 : } else if (qq == -1) {
1696 174777 : for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
1697 106913 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
1698 : } else {
1699 290973 : for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
1700 164191 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
1701 : }
1702 : }
1703 :
1704 : /* update L[k,] */
1705 : static void
1706 1170949 : ZRED(long k, long l, GEN x, GEN L, GEN B)
1707 : {
1708 1170949 : GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
1709 1170949 : if (!signe(q)) return;
1710 402314 : q = negi(q);
1711 402314 : Zupdate_row(k,l,q,L,B);
1712 402314 : gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
1713 : }
1714 :
1715 : /* Gram-Schmidt reduction, x a ZM */
1716 : static void
1717 1279495 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
1718 : {
1719 : long i, j;
1720 4130456 : for (j=1; j<=k; j++)
1721 : {
1722 2850961 : pari_sp av = avma;
1723 2850961 : GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
1724 6125547 : for (i=1; i<j; i++)
1725 : {
1726 3274586 : u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
1727 3274586 : u = diviiexact(u, gel(B,i));
1728 : }
1729 2850961 : gcoeff(L,k,j) = gc_INT(av, u);
1730 : }
1731 1279495 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
1732 1279495 : }
1733 :
1734 : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
1735 : * Very inefficient if y is not LLL-reduced of maximal rank */
1736 : static GEN
1737 112326 : ZC_reducemodmatrix_i(GEN v, GEN y)
1738 : {
1739 112326 : GEN B, L, x = shallowconcat(y, v);
1740 112326 : long k, lx = lg(x), nx = lx-1;
1741 :
1742 112326 : B = scalarcol_shallow(gen_1, lx);
1743 112326 : L = zeromatcopy(nx, nx);
1744 462032 : for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
1745 349706 : for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1746 112326 : return gel(x,nx);
1747 : }
1748 : GEN
1749 112326 : ZC_reducemodmatrix(GEN v, GEN y) {
1750 112326 : pari_sp av = avma;
1751 112326 : return gc_GEN(av, ZC_reducemodmatrix_i(v,y));
1752 : }
1753 :
1754 : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
1755 : * Very inefficient if y is not LLL-reduced of maximal rank */
1756 : static GEN
1757 239883 : ZM_reducemodmatrix_i(GEN v, GEN y)
1758 : {
1759 : GEN B, L, V;
1760 239883 : long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
1761 :
1762 239883 : V = cgetg(lv, t_MAT);
1763 239883 : B = scalarcol_shallow(gen_1, lx);
1764 239883 : L = zeromatcopy(nx, nx);
1765 640647 : for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
1766 768908 : for (j = 1; j < lg(v); j++)
1767 : {
1768 529025 : GEN x = shallowconcat(y, gel(v,j));
1769 529025 : ZincrementalGS(x, L, B, nx); /* overwrite last */
1770 1462594 : for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1771 529025 : gel(V,j) = gel(x,nx);
1772 : }
1773 239883 : return V;
1774 : }
1775 : GEN
1776 239883 : ZM_reducemodmatrix(GEN v, GEN y) {
1777 239883 : pari_sp av = avma;
1778 239883 : return gc_GEN(av, ZM_reducemodmatrix_i(v,y));
1779 : }
1780 :
1781 : GEN
1782 98016 : ZC_reducemodlll(GEN x,GEN y)
1783 : {
1784 98016 : pari_sp av = avma;
1785 98016 : GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
1786 98016 : return gc_GEN(av, z);
1787 : }
1788 : GEN
1789 0 : ZM_reducemodlll(GEN x,GEN y)
1790 : {
1791 0 : pari_sp av = avma;
1792 0 : GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
1793 0 : return gc_GEN(av, z);
1794 : }
|
