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
check_ZV(GEN x,long l)19 check_ZV(GEN x, long l)
20 {
21 long i;
22 for (i=1; i<l; i++)
23 if (typ(gel(x,i)) != t_INT) return 0;
24 return 1;
25 }
26 void
RgV_check_ZV(GEN A,const char * s)27 RgV_check_ZV(GEN A, const char *s)
28 {
29 if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
30 }
31 void
RgM_check_ZM(GEN A,const char * s)32 RgM_check_ZM(GEN A, const char *s)
33 {
34 long n = lg(A);
35 if (n != 1)
36 {
37 long j, m = lgcols(A);
38 for (j=1; j<n; j++)
39 if (!check_ZV(gel(A,j), m))
40 pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
41 }
42 }
43
44 static long
ZV_max_lg_i(GEN x,long m)45 ZV_max_lg_i(GEN x, long m)
46 {
47 long i, prec = 2;
48 for (i=1; i<m; i++) { long l = lgefint(gel(x,i)); if (l > prec) prec = l; }
49 return prec;
50 }
51
52 long
ZV_max_lg(GEN x)53 ZV_max_lg(GEN x)
54 { return ZV_max_lg_i(x, lg(x)); }
55
56 static long
ZM_max_lg_i(GEN x,long n,long m)57 ZM_max_lg_i(GEN x, long n, long m)
58 {
59 long prec = 2;
60 if (n != 1)
61 {
62 long j;
63 for (j=1; j<n; j++)
64 {
65 long l = ZV_max_lg_i(gel(x,j), m);
66 if (l > prec) prec = l;
67 }
68 }
69 return prec;
70 }
71
72 long
ZM_max_lg(GEN x)73 ZM_max_lg(GEN x)
74 {
75 long n = lg(x);
76 if (n==1) return 2;
77 return ZM_max_lg_i(x, n, lgcols(x));
78 }
79
80 GEN
ZM_supnorm(GEN x)81 ZM_supnorm(GEN x)
82 {
83 long i, j, h, lx = lg(x);
84 GEN s = gen_0;
85 if (lx == 1) return gen_1;
86 h = lgcols(x);
87 for (j=1; j<lx; j++)
88 {
89 GEN xj = gel(x,j);
90 for (i=1; i<h; i++)
91 {
92 GEN c = gel(xj,i);
93 if (abscmpii(c, s) > 0) s = c;
94 }
95 }
96 return absi(s);
97 }
98
99 /********************************************************************/
100 /** **/
101 /** MULTIPLICATION **/
102 /** **/
103 /********************************************************************/
104 /* x nonempty ZM, y a compatible nc (dimension > 0). */
105 static GEN
ZM_nc_mul_i(GEN x,GEN y,long c,long l)106 ZM_nc_mul_i(GEN x, GEN y, long c, long l)
107 {
108 long i, j;
109 pari_sp av;
110 GEN z = cgetg(l,t_COL), s;
111
112 for (i=1; i<l; i++)
113 {
114 av = avma; s = muliu(gcoeff(x,i,1),y[1]);
115 for (j=2; j<c; j++)
116 if (y[j]) s = addii(s, muliu(gcoeff(x,i,j),y[j]));
117 gel(z,i) = gerepileuptoint(av,s);
118 }
119 return z;
120 }
121
122 /* x ZV, y a compatible zc. */
123 GEN
ZV_zc_mul(GEN x,GEN y)124 ZV_zc_mul(GEN x, GEN y)
125 {
126 long j, l = lg(x);
127 pari_sp av = avma;
128 GEN s = mulis(gel(x,1),y[1]);
129 for (j=2; j<l; j++)
130 if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
131 return gerepileuptoint(av,s);
132 }
133
134 /* x nonempty ZM, y a compatible zc (dimension > 0). */
135 static GEN
ZM_zc_mul_i(GEN x,GEN y,long c,long l)136 ZM_zc_mul_i(GEN x, GEN y, long c, long l)
137 {
138 long i, j;
139 GEN z = cgetg(l,t_COL);
140
141 for (i=1; i<l; i++)
142 {
143 pari_sp av = avma;
144 GEN s = mulis(gcoeff(x,i,1),y[1]);
145 for (j=2; j<c; j++)
146 if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
147 gel(z,i) = gerepileuptoint(av,s);
148 }
149 return z;
150 }
151 GEN
ZM_zc_mul(GEN x,GEN y)152 ZM_zc_mul(GEN x, GEN y) {
153 long l = lg(x);
154 if (l == 1) return cgetg(1, t_COL);
155 return ZM_zc_mul_i(x,y, l, lgcols(x));
156 }
157
158 /* y nonempty ZM, x a compatible zv (dimension > 0). */
159 GEN
zv_ZM_mul(GEN x,GEN y)160 zv_ZM_mul(GEN x, GEN y) {
161 long i,j, lx = lg(x), ly = lg(y);
162 GEN z;
163 if (lx == 1) return zerovec(ly-1);
164 z = cgetg(ly,t_VEC);
165 for (j=1; j<ly; j++)
166 {
167 pari_sp av = avma;
168 GEN s = mulsi(x[1], gcoeff(y,1,j));
169 for (i=2; i<lx; i++)
170 if (x[i]) s = addii(s, mulsi(x[i], gcoeff(y,i,j)));
171 gel(z,j) = gerepileuptoint(av,s);
172 }
173 return z;
174 }
175
176 /* x ZM, y a compatible zm (dimension > 0). */
177 GEN
ZM_zm_mul(GEN x,GEN y)178 ZM_zm_mul(GEN x, GEN y)
179 {
180 long j, c, l = lg(x), ly = lg(y);
181 GEN z = cgetg(ly, t_MAT);
182 if (l == 1) return z;
183 c = lgcols(x);
184 for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
185 return z;
186 }
187 /* x ZM, y a compatible zn (dimension > 0). */
188 GEN
ZM_nm_mul(GEN x,GEN y)189 ZM_nm_mul(GEN x, GEN y)
190 {
191 long j, c, l = lg(x), ly = lg(y);
192 GEN z = cgetg(ly, t_MAT);
193 if (l == 1) return z;
194 c = lgcols(x);
195 for (j = 1; j < ly; j++) gel(z,j) = ZM_nc_mul_i(x, gel(y,j), l,c);
196 return z;
197 }
198
199 /* Strassen-Winograd algorithm */
200
201 /*
202 Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
203 as an (m x n)-matrix, padding the input with zeroes as necessary.
204 */
205 static GEN
add_slices(long m,long n,GEN A,long ma,long da,long na,long ea,GEN B,long mb,long db,long nb,long eb)206 add_slices(long m, long n,
207 GEN A, long ma, long da, long na, long ea,
208 GEN B, long mb, long db, long nb, long eb)
209 {
210 long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
211 GEN M = cgetg(n + 1, t_MAT), C;
212
213 for (j = 1; j <= min_e; j++) {
214 gel(M, j) = C = cgetg(m + 1, t_COL);
215 for (i = 1; i <= min_d; i++)
216 gel(C, i) = addii(gcoeff(A, ma + i, na + j),
217 gcoeff(B, mb + i, nb + j));
218 for (; i <= da; i++)
219 gel(C, i) = gcoeff(A, ma + i, na + j);
220 for (; i <= db; i++)
221 gel(C, i) = gcoeff(B, mb + i, nb + j);
222 for (; i <= m; i++)
223 gel(C, i) = gen_0;
224 }
225 for (; j <= ea; j++) {
226 gel(M, j) = C = cgetg(m + 1, t_COL);
227 for (i = 1; i <= da; i++)
228 gel(C, i) = gcoeff(A, ma + i, na + j);
229 for (; i <= m; i++)
230 gel(C, i) = gen_0;
231 }
232 for (; j <= eb; j++) {
233 gel(M, j) = C = cgetg(m + 1, t_COL);
234 for (i = 1; i <= db; i++)
235 gel(C, i) = gcoeff(B, mb + i, nb + j);
236 for (; i <= m; i++)
237 gel(C, i) = gen_0;
238 }
239 for (; j <= n; j++)
240 gel(M, j) = zerocol(m);
241 return M;
242 }
243
244 /*
245 Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
246 as an (m x n)-matrix, padding the input with zeroes as necessary.
247 */
248 static GEN
subtract_slices(long m,long n,GEN A,long ma,long da,long na,long ea,GEN B,long mb,long db,long nb,long eb)249 subtract_slices(long m, long n,
250 GEN A, long ma, long da, long na, long ea,
251 GEN B, long mb, long db, long nb, long eb)
252 {
253 long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
254 GEN M = cgetg(n + 1, t_MAT), C;
255
256 for (j = 1; j <= min_e; j++) {
257 gel(M, j) = C = cgetg(m + 1, t_COL);
258 for (i = 1; i <= min_d; i++)
259 gel(C, i) = subii(gcoeff(A, ma + i, na + j),
260 gcoeff(B, mb + i, nb + j));
261 for (; i <= da; i++)
262 gel(C, i) = gcoeff(A, ma + i, na + j);
263 for (; i <= db; i++)
264 gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
265 for (; i <= m; i++)
266 gel(C, i) = gen_0;
267 }
268 for (; j <= ea; j++) {
269 gel(M, j) = C = cgetg(m + 1, t_COL);
270 for (i = 1; i <= da; i++)
271 gel(C, i) = gcoeff(A, ma + i, na + j);
272 for (; i <= m; i++)
273 gel(C, i) = gen_0;
274 }
275 for (; j <= eb; j++) {
276 gel(M, j) = C = cgetg(m + 1, t_COL);
277 for (i = 1; i <= db; i++)
278 gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
279 for (; i <= m; i++)
280 gel(C, i) = gen_0;
281 }
282 for (; j <= n; j++)
283 gel(M, j) = zerocol(m);
284 return M;
285 }
286
287 static GEN ZM_mul_i(GEN x, GEN y, long l, long lx, long ly);
288
289 /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
290 static GEN
ZM_mul_sw(GEN A,GEN B,long m,long n,long p)291 ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
292 {
293 pari_sp av = avma;
294 long m1 = (m + 1)/2, m2 = m/2,
295 n1 = (n + 1)/2, n2 = n/2,
296 p1 = (p + 1)/2, p2 = p/2;
297 GEN A11, A12, A22, B11, B21, B22,
298 S1, S2, S3, S4, T1, T2, T3, T4,
299 M1, M2, M3, M4, M5, M6, M7,
300 V1, V2, V3, C11, C12, C21, C22, C;
301
302 T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
303 S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
304 M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
305 if (gc_needed(av, 1))
306 gerepileall(av, 2, &T2, &M2); /* destroy S1 */
307 T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
308 if (gc_needed(av, 1))
309 gerepileall(av, 2, &M2, &T3); /* destroy T2 */
310 S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
311 T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
312 M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
313 if (gc_needed(av, 1))
314 gerepileall(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
315 S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
316 if (gc_needed(av, 1))
317 gerepileall(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
318 A11 = matslice(A, 1, m1, 1, n1);
319 B11 = matslice(B, 1, n1, 1, p1);
320 M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
321 if (gc_needed(av, 1))
322 gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
323 A12 = matslice(A, 1, m1, n1 + 1, n);
324 B21 = matslice(B, n1 + 1, n, 1, p1);
325 M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
326 if (gc_needed(av, 1))
327 gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
328 C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
329 if (gc_needed(av, 1))
330 gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11); /* destroy M4 */
331 M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
332 S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
333 if (gc_needed(av, 1))
334 gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4); /* destroy S3 */
335 T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
336 if (gc_needed(av, 1))
337 gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4); /* destroy T3 */
338 V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
339 if (gc_needed(av, 1))
340 gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1); /* destroy M1, M5 */
341 B22 = matslice(B, n1 + 1, n, p1 + 1, p);
342 M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
343 if (gc_needed(av, 1))
344 gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6); /* destroy S4, B22 */
345 A22 = matslice(A, m1 + 1, m, n1 + 1, n);
346 M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
347 if (gc_needed(av, 1))
348 gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7); /* destroy A22, T4 */
349 V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
350 C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
351 if (gc_needed(av, 1))
352 gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12); /* destroy V3, M6 */
353 V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
354 if (gc_needed(av, 1))
355 gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2); /* destroy V1, M2 */
356 C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
357 if (gc_needed(av, 1))
358 gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21); /* destroy M7 */
359 C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
360 if (gc_needed(av, 1))
361 gerepileall(av, 4, &C11, &C12, &C21, &C22); /* destroy V2, M3 */
362 C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
363 return gerepilecopy(av, C);
364 }
365
366 /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
367 static GEN
ZMrow_ZC_mul_i(GEN x,GEN y,long i,long lx)368 ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
369 {
370 pari_sp av = avma;
371 GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
372 long k;
373 for (k = 2; k < lx; k++)
374 {
375 GEN t = mulii(gcoeff(x,i,k), gel(y,k));
376 if (t != ZERO) c = addii(c, t);
377 }
378 return gerepileuptoint(av, c);
379 }
380 GEN
ZMrow_ZC_mul(GEN x,GEN y,long i)381 ZMrow_ZC_mul(GEN x, GEN y, long i)
382 { return ZMrow_ZC_mul_i(x, y, i, lg(x)); }
383
384 /* return x * y, 1 < lx = lg(x), l = lgcols(x) */
385 static GEN
ZM_ZC_mul_i(GEN x,GEN y,long lx,long l)386 ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
387 {
388 GEN z = cgetg(l,t_COL);
389 long i;
390 for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
391 return z;
392 }
393
394 static GEN
ZM_mul_classical(GEN x,GEN y,long l,long lx,long ly)395 ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
396 {
397 long j;
398 GEN z = cgetg(ly, t_MAT);
399 for (j = 1; j < ly; j++)
400 gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
401 return z;
402 }
403
404 static GEN
ZM_mul_slice(GEN A,GEN B,GEN P,GEN * mod)405 ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
406 {
407 pari_sp av = avma;
408 long i, n = lg(P)-1;
409 GEN H, T;
410 if (n == 1)
411 {
412 ulong p = uel(P,1);
413 GEN a = ZM_to_Flm(A, p);
414 GEN b = ZM_to_Flm(B, p);
415 GEN Hp = gerepileupto(av, Flm_to_ZM(Flm_mul(a, b, p)));
416 *mod = utoi(p); return Hp;
417 }
418 T = ZV_producttree(P);
419 A = ZM_nv_mod_tree(A, P, T);
420 B = ZM_nv_mod_tree(B, P, T);
421 H = cgetg(n+1, t_VEC);
422 for(i=1; i <= n; i++)
423 gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
424 H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
425 *mod = gmael(T, lg(T)-1, 1);
426 gerepileall(av, 2, &H, mod);
427 return H;
428 }
429
430 GEN
ZM_mul_worker(GEN P,GEN A,GEN B)431 ZM_mul_worker(GEN P, GEN A, GEN B)
432 {
433 GEN V = cgetg(3, t_VEC);
434 gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
435 return V;
436 }
437
438 static long
ZC_expi(GEN x)439 ZC_expi(GEN x)
440 {
441 long i, l = lg(x), m=-1;
442 for(i = 1; i < l; i++)
443 {
444 GEN xi = gel(x,i);
445 long e = signe(xi) ? expi(xi): -1;
446 if (e > m) m = e;
447 }
448 return m;
449 }
450
451 static long
ZM_expi(GEN x)452 ZM_expi(GEN x)
453 {
454 long i, l = lg(x), m=-1;
455 for(i = 1; i < l; i++)
456 {
457 long e = ZC_expi(gel(x,i));
458 if (e > m) m = e;
459 }
460 return m;
461 }
462
463 static GEN
ZM_mul_fast(GEN A,GEN B,long lA,long lB)464 ZM_mul_fast(GEN A, GEN B, long lA, long lB)
465 {
466 pari_sp av = avma;
467 forprime_t S;
468 GEN H, worker;
469 long h, mA = ZM_expi(A), mB = ZM_expi(B);
470 if (mA==-1 || mB==-1) return zeromat(nbrows(A),lB-1);
471 h = 3 + mA + mB + expu(lA-1);
472 init_modular_big(&S);
473 worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
474 H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
475 nmV_chinese_center, FpM_center);
476 return gerepileupto(av, H);
477 }
478
479 /* Strassen-Winograd used for dim >= ZM_sw_bound */
480 static GEN
ZM_mul_i(GEN x,GEN y,long l,long lx,long ly)481 ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
482 {
483 if (lx > 70 && ly > 70 && l > 70)
484 return ZM_mul_fast(x, y, lx, ly);
485 else
486 {
487 long s = maxss(ZM_max_lg_i(x,lx,l), ZM_max_lg_i(y,ly,lx));
488 long ZM_sw_bound = s > 60 ? 2: s > 25 ? 4: s>15 ? 8 : s > 8 ? 16 : 32;
489 if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
490 return ZM_mul_classical(x, y, l, lx, ly);
491 else
492 return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
493 }
494 }
495
496 GEN
ZM_mul(GEN x,GEN y)497 ZM_mul(GEN x, GEN y)
498 {
499 long lx=lg(x), ly=lg(y);
500 if (ly==1) return cgetg(1,t_MAT);
501 if (lx==1) return zeromat(0, ly-1);
502 return ZM_mul_i(x, y, lgcols(x), lx, ly);
503 }
504
505 GEN
QM_mul(GEN x,GEN y)506 QM_mul(GEN x, GEN y)
507 {
508 GEN dx, nx = Q_primitive_part(x, &dx);
509 GEN dy, ny = Q_primitive_part(y, &dy);
510 GEN z = ZM_mul(nx, ny);
511 if (dx || dy)
512 {
513 GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
514 if (!gequal1(d)) z = ZM_Q_mul(z, d);
515 }
516 return z;
517 }
518
519 GEN
QM_sqr(GEN x)520 QM_sqr(GEN x)
521 {
522 GEN dx, nx = Q_primitive_part(x, &dx);
523 GEN z = ZM_sqr(nx);
524 if (dx)
525 z = ZM_Q_mul(z, gsqr(dx));
526 return z;
527 }
528
529 GEN
QM_QC_mul(GEN x,GEN y)530 QM_QC_mul(GEN x, GEN y)
531 {
532 GEN dx, nx = Q_primitive_part(x, &dx);
533 GEN dy, ny = Q_primitive_part(y, &dy);
534 GEN z = ZM_ZC_mul(nx, ny);
535 if (dx || dy)
536 {
537 GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
538 if (!gequal1(d)) z = ZC_Q_mul(z, d);
539 }
540 return z;
541 }
542
543 /* assume result is symmetric */
544 GEN
ZM_multosym(GEN x,GEN y)545 ZM_multosym(GEN x, GEN y)
546 {
547 long j, lx, ly = lg(y);
548 GEN M;
549 if (ly == 1) return cgetg(1,t_MAT);
550 lx = lg(x); /* = lgcols(y) */
551 if (lx == 1) return cgetg(1,t_MAT);
552 /* ly = lgcols(x) */
553 M = cgetg(ly, t_MAT);
554 for (j=1; j<ly; j++)
555 {
556 GEN z = cgetg(ly,t_COL), yj = gel(y,j);
557 long i;
558 for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
559 for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
560 gel(M,j) = z;
561 }
562 return M;
563 }
564
565 /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
566 GEN
ZM_mul_diag(GEN m,GEN d)567 ZM_mul_diag(GEN m, GEN d)
568 {
569 long j, l;
570 GEN y = cgetg_copy(m, &l);
571 for (j=1; j<l; j++)
572 {
573 GEN c = gel(d,j);
574 gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
575 }
576 return y;
577 }
578 /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
579 GEN
ZM_diag_mul(GEN d,GEN m)580 ZM_diag_mul(GEN d, GEN m)
581 {
582 long i, j, l = lg(d), lm = lg(m);
583 GEN y = cgetg(lm, t_MAT);
584 for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
585 for (i=1; i<l; i++)
586 {
587 GEN c = gel(d,i);
588 if (equali1(c))
589 for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
590 else
591 for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
592 }
593 return y;
594 }
595
596 /* assume lx > 1 is lg(x) = lg(y) */
597 static GEN
ZV_dotproduct_i(GEN x,GEN y,long lx)598 ZV_dotproduct_i(GEN x, GEN y, long lx)
599 {
600 pari_sp av = avma;
601 GEN c = mulii(gel(x,1), gel(y,1));
602 long i;
603 for (i = 2; i < lx; i++)
604 {
605 GEN t = mulii(gel(x,i), gel(y,i));
606 if (t != gen_0) c = addii(c, t);
607 }
608 return gerepileuptoint(av, c);
609 }
610
611 /* x~ * y, assuming result is symmetric */
612 GEN
ZM_transmultosym(GEN x,GEN y)613 ZM_transmultosym(GEN x, GEN y)
614 {
615 long i, j, l, ly = lg(y);
616 GEN M;
617 if (ly == 1) return cgetg(1,t_MAT);
618 /* lg(x) = ly */
619 l = lgcols(y); /* = lgcols(x) */
620 M = cgetg(ly, t_MAT);
621 for (i=1; i<ly; i++)
622 {
623 GEN xi = gel(x,i), c = cgetg(ly,t_COL);
624 gel(M,i) = c;
625 for (j=1; j<i; j++)
626 gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
627 gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
628 }
629 return M;
630 }
631 /* x~ * y */
632 GEN
ZM_transmul(GEN x,GEN y)633 ZM_transmul(GEN x, GEN y)
634 {
635 long i, j, l, lx, ly = lg(y);
636 GEN M;
637 if (ly == 1) return cgetg(1,t_MAT);
638 lx = lg(x);
639 l = lgcols(y);
640 if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
641 M = cgetg(ly, t_MAT);
642 for (i=1; i<ly; i++)
643 {
644 GEN yi = gel(y,i), c = cgetg(lx,t_COL);
645 gel(M,i) = c;
646 for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
647 }
648 return M;
649 }
650
651 static GEN
ZM_sqr_i(GEN x,long l)652 ZM_sqr_i(GEN x, long l)
653 {
654 if (l > 70)
655 return ZM_mul_fast(x, x, l, l);
656 else
657 {
658 long s = ZM_max_lg_i(x,l,l);
659 long ZM_sw_bound = s > 60 ? 2: s > 25 ? 4: s>15 ? 8 : s > 8 ? 16 : 32;
660 if (l <= ZM_sw_bound)
661 return ZM_mul_classical(x, x, l, l, l);
662 else
663 return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
664 }
665 }
666
667 GEN
ZM_sqr(GEN x)668 ZM_sqr(GEN x)
669 {
670 long lx=lg(x);
671 if (lx==1) return cgetg(1,t_MAT);
672 return ZM_sqr_i(x, lx);
673 }
674 GEN
ZM_ZC_mul(GEN x,GEN y)675 ZM_ZC_mul(GEN x, GEN y)
676 {
677 long lx = lg(x);
678 return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
679 }
680
681 GEN
ZC_Z_div(GEN x,GEN c)682 ZC_Z_div(GEN x, GEN c)
683 { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
684
685 GEN
ZM_Z_div(GEN x,GEN c)686 ZM_Z_div(GEN x, GEN c)
687 { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
688
689 GEN
ZC_Q_mul(GEN A,GEN z)690 ZC_Q_mul(GEN A, GEN z)
691 {
692 pari_sp av = avma;
693 long i, l = lg(A);
694 GEN d, n, Ad, B, u;
695 if (typ(z)==t_INT) return ZC_Z_mul(A,z);
696 n = gel(z, 1); d = gel(z, 2);
697 Ad = FpC_red(A, d);
698 u = gcdii(d, FpV_factorback(Ad, NULL, d));
699 B = cgetg(l, t_COL);
700 if (equali1(u))
701 {
702 for(i=1; i<l; i++)
703 gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
704 } else
705 {
706 for(i=1; i<l; i++)
707 {
708 GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
709 if (equalii(d, di))
710 gel(B, i) = ni;
711 else
712 gel(B, i) = mkfrac(ni, diviiexact(d, di));
713 }
714 }
715 return gerepilecopy(av, B);
716 }
717
718 GEN
ZM_Q_mul(GEN x,GEN z)719 ZM_Q_mul(GEN x, GEN z)
720 {
721 if (typ(z)==t_INT) return ZM_Z_mul(x,z);
722 pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
723 }
724
725 long
zv_dotproduct(GEN x,GEN y)726 zv_dotproduct(GEN x, GEN y)
727 {
728 long i, lx = lg(x);
729 ulong c;
730 if (lx == 1) return 0;
731 c = uel(x,1)*uel(y,1);
732 for (i = 2; i < lx; i++)
733 c += uel(x,i)*uel(y,i);
734 return c;
735 }
736
737 GEN
ZV_ZM_mul(GEN x,GEN y)738 ZV_ZM_mul(GEN x, GEN y)
739 {
740 long i, lx = lg(x), ly = lg(y);
741 GEN z;
742 if (lx == 1) return zerovec(ly-1);
743 z = cgetg(ly, t_VEC);
744 for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
745 return z;
746 }
747
748 GEN
ZC_ZV_mul(GEN x,GEN y)749 ZC_ZV_mul(GEN x, GEN y)
750 {
751 long i,j, lx=lg(x), ly=lg(y);
752 GEN z;
753 if (ly==1) return cgetg(1,t_MAT);
754 z = cgetg(ly,t_MAT);
755 for (j=1; j < ly; j++)
756 {
757 gel(z,j) = cgetg(lx,t_COL);
758 for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
759 }
760 return z;
761 }
762
763 GEN
ZV_dotsquare(GEN x)764 ZV_dotsquare(GEN x)
765 {
766 long i, lx;
767 pari_sp av;
768 GEN z;
769 lx = lg(x);
770 if (lx == 1) return gen_0;
771 av = avma; z = sqri(gel(x,1));
772 for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
773 return gerepileuptoint(av,z);
774 }
775
776 GEN
ZV_dotproduct(GEN x,GEN y)777 ZV_dotproduct(GEN x,GEN y)
778 {
779 long lx;
780 if (x == y) return ZV_dotsquare(x);
781 lx = lg(x);
782 if (lx == 1) return gen_0;
783 return ZV_dotproduct_i(x, y, lx);
784 }
785
786 static GEN
_ZM_mul(void * data,GEN x,GEN y)787 _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
788 { (void)data; return ZM_mul(x,y); }
789 static GEN
_ZM_sqr(void * data,GEN x)790 _ZM_sqr(void *data /*ignored*/, GEN x)
791 { (void)data; return ZM_sqr(x); }
792 GEN
ZM_pow(GEN x,GEN n)793 ZM_pow(GEN x, GEN n)
794 {
795 pari_sp av = avma;
796 if (!signe(n)) return matid(lg(x)-1);
797 return gerepilecopy(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
798 }
799 GEN
ZM_powu(GEN x,ulong n)800 ZM_powu(GEN x, ulong n)
801 {
802 pari_sp av = avma;
803 if (!n) return matid(lg(x)-1);
804 return gerepilecopy(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
805 }
806 /********************************************************************/
807 /** **/
808 /** ADD, SUB **/
809 /** **/
810 /********************************************************************/
811 static GEN
ZC_add_i(GEN x,GEN y,long lx)812 ZC_add_i(GEN x, GEN y, long lx)
813 {
814 GEN A = cgetg(lx, t_COL);
815 long i;
816 for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
817 return A;
818 }
819 GEN
ZC_add(GEN x,GEN y)820 ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
821 GEN
ZC_Z_add(GEN x,GEN y)822 ZC_Z_add(GEN x, GEN y)
823 {
824 long k, lx = lg(x);
825 GEN z = cgetg(lx, t_COL);
826 if (lx == 1) pari_err_TYPE2("+",x,y);
827 gel(z,1) = addii(y,gel(x,1));
828 for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
829 return z;
830 }
831
832 static GEN
ZC_sub_i(GEN x,GEN y,long lx)833 ZC_sub_i(GEN x, GEN y, long lx)
834 {
835 long i;
836 GEN A = cgetg(lx, t_COL);
837 for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
838 return A;
839 }
840 GEN
ZC_sub(GEN x,GEN y)841 ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
842 GEN
ZC_Z_sub(GEN x,GEN y)843 ZC_Z_sub(GEN x, GEN y)
844 {
845 long k, lx = lg(x);
846 GEN z = cgetg(lx, t_COL);
847 if (lx == 1) pari_err_TYPE2("+",x,y);
848 gel(z,1) = subii(gel(x,1), y);
849 for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
850 return z;
851 }
852 GEN
Z_ZC_sub(GEN a,GEN x)853 Z_ZC_sub(GEN a, GEN x)
854 {
855 long k, lx = lg(x);
856 GEN z = cgetg(lx, t_COL);
857 if (lx == 1) pari_err_TYPE2("-",a,x);
858 gel(z,1) = subii(a, gel(x,1));
859 for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
860 return z;
861 }
862
863 GEN
ZM_add(GEN x,GEN y)864 ZM_add(GEN x, GEN y)
865 {
866 long lx = lg(x), l, j;
867 GEN z;
868 if (lx == 1) return cgetg(1, t_MAT);
869 z = cgetg(lx, t_MAT); l = lgcols(x);
870 for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
871 return z;
872 }
873 GEN
ZM_sub(GEN x,GEN y)874 ZM_sub(GEN x, GEN y)
875 {
876 long lx = lg(x), l, j;
877 GEN z;
878 if (lx == 1) return cgetg(1, t_MAT);
879 z = cgetg(lx, t_MAT); l = lgcols(x);
880 for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
881 return z;
882 }
883 /********************************************************************/
884 /** **/
885 /** LINEAR COMBINATION **/
886 /** **/
887 /********************************************************************/
888 /* return X/c assuming division is exact */
889 GEN
ZC_Z_divexact(GEN x,GEN c)890 ZC_Z_divexact(GEN x, GEN c)
891 { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
892 GEN
ZC_u_divexact(GEN x,ulong c)893 ZC_u_divexact(GEN x, ulong c)
894 { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
895
896 GEN
ZM_Z_divexact(GEN x,GEN c)897 ZM_Z_divexact(GEN x, GEN c)
898 { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
899
900 GEN
ZC_Z_mul(GEN x,GEN c)901 ZC_Z_mul(GEN x, GEN c)
902 {
903 if (!signe(c)) return zerocol(lg(x)-1);
904 if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
905 pari_APPLY_type(t_COL, mulii(gel(x,i), c))
906 }
907
908 GEN
ZC_z_mul(GEN x,long c)909 ZC_z_mul(GEN x, long c)
910 {
911 if (!c) return zerocol(lg(x)-1);
912 if (c == 1) return ZC_copy(x);
913 if (c ==-1) return ZC_neg(x);
914 pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
915 }
916
917 GEN
zv_z_mul(GEN x,long n)918 zv_z_mul(GEN x, long n)
919 { pari_APPLY_long(x[i]*n) }
920
921 /* return a ZM */
922 GEN
nm_Z_mul(GEN X,GEN c)923 nm_Z_mul(GEN X, GEN c)
924 {
925 long i, j, h, l = lg(X), s = signe(c);
926 GEN A;
927 if (l == 1) return cgetg(1, t_MAT);
928 h = lgcols(X);
929 if (!s) return zeromat(h-1, l-1);
930 if (is_pm1(c)) {
931 if (s > 0) return Flm_to_ZM(X);
932 X = Flm_to_ZM(X); ZM_togglesign(X); return X;
933 }
934 A = cgetg(l, t_MAT);
935 for (j = 1; j < l; j++)
936 {
937 GEN a = cgetg(h, t_COL), x = gel(X, j);
938 for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
939 gel(A,j) = a;
940 }
941 return A;
942 }
943 GEN
ZM_Z_mul(GEN X,GEN c)944 ZM_Z_mul(GEN X, GEN c)
945 {
946 long i, j, h, l = lg(X);
947 GEN A;
948 if (l == 1) return cgetg(1, t_MAT);
949 h = lgcols(X);
950 if (!signe(c)) return zeromat(h-1, l-1);
951 if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
952 A = cgetg(l, t_MAT);
953 for (j = 1; j < l; j++)
954 {
955 GEN a = cgetg(h, t_COL), x = gel(X, j);
956 for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
957 gel(A,j) = a;
958 }
959 return A;
960 }
961 void
ZC_lincomb1_inplace_i(GEN X,GEN Y,GEN v,long n)962 ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
963 {
964 long i;
965 for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
966 }
967 /* X <- X + v Y (elementary col operation) */
968 void
ZC_lincomb1_inplace(GEN X,GEN Y,GEN v)969 ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
970 {
971 if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
972 }
973 void
Flc_lincomb1_inplace(GEN X,GEN Y,ulong v,ulong q)974 Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
975 {
976 long i;
977 if (!v) return; /* v = 0 */
978 for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
979 }
980
981 /* X + v Y, wasteful if (v = 0) */
982 static GEN
ZC_lincomb1(GEN v,GEN x,GEN y)983 ZC_lincomb1(GEN v, GEN x, GEN y)
984 { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
985
986 /* -X + vY */
987 static GEN
ZC_lincomb_1(GEN v,GEN x,GEN y)988 ZC_lincomb_1(GEN v, GEN x, GEN y)
989 { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
990
991 /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
992 GEN
ZC_lincomb(GEN u,GEN v,GEN X,GEN Y)993 ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
994 {
995 long su, sv;
996 GEN A;
997
998 su = signe(u); if (!su) return ZC_Z_mul(Y, v);
999 sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
1000 if (is_pm1(v))
1001 {
1002 if (is_pm1(u))
1003 {
1004 if (su != sv) A = ZC_sub(X, Y);
1005 else A = ZC_add(X, Y);
1006 if (su < 0) ZV_togglesign(A); /* in place but was created above */
1007 }
1008 else
1009 {
1010 if (sv > 0) A = ZC_lincomb1 (u, Y, X);
1011 else A = ZC_lincomb_1(u, Y, X);
1012 }
1013 }
1014 else if (is_pm1(u))
1015 {
1016 if (su > 0) A = ZC_lincomb1 (v, X, Y);
1017 else A = ZC_lincomb_1(v, X, Y);
1018 }
1019 else
1020 { /* not cgetg_copy: x may be a t_VEC */
1021 long i, lx = lg(X);
1022 A = cgetg(lx,t_COL);
1023 for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
1024 }
1025 return A;
1026 }
1027
1028 /********************************************************************/
1029 /** **/
1030 /** CONVERSIONS **/
1031 /** **/
1032 /********************************************************************/
1033 GEN
ZV_to_nv(GEN x)1034 ZV_to_nv(GEN x)
1035 { pari_APPLY_ulong(itou(gel(x,i))) }
1036
1037 GEN
zm_to_ZM(GEN x)1038 zm_to_ZM(GEN x)
1039 { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
1040
1041 GEN
zmV_to_ZMV(GEN x)1042 zmV_to_ZMV(GEN x)
1043 { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
1044
1045 /* same as Flm_to_ZM but do not assume positivity */
1046 GEN
ZM_to_zm(GEN x)1047 ZM_to_zm(GEN x)
1048 { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
1049
1050 GEN
zv_to_Flv(GEN x,ulong p)1051 zv_to_Flv(GEN x, ulong p)
1052 { pari_APPLY_ulong(umodsu(x[i], p)) }
1053
1054 GEN
zm_to_Flm(GEN x,ulong p)1055 zm_to_Flm(GEN x, ulong p)
1056 { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
1057
1058 GEN
ZMV_to_zmV(GEN x)1059 ZMV_to_zmV(GEN x)
1060 { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
1061
1062 /********************************************************************/
1063 /** **/
1064 /** COPY, NEGATION **/
1065 /** **/
1066 /********************************************************************/
1067 GEN
ZC_copy(GEN x)1068 ZC_copy(GEN x)
1069 {
1070 long i, lx = lg(x);
1071 GEN y = cgetg(lx, t_COL);
1072 for (i=1; i<lx; i++)
1073 {
1074 GEN c = gel(x,i);
1075 gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
1076 }
1077 return y;
1078 }
1079
1080 GEN
ZM_copy(GEN x)1081 ZM_copy(GEN x)
1082 { pari_APPLY_same(ZC_copy(gel(x,i))) }
1083
1084 void
ZV_neg_inplace(GEN M)1085 ZV_neg_inplace(GEN M)
1086 {
1087 long l = lg(M);
1088 while (--l > 0) gel(M,l) = negi(gel(M,l));
1089 }
1090 GEN
ZC_neg(GEN x)1091 ZC_neg(GEN x)
1092 { pari_APPLY_type(t_COL, negi(gel(x,i))) }
1093 GEN
zv_neg(GEN x)1094 zv_neg(GEN x)
1095 { pari_APPLY_long(-x[i]) }
1096 GEN
zv_neg_inplace(GEN M)1097 zv_neg_inplace(GEN M)
1098 {
1099 long l = lg(M);
1100 while (--l > 0) M[l] = -M[l];
1101 return M;
1102 }
1103 GEN
ZM_neg(GEN x)1104 ZM_neg(GEN x)
1105 { pari_APPLY_same(ZC_neg(gel(x,i))) }
1106
1107 void
ZV_togglesign(GEN M)1108 ZV_togglesign(GEN M)
1109 {
1110 long l = lg(M);
1111 while (--l > 0) togglesign_safe(&gel(M,l));
1112 }
1113 void
ZM_togglesign(GEN M)1114 ZM_togglesign(GEN M)
1115 {
1116 long l = lg(M);
1117 while (--l > 0) ZV_togglesign(gel(M,l));
1118 }
1119
1120 /********************************************************************/
1121 /** **/
1122 /** "DIVISION" mod HNF **/
1123 /** **/
1124 /********************************************************************/
1125 /* Reduce ZC x modulo ZM y in HNF, may return x itself (not a copy) */
1126 GEN
ZC_hnfremdiv(GEN x,GEN y,GEN * Q)1127 ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
1128 {
1129 long i, l = lg(x);
1130 GEN q;
1131
1132 if (Q) *Q = cgetg(l,t_COL);
1133 if (l == 1) return cgetg(1,t_COL);
1134 for (i = l-1; i>0; i--)
1135 {
1136 q = diviiround(gel(x,i), gcoeff(y,i,i));
1137 if (signe(q)) {
1138 togglesign(q);
1139 x = ZC_lincomb(gen_1, q, x, gel(y,i));
1140 }
1141 if (Q) gel(*Q, i) = q;
1142 }
1143 return x;
1144 }
1145
1146 /* x = y Q + R, may return some columns of x (not copies) */
1147 GEN
ZM_hnfdivrem(GEN x,GEN y,GEN * Q)1148 ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
1149 {
1150 long lx = lg(x), i;
1151 GEN R = cgetg(lx, t_MAT);
1152 if (Q)
1153 {
1154 GEN q = cgetg(lx, t_MAT); *Q = q;
1155 for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
1156 }
1157 else
1158 for (i=1; i<lx; i++)
1159 {
1160 pari_sp av = avma;
1161 GEN z = ZC_hnfrem(gel(x,i),y);
1162 gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
1163 }
1164 return R;
1165 }
1166
1167 /********************************************************************/
1168 /** **/
1169 /** TESTS **/
1170 /** **/
1171 /********************************************************************/
1172 int
zv_equal0(GEN V)1173 zv_equal0(GEN V)
1174 {
1175 long l = lg(V);
1176 while (--l > 0)
1177 if (V[l]) return 0;
1178 return 1;
1179 }
1180
1181 int
ZV_equal0(GEN V)1182 ZV_equal0(GEN V)
1183 {
1184 long l = lg(V);
1185 while (--l > 0)
1186 if (signe(gel(V,l))) return 0;
1187 return 1;
1188 }
1189 int
ZMrow_equal0(GEN V,long i)1190 ZMrow_equal0(GEN V, long i)
1191 {
1192 long l = lg(V);
1193 while (--l > 0)
1194 if (signe(gcoeff(V,i,l))) return 0;
1195 return 1;
1196 }
1197
1198 static int
ZV_equal_lg(GEN V,GEN W,long l)1199 ZV_equal_lg(GEN V, GEN W, long l)
1200 {
1201 while (--l > 0)
1202 if (!equalii(gel(V,l), gel(W,l))) return 0;
1203 return 1;
1204 }
1205 int
ZV_equal(GEN V,GEN W)1206 ZV_equal(GEN V, GEN W)
1207 {
1208 long l = lg(V);
1209 if (lg(W) != l) return 0;
1210 return ZV_equal_lg(V, W, l);
1211 }
1212 int
ZM_equal(GEN A,GEN B)1213 ZM_equal(GEN A, GEN B)
1214 {
1215 long i, m, l = lg(A);
1216 if (lg(B) != l) return 0;
1217 if (l == 1) return 1;
1218 m = lgcols(A);
1219 if (lgcols(B) != m) return 0;
1220 for (i = 1; i < l; i++)
1221 if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
1222 return 1;
1223 }
1224 int
ZM_equal0(GEN A)1225 ZM_equal0(GEN A)
1226 {
1227 long i, j, m, l = lg(A);
1228 if (l == 1) return 1;
1229 m = lgcols(A);
1230 for (j = 1; j < l; j++)
1231 for (i = 1; i < m; i++)
1232 if (signe(gcoeff(A,i,j))) return 0;
1233 return 1;
1234 }
1235 int
zv_equal(GEN V,GEN W)1236 zv_equal(GEN V, GEN W)
1237 {
1238 long l = lg(V);
1239 if (lg(W) != l) return 0;
1240 while (--l > 0)
1241 if (V[l] != W[l]) return 0;
1242 return 1;
1243 }
1244
1245 int
zvV_equal(GEN V,GEN W)1246 zvV_equal(GEN V, GEN W)
1247 {
1248 long l = lg(V);
1249 if (lg(W) != l) return 0;
1250 while (--l > 0)
1251 if (!zv_equal(gel(V,l),gel(W,l))) return 0;
1252 return 1;
1253 }
1254
1255 int
ZM_ishnf(GEN x)1256 ZM_ishnf(GEN x)
1257 {
1258 long i,j, lx = lg(x);
1259 for (i=1; i<lx; i++)
1260 {
1261 GEN xii = gcoeff(x,i,i);
1262 if (signe(xii) <= 0) return 0;
1263 for (j=1; j<i; j++)
1264 if (signe(gcoeff(x,i,j))) return 0;
1265 for (j=i+1; j<lx; j++)
1266 {
1267 GEN xij = gcoeff(x,i,j);
1268 if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
1269 }
1270 }
1271 return 1;
1272 }
1273 int
ZM_isidentity(GEN x)1274 ZM_isidentity(GEN x)
1275 {
1276 long i,j, lx = lg(x);
1277
1278 if (lx == 1) return 1;
1279 if (lx != lgcols(x)) return 0;
1280 for (j=1; j<lx; j++)
1281 {
1282 GEN c = gel(x,j);
1283 for (i=1; i<j; )
1284 if (signe(gel(c,i++))) return 0;
1285 /* i = j */
1286 if (!equali1(gel(c,i++))) return 0;
1287 for ( ; i<lx; )
1288 if (signe(gel(c,i++))) return 0;
1289 }
1290 return 1;
1291 }
1292 int
ZM_isdiagonal(GEN x)1293 ZM_isdiagonal(GEN x)
1294 {
1295 long i,j, lx = lg(x);
1296 if (lx == 1) return 1;
1297 if (lx != lgcols(x)) return 0;
1298
1299 for (j=1; j<lx; j++)
1300 {
1301 GEN c = gel(x,j);
1302 for (i=1; i<j; i++)
1303 if (signe(gel(c,i))) return 0;
1304 for (i++; i<lx; i++)
1305 if (signe(gel(c,i))) return 0;
1306 }
1307 return 1;
1308 }
1309 int
ZM_isscalar(GEN x,GEN s)1310 ZM_isscalar(GEN x, GEN s)
1311 {
1312 long i, j, lx = lg(x);
1313
1314 if (lx == 1) return 1;
1315 if (!s) s = gcoeff(x,1,1);
1316 if (equali1(s)) return ZM_isidentity(x);
1317 if (lx != lgcols(x)) return 0;
1318 for (j=1; j<lx; j++)
1319 {
1320 GEN c = gel(x,j);
1321 for (i=1; i<j; )
1322 if (signe(gel(c,i++))) return 0;
1323 /* i = j */
1324 if (!equalii(gel(c,i++), s)) return 0;
1325 for ( ; i<lx; )
1326 if (signe(gel(c,i++))) return 0;
1327 }
1328 return 1;
1329 }
1330
1331 long
ZC_is_ei(GEN x)1332 ZC_is_ei(GEN x)
1333 {
1334 long i, j = 0, l = lg(x);
1335 for (i = 1; i < l; i++)
1336 {
1337 GEN c = gel(x,i);
1338 long s = signe(c);
1339 if (!s) continue;
1340 if (s < 0 || !is_pm1(c) || j) return 0;
1341 j = i;
1342 }
1343 return j;
1344 }
1345
1346 /********************************************************************/
1347 /** **/
1348 /** MISCELLANEOUS **/
1349 /** **/
1350 /********************************************************************/
1351 /* assume lg(x) = lg(y), x,y in Z^n */
1352 int
ZV_cmp(GEN x,GEN y)1353 ZV_cmp(GEN x, GEN y)
1354 {
1355 long fl,i, lx = lg(x);
1356 for (i=1; i<lx; i++)
1357 if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
1358 return 0;
1359 }
1360 /* assume lg(x) = lg(y), x,y in Z^n */
1361 int
ZV_abscmp(GEN x,GEN y)1362 ZV_abscmp(GEN x, GEN y)
1363 {
1364 long fl,i, lx = lg(x);
1365 for (i=1; i<lx; i++)
1366 if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
1367 return 0;
1368 }
1369
1370 long
zv_content(GEN x)1371 zv_content(GEN x)
1372 {
1373 long i, s, l = lg(x);
1374 if (l == 1) return 0;
1375 s = labs(x[1]);
1376 for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
1377 return s;
1378 }
1379 GEN
ZV_content(GEN x)1380 ZV_content(GEN x)
1381 {
1382 long i, l = lg(x);
1383 pari_sp av = avma;
1384 GEN c;
1385 if (l == 1) return gen_0;
1386 if (l == 2) return absi(gel(x,1));
1387 c = gel(x,1);
1388 for (i = 2; i < l; i++)
1389 {
1390 c = gcdii(c, gel(x,i));
1391 if (is_pm1(c)) { set_avma(av); return gen_1; }
1392 }
1393 return gerepileuptoint(av, c);
1394 }
1395
1396 GEN
ZM_det_triangular(GEN mat)1397 ZM_det_triangular(GEN mat)
1398 {
1399 pari_sp av;
1400 long i,l = lg(mat);
1401 GEN s;
1402
1403 if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
1404 av = avma; s = gcoeff(mat,1,1);
1405 for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
1406 return gerepileuptoint(av,s);
1407 }
1408
1409 /* assumes no overflow */
1410 long
zv_prod(GEN v)1411 zv_prod(GEN v)
1412 {
1413 long n, i, l = lg(v);
1414 if (l == 1) return 1;
1415 n = v[1]; for (i = 2; i < l; i++) n *= v[i];
1416 return n;
1417 }
1418
1419 static GEN
_mulii(void * E,GEN a,GEN b)1420 _mulii(void *E, GEN a, GEN b)
1421 { (void) E; return mulii(a, b); }
1422
1423 /* product of ulongs */
1424 GEN
zv_prod_Z(GEN v)1425 zv_prod_Z(GEN v)
1426 {
1427 pari_sp av = avma;
1428 long k, m, n = lg(v)-1;
1429 GEN V;
1430 switch(n) {
1431 case 0: return gen_1;
1432 case 1: return utoi(v[1]);
1433 case 2: return muluu(v[1], v[2]);
1434 }
1435 m = n >> 1;
1436 V = cgetg(m + (odd(n)? 2: 1), t_VEC);
1437 for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
1438 if (odd(n)) gel(V,k) = utoipos(v[n]);
1439 return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
1440 }
1441 GEN
vecsmall_prod(GEN v)1442 vecsmall_prod(GEN v)
1443 {
1444 pari_sp av = avma;
1445 long k, m, n = lg(v)-1;
1446 GEN V;
1447 switch (n) {
1448 case 0: return gen_1;
1449 case 1: return stoi(v[1]);
1450 case 2: return mulss(v[1], v[2]);
1451 }
1452 m = n >> 1;
1453 V = cgetg(m + (odd(n)? 2: 1), t_VEC);
1454 for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
1455 if (odd(n)) gel(V,k) = stoi(v[n]);
1456 return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
1457 }
1458
1459 GEN
ZV_prod(GEN v)1460 ZV_prod(GEN v)
1461 {
1462 pari_sp av = avma;
1463 long i, l = lg(v);
1464 GEN n;
1465 if (l == 1) return gen_1;
1466 if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
1467 n = gel(v,1);
1468 if (l == 2) return icopy(n);
1469 for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
1470 return gerepileuptoint(av, n);
1471 }
1472 /* assumes no overflow */
1473 long
zv_sum(GEN v)1474 zv_sum(GEN v)
1475 {
1476 long n, i, l = lg(v);
1477 if (l == 1) return 0;
1478 n = v[1]; for (i = 2; i < l; i++) n += v[i];
1479 return n;
1480 }
1481 /* assumes no overflow and 0 <= n <= #v */
1482 long
zv_sumpart(GEN v,long n)1483 zv_sumpart(GEN v, long n)
1484 {
1485 long i, p;
1486 if (!n) return 0;
1487 p = v[1]; for (i = 2; i <= n; i++) p += v[i];
1488 return p;
1489 }
1490 GEN
ZV_sum(GEN v)1491 ZV_sum(GEN v)
1492 {
1493 pari_sp av = avma;
1494 long i, l = lg(v);
1495 GEN n;
1496 if (l == 1) return gen_0;
1497 n = gel(v,1);
1498 if (l == 2) return icopy(n);
1499 for (i = 2; i < l; i++) n = addii(n, gel(v,i));
1500 return gerepileuptoint(av, n);
1501 }
1502
1503 /********************************************************************/
1504 /** **/
1505 /** GRAM SCHMIDT REDUCTION (integer matrices) **/
1506 /** **/
1507 /********************************************************************/
1508
1509 /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
1510 static void
Zupdate_row(long k,long l,GEN q,GEN L,GEN B)1511 Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
1512 {
1513 long i, qq = itos_or_0(q);
1514 if (!qq)
1515 {
1516 for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
1517 gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
1518 return;
1519 }
1520 if (qq == 1) {
1521 for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
1522 gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
1523 } else if (qq == -1) {
1524 for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
1525 gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
1526 } else {
1527 for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
1528 gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
1529 }
1530 }
1531
1532 /* update L[k,] */
1533 static void
ZRED(long k,long l,GEN x,GEN L,GEN B)1534 ZRED(long k, long l, GEN x, GEN L, GEN B)
1535 {
1536 GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
1537 if (!signe(q)) return;
1538 q = negi(q);
1539 Zupdate_row(k,l,q,L,B);
1540 gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
1541 }
1542
1543 /* Gram-Schmidt reduction, x a ZM */
1544 static void
ZincrementalGS(GEN x,GEN L,GEN B,long k)1545 ZincrementalGS(GEN x, GEN L, GEN B, long k)
1546 {
1547 long i, j;
1548 for (j=1; j<=k; j++)
1549 {
1550 pari_sp av = avma;
1551 GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
1552 for (i=1; i<j; i++)
1553 {
1554 u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
1555 u = diviiexact(u, gel(B,i));
1556 }
1557 gcoeff(L,k,j) = gerepileuptoint(av, u);
1558 }
1559 gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
1560 }
1561
1562 /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
1563 * Very inefficient if y is not LLL-reduced of maximal rank */
1564 static GEN
ZC_reducemodmatrix_i(GEN v,GEN y)1565 ZC_reducemodmatrix_i(GEN v, GEN y)
1566 {
1567 GEN B, L, x = shallowconcat(y, v);
1568 long k, lx = lg(x), nx = lx-1;
1569
1570 B = scalarcol_shallow(gen_1, lx);
1571 L = zeromatcopy(nx, nx);
1572 for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
1573 for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1574 return gel(x,nx);
1575 }
1576 GEN
ZC_reducemodmatrix(GEN v,GEN y)1577 ZC_reducemodmatrix(GEN v, GEN y) {
1578 pari_sp av = avma;
1579 return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
1580 }
1581
1582 /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
1583 * Very inefficient if y is not LLL-reduced of maximal rank */
1584 static GEN
ZM_reducemodmatrix_i(GEN v,GEN y)1585 ZM_reducemodmatrix_i(GEN v, GEN y)
1586 {
1587 GEN B, L, V;
1588 long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
1589
1590 V = cgetg(lv, t_MAT);
1591 B = scalarcol_shallow(gen_1, lx);
1592 L = zeromatcopy(nx, nx);
1593 for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
1594 for (j = 1; j < lg(v); j++)
1595 {
1596 GEN x = shallowconcat(y, gel(v,j));
1597 ZincrementalGS(x, L, B, nx); /* overwrite last */
1598 for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1599 gel(V,j) = gel(x,nx);
1600 }
1601 return V;
1602 }
1603 GEN
ZM_reducemodmatrix(GEN v,GEN y)1604 ZM_reducemodmatrix(GEN v, GEN y) {
1605 pari_sp av = avma;
1606 return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
1607 }
1608
1609 GEN
ZC_reducemodlll(GEN x,GEN y)1610 ZC_reducemodlll(GEN x,GEN y)
1611 {
1612 pari_sp av = avma;
1613 GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
1614 return gerepilecopy(av, z);
1615 }
1616 GEN
ZM_reducemodlll(GEN x,GEN y)1617 ZM_reducemodlll(GEN x,GEN y)
1618 {
1619 pari_sp av = avma;
1620 GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
1621 return gerepilecopy(av, z);
1622 }
1623