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