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 1723195 : check_ZV(GEN x, long l)
20 : {
21 : long i;
22 10247414 : for (i=1; i<l; i++)
23 8524310 : if (typ(gel(x,i)) != t_INT) return 0;
24 1723104 : return 1;
25 : }
26 : void
27 1415239 : RgV_check_ZV(GEN A, const char *s)
28 : {
29 1415239 : if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
30 1415233 : }
31 : void
32 599031 : RgM_check_ZM(GEN A, const char *s)
33 : {
34 599031 : long n = lg(A);
35 599031 : if (n != 1)
36 : {
37 598884 : long j, m = lgcols(A);
38 2321988 : for (j=1; j<n; j++)
39 1723194 : if (!check_ZV(gel(A,j), m))
40 91 : pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
41 : }
42 598941 : }
43 :
44 : static long
45 68386521 : ZV_max_lg_i(GEN x, long m)
46 : {
47 68386521 : long i, prec = 2;
48 653955392 : for (i = 1; i < m; i++) prec = maxss(prec, lgefint(gel(x,i)));
49 68386898 : 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 17824259 : ZM_max_lg_i(GEN x, long n, long m)
56 : {
57 17824259 : long j, prec = 2;
58 86199805 : for (j = 1; j < n; j++) prec = maxss(prec, ZV_max_lg_i(gel(x,j), m));
59 17824580 : return prec;
60 : }
61 : long
62 1594571 : ZM_max_lg(GEN x)
63 : {
64 1594571 : long n = lg(x);
65 1594571 : if (n == 1) return 2;
66 1594571 : return ZM_max_lg_i(x, n, lgcols(x));
67 : }
68 :
69 : GEN
70 3784 : ZM_supnorm(GEN x)
71 : {
72 3784 : long i, j, h, lx = lg(x);
73 3784 : GEN s = gen_0;
74 3784 : if (lx == 1) return gen_1;
75 3784 : h = lgcols(x);
76 23211 : for (j=1; j<lx; j++)
77 : {
78 19427 : GEN xj = gel(x,j);
79 272114 : for (i=1; i<h; i++)
80 : {
81 252687 : GEN c = gel(xj,i);
82 252687 : if (abscmpii(c, s) > 0) s = c;
83 : }
84 : }
85 3784 : 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 23054836 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
126 : {
127 : long i, j;
128 23054836 : GEN z = cgetg(l,t_COL);
129 :
130 133470396 : for (i=1; i<l; i++)
131 : {
132 110418424 : pari_sp av = avma;
133 110418424 : GEN s = mulis(gcoeff(x,i,1),y[1]);
134 1195963634 : for (j=2; j<c; j++)
135 1085551083 : if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
136 110412551 : gel(z,i) = gerepileuptoint(av,s);
137 : }
138 23051972 : return z;
139 : }
140 : GEN
141 21344912 : ZM_zc_mul(GEN x, GEN y) {
142 21344912 : long l = lg(x);
143 21344912 : if (l == 1) return cgetg(1, t_COL);
144 21344912 : 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 806810 : ZM_zm_mul(GEN x, GEN y)
168 : {
169 806810 : long j, c, l = lg(x), ly = lg(y);
170 806810 : GEN z = cgetg(ly, t_MAT);
171 806810 : if (l == 1) return z;
172 806803 : c = lgcols(x);
173 2516946 : for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
174 806803 : 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 274352 : 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 274352 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
198 274352 : GEN M = cgetg(n + 1, t_MAT), C;
199 :
200 2607907 : for (j = 1; j <= min_e; j++) {
201 2333555 : gel(M, j) = C = cgetg(m + 1, t_COL);
202 41208305 : for (i = 1; i <= min_d; i++)
203 38874750 : gel(C, i) = addii(gcoeff(A, ma + i, na + j),
204 38874750 : gcoeff(B, mb + i, nb + j));
205 2362960 : for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
206 2333555 : for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
207 2333555 : for (; i <= m; i++) gel(C, i) = gen_0;
208 : }
209 293470 : for (; j <= ea; j++) {
210 19118 : gel(M, j) = C = cgetg(m + 1, t_COL);
211 83048 : for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
212 19118 : for (; i <= m; i++) gel(C, i) = gen_0;
213 : }
214 274352 : 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 274352 : for (; j <= n; j++) gel(M, j) = zerocol(m);
220 274352 : 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 240058 : 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 240058 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
231 240058 : GEN M = cgetg(n + 1, t_MAT), C;
232 :
233 2271095 : for (j = 1; j <= min_e; j++) {
234 2031037 : gel(M, j) = C = cgetg(m + 1, t_COL);
235 36710754 : for (i = 1; i <= min_d; i++)
236 34679717 : gel(C, i) = subii(gcoeff(A, ma + i, na + j),
237 34679717 : gcoeff(B, mb + i, nb + j));
238 2050251 : for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
239 2062896 : for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
240 2031037 : for (; i <= m; i++) gel(C, i) = gen_0;
241 : }
242 240058 : 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 254964 : for (; j <= eb; j++) {
248 14906 : gel(M, j) = C = cgetg(m + 1, t_COL);
249 46788 : for (i = 1; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
250 14906 : for (; i <= m; i++) gel(C, i) = gen_0;
251 : }
252 254964 : for (; j <= n; j++) gel(M, j) = zerocol(m);
253 240058 : 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 34294 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
261 : {
262 34294 : pari_sp av = avma;
263 34294 : long m1 = (m + 1)/2, m2 = m/2,
264 34294 : n1 = (n + 1)/2, n2 = n/2,
265 34294 : 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 34294 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
272 34294 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
273 34294 : M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
274 34294 : if (gc_needed(av, 1))
275 0 : gerepileall(av, 2, &T2, &M2); /* destroy S1 */
276 34294 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
277 34294 : if (gc_needed(av, 1))
278 0 : gerepileall(av, 2, &M2, &T3); /* destroy T2 */
279 34294 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
280 34294 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
281 34294 : M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
282 34294 : if (gc_needed(av, 1))
283 0 : gerepileall(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
284 34294 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
285 34294 : if (gc_needed(av, 1))
286 0 : gerepileall(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
287 34294 : A11 = matslice(A, 1, m1, 1, n1);
288 34294 : B11 = matslice(B, 1, n1, 1, p1);
289 34294 : M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
290 34294 : if (gc_needed(av, 1))
291 0 : gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
292 34294 : A12 = matslice(A, 1, m1, n1 + 1, n);
293 34294 : B21 = matslice(B, n1 + 1, n, 1, p1);
294 34294 : M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
295 34294 : if (gc_needed(av, 1))
296 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
297 34294 : C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
298 34294 : if (gc_needed(av, 1))
299 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11); /* destroy M4 */
300 34294 : M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
301 34294 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
302 34294 : if (gc_needed(av, 1))
303 5 : gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4); /* destroy S3 */
304 34294 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
305 34294 : if (gc_needed(av, 1))
306 0 : gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4); /* destroy T3 */
307 34294 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
308 34294 : if (gc_needed(av, 1))
309 1 : gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1); /* destroy M1, M5 */
310 34294 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
311 34294 : M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
312 34294 : if (gc_needed(av, 1))
313 6 : gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6); /* destroy S4, B22 */
314 34294 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
315 34294 : M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
316 34294 : if (gc_needed(av, 1))
317 0 : gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7); /* destroy A22, T4 */
318 34294 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
319 34294 : C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
320 34294 : if (gc_needed(av, 1))
321 12 : gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12); /* destroy V3, M6 */
322 34294 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
323 34294 : if (gc_needed(av, 1))
324 0 : gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2); /* destroy V1, M2 */
325 34294 : C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
326 34294 : if (gc_needed(av, 1))
327 6 : gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21); /* destroy M7 */
328 34294 : C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
329 34294 : if (gc_needed(av, 1))
330 0 : gerepileall(av, 4, &C11, &C12, &C21, &C22); /* destroy V2, M3 */
331 34294 : C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
332 34294 : return gerepilecopy(av, C);
333 : }
334 :
335 : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
336 : static GEN
337 444632540 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
338 : {
339 444632540 : pari_sp av = avma;
340 444632540 : GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
341 : long k;
342 4880370806 : for (k = 2; k < lx; k++)
343 : {
344 4436147661 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
345 4434047926 : if (t != ZERO) c = addii(c, t);
346 : }
347 444223145 : return gerepileuptoint(av, c);
348 : }
349 : GEN
350 150688514 : ZMrow_ZC_mul(GEN x, GEN y, long i)
351 150688514 : { 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 51528431 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
356 : {
357 51528431 : GEN z = cgetg(l,t_COL);
358 : long i;
359 345489362 : for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
360 51500439 : return z;
361 : }
362 :
363 : static GEN
364 8092590 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
365 : {
366 : long j;
367 8092590 : GEN z = cgetg(ly, t_MAT);
368 36692520 : for (j = 1; j < ly; j++)
369 28601355 : gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
370 8091165 : return z;
371 : }
372 :
373 : static GEN
374 1290 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
375 : {
376 1290 : pari_sp av = avma;
377 1290 : long i, n = lg(P)-1;
378 : GEN H, T;
379 1290 : 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 1290 : T = ZV_producttree(P);
388 1290 : A = ZM_nv_mod_tree(A, P, T);
389 1290 : B = ZM_nv_mod_tree(B, P, T);
390 1290 : H = cgetg(n+1, t_VEC);
391 9882 : for(i=1; i <= n; i++)
392 8592 : gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
393 1290 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
394 1290 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
395 : }
396 :
397 : GEN
398 1290 : ZM_mul_worker(GEN P, GEN A, GEN B)
399 : {
400 1290 : GEN V = cgetg(3, t_VEC);
401 1290 : gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
402 1290 : return V;
403 : }
404 :
405 : static GEN
406 919 : ZM_mul_fast(GEN A, GEN B, long lA, long lB, long sA, long sB)
407 : {
408 919 : pari_sp av = avma;
409 : forprime_t S;
410 : GEN H, worker;
411 : long h;
412 919 : if (sA == 2 || sB == 2) return zeromat(nbrows(A),lB-1);
413 905 : h = 1 + (sA + sB - 4) * BITS_IN_LONG + expu(lA-1);
414 905 : init_modular_big(&S);
415 905 : worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
416 905 : H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
417 : nmV_chinese_center, FpM_center);
418 905 : 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 8126854 : sw_bound(long s)
425 8126854 : { return s > 60 ? 2: s > 25 ? 4: s > 15 ? 8 : s > 8 ? 16 : 32; }
426 :
427 : static GEN
428 15307810 : ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
429 : {
430 : long sx, sy, B;
431 15307810 : if (lx==3 && l==3 && ly==3) return ZM2_mul(x,y);
432 8101932 : sx = ZM_max_lg_i(x, lx, l);
433 8102549 : sy = ZM_max_lg_i(y, ly, lx);
434 8102630 : if ((lx > 70 && ly > 70 && l > 70) && (sx <= 10 * sy && sy <= 10 * sx))
435 919 : return ZM_mul_fast(x, y, lx, ly, sx, sy);
436 :
437 8101711 : B = sw_bound(minss(sx, sy));
438 8101688 : if (l <= B || lx <= B || ly <= B)
439 8067424 : return ZM_mul_classical(x, y, l, lx, ly);
440 : else
441 34264 : return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
442 : }
443 :
444 : GEN
445 15203665 : ZM_mul(GEN x, GEN y)
446 : {
447 15203665 : long lx = lg(x), ly = lg(y);
448 15203665 : if (ly == 1) return cgetg(1,t_MAT);
449 15068985 : if (lx == 1) return zeromat(0, ly-1);
450 15067405 : return ZM_mul_i(x, y, lgcols(x), lx, ly);
451 : }
452 :
453 : GEN
454 538453 : QM_mul(GEN x, GEN y)
455 : {
456 538453 : GEN dx, nx = Q_primitive_part(x, &dx);
457 538453 : GEN dy, ny = Q_primitive_part(y, &dy);
458 538453 : GEN z = ZM_mul(nx, ny);
459 538453 : if (dx || dy)
460 : {
461 452820 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
462 452820 : if (!gequal1(d)) z = ZM_Q_mul(z, d);
463 : }
464 538453 : 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 132926 : QM_QC_mul(GEN x, GEN y)
479 : {
480 132926 : GEN dx, nx = Q_primitive_part(x, &dx);
481 132926 : GEN dy, ny = Q_primitive_part(y, &dy);
482 132926 : GEN z = ZM_ZC_mul(nx, ny);
483 132926 : if (dx || dy)
484 : {
485 132898 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
486 132898 : if (!gequal1(d)) z = ZC_Q_mul(z, d);
487 : }
488 132926 : 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 530461 : ZM_diag_mul(GEN d, GEN m)
529 : {
530 530461 : long i, j, l = lg(d), lm = lg(m);
531 530461 : GEN y = cgetg(lm, t_MAT);
532 1513696 : for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
533 1639876 : for (i=1; i<l; i++)
534 : {
535 1109744 : GEN c = gel(d,i);
536 1109744 : if (equali1(c))
537 236155 : for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
538 : else
539 3388730 : for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
540 : }
541 530132 : return y;
542 : }
543 :
544 : /* assume lx > 1 is lg(x) = lg(y) */
545 : static GEN
546 20888467 : ZV_dotproduct_i(GEN x, GEN y, long lx)
547 : {
548 20888467 : pari_sp av = avma;
549 20888467 : GEN c = mulii(gel(x,1), gel(y,1));
550 : long i;
551 147501445 : for (i = 2; i < lx; i++)
552 : {
553 126614779 : GEN t = mulii(gel(x,i), gel(y,i));
554 126614615 : if (t != gen_0) c = addii(c, t);
555 : }
556 20886666 : return gerepileuptoint(av, c);
557 : }
558 :
559 : /* x~ * y, assuming result is symmetric */
560 : GEN
561 526314 : ZM_transmultosym(GEN x, GEN y)
562 : {
563 526314 : long i, j, l, ly = lg(y);
564 : GEN M;
565 526314 : if (ly == 1) return cgetg(1,t_MAT);
566 : /* lg(x) = ly */
567 526314 : l = lgcols(y); /* = lgcols(x) */
568 526314 : M = cgetg(ly, t_MAT);
569 2678336 : for (i=1; i<ly; i++)
570 : {
571 2152022 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
572 2152022 : gel(M,i) = c;
573 7018586 : for (j=1; j<i; j++)
574 4866564 : gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
575 2152022 : gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
576 : }
577 526314 : return M;
578 : }
579 : /* x~ * y */
580 : GEN
581 2374 : ZM_transmul(GEN x, GEN y)
582 : {
583 2374 : long i, j, l, lx, ly = lg(y);
584 : GEN M;
585 2374 : if (ly == 1) return cgetg(1,t_MAT);
586 2374 : lx = lg(x);
587 2374 : l = lgcols(y);
588 2374 : if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
589 2374 : M = cgetg(ly, t_MAT);
590 7248 : for (i=1; i<ly; i++)
591 : {
592 4874 : GEN yi = gel(y,i), c = cgetg(lx,t_COL);
593 4874 : gel(M,i) = c;
594 12569 : for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
595 : }
596 2374 : 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 23014861 : ZM_ZC_mul(GEN x, GEN y)
619 : {
620 23014861 : long lx = lg(x);
621 23014861 : return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
622 : }
623 :
624 : GEN
625 1446932 : ZC_Z_div(GEN x, GEN c)
626 7120621 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
627 :
628 : GEN
629 14837 : ZM_Z_div(GEN x, GEN c)
630 168167 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
631 :
632 : GEN
633 2629652 : ZC_Q_mul(GEN A, GEN z)
634 : {
635 2629652 : pari_sp av = avma;
636 2629652 : long i, l = lg(A);
637 : GEN d, n, Ad, B, u;
638 2629652 : if (typ(z)==t_INT) return ZC_Z_mul(A,z);
639 2627139 : n = gel(z, 1); d = gel(z, 2);
640 2627139 : Ad = FpC_red(A, d);
641 2627082 : u = gcdii(d, FpV_factorback(Ad, NULL, d));
642 2627068 : B = cgetg(l, t_COL);
643 2627085 : if (equali1(u))
644 : {
645 295317 : for(i=1; i<l; i++)
646 236557 : gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
647 : } else
648 : {
649 17631077 : for(i=1; i<l; i++)
650 : {
651 15062844 : GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
652 15062723 : if (equalii(d, di))
653 10420598 : gel(B, i) = ni;
654 : else
655 4642147 : gel(B, i) = mkfrac(ni, diviiexact(d, di));
656 : }
657 : }
658 2626993 : return gerepilecopy(av, B);
659 : }
660 :
661 : GEN
662 1095530 : ZM_Q_mul(GEN x, GEN z)
663 : {
664 1095530 : if (typ(z)==t_INT) return ZM_Z_mul(x,z);
665 3162372 : 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 10119216 : ZV_dotsquare(GEN x)
708 : {
709 : long i, lx;
710 : pari_sp av;
711 : GEN z;
712 10119216 : lx = lg(x);
713 10119216 : if (lx == 1) return gen_0;
714 10119216 : av = avma; z = sqri(gel(x,1));
715 33334084 : for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
716 10116283 : return gerepileuptoint(av,z);
717 : }
718 :
719 : GEN
720 19740843 : ZV_dotproduct(GEN x,GEN y)
721 : {
722 : long lx;
723 19740843 : if (x == y) return ZV_dotsquare(x);
724 13230314 : lx = lg(x);
725 13230314 : if (lx == 1) return gen_0;
726 13230314 : 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 34237944 : ZC_add_i(GEN x, GEN y, long lx)
756 : {
757 34237944 : GEN A = cgetg(lx, t_COL);
758 : long i;
759 472043560 : for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
760 34224582 : return A;
761 : }
762 : GEN
763 24926027 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
764 : GEN
765 361337 : ZC_Z_add(GEN x, GEN y)
766 : {
767 361337 : long k, lx = lg(x);
768 361337 : GEN z = cgetg(lx, t_COL);
769 361336 : if (lx == 1) pari_err_TYPE2("+",x,y);
770 361336 : gel(z,1) = addii(y,gel(x,1));
771 2479872 : for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
772 361355 : return z;
773 : }
774 :
775 : static GEN
776 8933827 : ZC_sub_i(GEN x, GEN y, long lx)
777 : {
778 : long i;
779 8933827 : GEN A = cgetg(lx, t_COL);
780 55977652 : for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
781 8932629 : return A;
782 : }
783 : GEN
784 8642823 : 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 444912 : Z_ZC_sub(GEN a, GEN x)
797 : {
798 444912 : long k, lx = lg(x);
799 444912 : GEN z = cgetg(lx, t_COL);
800 444929 : if (lx == 1) pari_err_TYPE2("-",a,x);
801 444929 : gel(z,1) = subii(a, gel(x,1));
802 1212737 : for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
803 444952 : return z;
804 : }
805 :
806 : GEN
807 818087 : ZM_add(GEN x, GEN y)
808 : {
809 818087 : long lx = lg(x), l, j;
810 : GEN z;
811 818087 : if (lx == 1) return cgetg(1, t_MAT);
812 738819 : z = cgetg(lx, t_MAT); l = lgcols(x);
813 10049880 : for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
814 738817 : return z;
815 : }
816 : GEN
817 113046 : ZM_sub(GEN x, GEN y)
818 : {
819 113046 : long lx = lg(x), l, j;
820 : GEN z;
821 113046 : if (lx == 1) return cgetg(1, t_MAT);
822 113046 : z = cgetg(lx, t_MAT); l = lgcols(x);
823 404020 : for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
824 113045 : return z;
825 : }
826 : /********************************************************************/
827 : /** **/
828 : /** LINEAR COMBINATION **/
829 : /** **/
830 : /********************************************************************/
831 : /* return X/c assuming division is exact */
832 : GEN
833 9212682 : ZC_Z_divexact(GEN x, GEN c)
834 77907399 : { 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 2067930 : ZM_Z_divexact(GEN x, GEN c)
841 9223537 : { 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 33154284 : ZC_Z_mul(GEN x, GEN c)
849 : {
850 33154284 : if (!signe(c)) return zerocol(lg(x)-1);
851 31909875 : if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
852 249913564 : pari_APPLY_type(t_COL, mulii(gel(x,i), c))
853 : }
854 :
855 : GEN
856 47563 : ZC_z_mul(GEN x, long c)
857 : {
858 47563 : if (!c) return zerocol(lg(x)-1);
859 40013 : if (c == 1) return ZC_copy(x);
860 35453 : if (c ==-1) return ZC_neg(x);
861 441132 : pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
862 : }
863 :
864 : GEN
865 28529 : zv_z_mul(GEN x, long n)
866 431552 : { 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 3248172 : ZM_Z_mul(GEN X, GEN c)
892 : {
893 3248172 : long i, j, h, l = lg(X);
894 : GEN A;
895 3248172 : if (l == 1) return cgetg(1, t_MAT);
896 3248172 : h = lgcols(X);
897 3248168 : if (!signe(c)) return zeromat(h-1, l-1);
898 3203103 : if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
899 2919517 : A = cgetg(l, t_MAT);
900 18084035 : for (j = 1; j < l; j++)
901 : {
902 15165754 : GEN a = cgetg(h, t_COL), x = gel(X, j);
903 229517801 : for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
904 15164512 : gel(A,j) = a;
905 : }
906 2918281 : return A;
907 : }
908 : void
909 90017350 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
910 : {
911 : long i;
912 2038929935 : for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
913 89220946 : }
914 : /* X <- X + v Y (elementary col operation) */
915 : void
916 79542626 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
917 : {
918 79542626 : if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
919 : }
920 : void
921 29271253 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
922 : {
923 : long i;
924 29271253 : if (!v) return; /* v = 0 */
925 745078356 : 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 17203972 : ZC_lincomb1(GEN v, GEN x, GEN y)
931 127234601 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
932 :
933 : /* -X + vY */
934 : static GEN
935 603214 : ZC_lincomb_1(GEN v, GEN x, GEN y)
936 4079253 : { 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 33904231 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
941 : {
942 : long su, sv;
943 : GEN A;
944 :
945 33904231 : su = signe(u); if (!su) return ZC_Z_mul(Y, v);
946 33903273 : sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
947 33476201 : if (is_pm1(v))
948 : {
949 10547810 : if (is_pm1(u))
950 : {
951 9379381 : if (su != sv) A = ZC_sub(X, Y);
952 2479334 : else A = ZC_add(X, Y);
953 9378794 : if (su < 0) ZV_togglesign(A); /* in place but was created above */
954 : }
955 : else
956 : {
957 1168400 : if (sv > 0) A = ZC_lincomb1 (u, Y, X);
958 471698 : else A = ZC_lincomb_1(u, Y, X);
959 : }
960 : }
961 22929702 : else if (is_pm1(u))
962 : {
963 16638338 : if (su > 0) A = ZC_lincomb1 (v, X, Y);
964 130993 : else A = ZC_lincomb_1(v, X, Y);
965 : }
966 : else
967 : { /* not cgetg_copy: x may be a t_VEC */
968 6294943 : long i, lx = lg(X);
969 6294943 : A = cgetg(lx,t_COL);
970 40420522 : for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
971 : }
972 33469010 : return A;
973 : }
974 :
975 : /********************************************************************/
976 : /** **/
977 : /** CONVERSIONS **/
978 : /** **/
979 : /********************************************************************/
980 : GEN
981 393789 : ZV_to_nv(GEN x)
982 739692 : { pari_APPLY_ulong(itou(gel(x,i))) }
983 :
984 : GEN
985 133485 : zm_to_ZM(GEN x)
986 744964 : { 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 12615181 : ZC_copy(GEN x)
1016 : {
1017 12615181 : long i, lx = lg(x);
1018 12615181 : GEN y = cgetg(lx, t_COL);
1019 81807744 : for (i=1; i<lx; i++)
1020 : {
1021 69192802 : GEN c = gel(x,i);
1022 69192802 : gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
1023 : }
1024 12614942 : return y;
1025 : }
1026 :
1027 : GEN
1028 512716 : ZM_copy(GEN x)
1029 4308892 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
1030 :
1031 : void
1032 340178 : ZV_neg_inplace(GEN M)
1033 : {
1034 340178 : long l = lg(M);
1035 1252026 : while (--l > 0) gel(M,l) = negi(gel(M,l));
1036 340200 : }
1037 : GEN
1038 3523854 : ZC_neg(GEN x)
1039 27107321 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
1040 :
1041 : GEN
1042 51336 : zv_neg(GEN x)
1043 658452 : { pari_APPLY_long(-x[i]) }
1044 : GEN
1045 126 : zv_neg_inplace(GEN M)
1046 : {
1047 126 : long l = lg(M);
1048 534 : 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 319744 : ZM_neg(GEN x)
1056 1310562 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
1057 :
1058 : void
1059 4988902 : ZV_togglesign(GEN M)
1060 : {
1061 4988902 : long l = lg(M);
1062 74560597 : while (--l > 0) togglesign_safe(&gel(M,l));
1063 4989025 : }
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 11020721 : ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
1079 : {
1080 11020721 : long i, l = lg(x);
1081 : GEN q;
1082 :
1083 11020721 : if (Q) *Q = cgetg(l,t_COL);
1084 11020721 : if (l == 1) return cgetg(1,t_COL);
1085 58465391 : for (i = l-1; i>0; i--)
1086 : {
1087 47447588 : q = diviiround(gel(x,i), gcoeff(y,i,i));
1088 47446117 : if (signe(q)) {
1089 24488096 : togglesign(q);
1090 24488329 : x = ZC_lincomb(gen_1, q, x, gel(y,i));
1091 : }
1092 47444670 : if (Q) gel(*Q, i) = q;
1093 : }
1094 11017803 : return x;
1095 : }
1096 :
1097 : /* x = y Q + R, may return some columns of x (not copies) */
1098 : GEN
1099 453542 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
1100 : {
1101 453542 : long lx = lg(x), i;
1102 453542 : GEN R = cgetg(lx, t_MAT);
1103 453550 : if (Q)
1104 : {
1105 127200 : GEN q = cgetg(lx, t_MAT); *Q = q;
1106 187983 : for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
1107 : }
1108 : else
1109 822652 : for (i=1; i<lx; i++)
1110 : {
1111 496299 : pari_sp av = avma;
1112 496299 : GEN z = ZC_hnfrem(gel(x,i),y);
1113 496283 : gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
1114 : }
1115 453555 : return R;
1116 : }
1117 :
1118 : /********************************************************************/
1119 : /** **/
1120 : /** TESTS **/
1121 : /** **/
1122 : /********************************************************************/
1123 : int
1124 23099474 : zv_equal0(GEN V)
1125 : {
1126 23099474 : long l = lg(V);
1127 37561040 : while (--l > 0)
1128 30787797 : if (V[l]) return 0;
1129 6773243 : return 1;
1130 : }
1131 :
1132 : int
1133 1693208 : ZV_equal0(GEN V)
1134 : {
1135 1693208 : long l = lg(V);
1136 4597394 : while (--l > 0)
1137 4394224 : if (signe(gel(V,l))) return 0;
1138 203170 : return 1;
1139 : }
1140 : int
1141 15974 : ZMrow_equal0(GEN V, long i)
1142 : {
1143 15974 : long l = lg(V);
1144 24891 : while (--l > 0)
1145 21429 : if (signe(gcoeff(V,i,l))) return 0;
1146 3462 : return 1;
1147 : }
1148 :
1149 : static int
1150 4462466 : ZV_equal_lg(GEN V, GEN W, long l)
1151 : {
1152 19239418 : while (--l > 0)
1153 15179813 : if (!equalii(gel(V,l), gel(W,l))) return 0;
1154 4059605 : return 1;
1155 : }
1156 : int
1157 238358 : ZV_equal(GEN V, GEN W)
1158 : {
1159 238358 : long l = lg(V);
1160 238358 : if (lg(W) != l) return 0;
1161 238358 : return ZV_equal_lg(V, W, l);
1162 : }
1163 : int
1164 1767985 : ZM_equal(GEN A, GEN B)
1165 : {
1166 1767985 : long i, m, l = lg(A);
1167 1767985 : if (lg(B) != l) return 0;
1168 1767627 : if (l == 1) return 1;
1169 1767627 : m = lgcols(A);
1170 1767626 : if (lgcols(B) != m) return 0;
1171 5766108 : for (i = 1; i < l; i++)
1172 4224109 : if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
1173 1541999 : return 1;
1174 : }
1175 : int
1176 70217 : ZM_equal0(GEN A)
1177 : {
1178 70217 : long i, j, m, l = lg(A);
1179 70217 : if (l == 1) return 1;
1180 70217 : m = lgcols(A);
1181 121019 : for (j = 1; j < l; j++)
1182 2616058 : for (i = 1; i < m; i++)
1183 2565256 : if (signe(gcoeff(A,i,j))) return 0;
1184 15541 : return 1;
1185 : }
1186 : int
1187 30764946 : zv_equal(GEN V, GEN W)
1188 : {
1189 30764946 : long l = lg(V);
1190 30764946 : if (lg(W) != l) return 0;
1191 261905846 : while (--l > 0)
1192 232214731 : if (V[l] != W[l]) return 0;
1193 29691115 : 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 748494 : ZM_ishnf(GEN x)
1208 : {
1209 748494 : long i,j, lx = lg(x);
1210 2266804 : for (i=1; i<lx; i++)
1211 : {
1212 1621538 : GEN xii = gcoeff(x,i,i);
1213 1621538 : if (signe(xii) <= 0) return 0;
1214 3165460 : for (j=1; j<i; j++)
1215 1574162 : if (signe(gcoeff(x,i,j))) return 0;
1216 3221145 : for (j=i+1; j<lx; j++)
1217 : {
1218 1702835 : GEN xij = gcoeff(x,i,j);
1219 1702835 : if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
1220 : }
1221 : }
1222 645266 : return 1;
1223 : }
1224 : int
1225 477833 : ZM_isidentity(GEN x)
1226 : {
1227 477833 : long i,j, lx = lg(x);
1228 :
1229 477833 : if (lx == 1) return 1;
1230 477826 : if (lx != lgcols(x)) return 0;
1231 2348533 : for (j=1; j<lx; j++)
1232 : {
1233 1870763 : GEN c = gel(x,j);
1234 5924740 : for (i=1; i<j; )
1235 4053998 : if (signe(gel(c,i++))) return 0;
1236 : /* i = j */
1237 1870742 : if (!equali1(gel(c,i++))) return 0;
1238 5924754 : for ( ; i<lx; )
1239 4054033 : if (signe(gel(c,i++))) return 0;
1240 : }
1241 477770 : return 1;
1242 : }
1243 : int
1244 555397 : ZM_isdiagonal(GEN x)
1245 : {
1246 555397 : long i,j, lx = lg(x);
1247 555397 : if (lx == 1) return 1;
1248 555397 : if (lx != lgcols(x)) return 0;
1249 :
1250 1434074 : for (j=1; j<lx; j++)
1251 : {
1252 1205032 : GEN c = gel(x,j);
1253 1687700 : for (i=1; i<j; i++)
1254 809006 : if (signe(gel(c,i))) return 0;
1255 2007392 : for (i++; i<lx; i++)
1256 1128733 : if (signe(gel(c,i))) return 0;
1257 : }
1258 229042 : return 1;
1259 : }
1260 : int
1261 106678 : ZM_isscalar(GEN x, GEN s)
1262 : {
1263 106678 : long i, j, lx = lg(x);
1264 :
1265 106678 : if (lx == 1) return 1;
1266 106678 : if (!s) s = gcoeff(x,1,1);
1267 106678 : if (equali1(s)) return ZM_isidentity(x);
1268 105425 : if (lx != lgcols(x)) return 0;
1269 555866 : for (j=1; j<lx; j++)
1270 : {
1271 467706 : GEN c = gel(x,j);
1272 2384584 : for (i=1; i<j; )
1273 1931924 : if (signe(gel(c,i++))) return 0;
1274 : /* i = j */
1275 452660 : if (!equalii(gel(c,i++), s)) return 0;
1276 2428013 : for ( ; i<lx; )
1277 1977572 : if (signe(gel(c,i++))) return 0;
1278 : }
1279 88160 : 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 2963245 : ZV_cmp(GEN x, GEN y)
1305 : {
1306 2963245 : long fl,i, lx = lg(x);
1307 5694906 : for (i=1; i<lx; i++)
1308 4575009 : if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
1309 1119897 : 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 68161 : for (i=1; i<lx; i++)
1317 68039 : if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
1318 122 : return 0;
1319 : }
1320 :
1321 : long
1322 23092930 : zv_content(GEN x)
1323 : {
1324 23092930 : long i, s, l = lg(x);
1325 23092930 : if (l == 1) return 0;
1326 23092923 : s = labs(x[1]);
1327 51268321 : for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
1328 23093814 : return s;
1329 : }
1330 : GEN
1331 406299 : ZV_content(GEN x)
1332 : {
1333 406299 : long i, l = lg(x);
1334 406299 : pari_sp av = avma;
1335 : GEN c;
1336 406299 : if (l == 1) return gen_0;
1337 406299 : if (l == 2) return absi(gel(x,1));
1338 332057 : c = gel(x,1);
1339 895281 : for (i = 2; i < l; i++)
1340 : {
1341 665957 : c = gcdii(c, gel(x,i));
1342 665952 : if (is_pm1(c)) { set_avma(av); return gen_1; }
1343 : }
1344 229324 : return gerepileuptoint(av, c);
1345 : }
1346 :
1347 : GEN
1348 3861950 : ZM_det_triangular(GEN mat)
1349 : {
1350 : pari_sp av;
1351 3861950 : long i,l = lg(mat);
1352 : GEN s;
1353 :
1354 3861950 : if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
1355 3462076 : av = avma; s = gcoeff(mat,1,1);
1356 9389123 : for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
1357 3461614 : return gerepileuptoint(av,s);
1358 : }
1359 :
1360 : /* assumes no overflow */
1361 : long
1362 939622 : zv_prod(GEN v)
1363 : {
1364 939622 : long n, i, l = lg(v);
1365 939622 : if (l == 1) return 1;
1366 951796 : n = v[1]; for (i = 2; i < l; i++) n *= v[i];
1367 763593 : return n;
1368 : }
1369 :
1370 : static GEN
1371 355912400 : _mulii(void *E, GEN a, GEN b)
1372 355912400 : { (void) E; return mulii(a, b); }
1373 :
1374 : /* product of ulongs */
1375 : GEN
1376 1861626 : zv_prod_Z(GEN v)
1377 : {
1378 1861626 : pari_sp av = avma;
1379 1861626 : long k, m, n = lg(v)-1;
1380 : GEN V;
1381 1861626 : switch(n) {
1382 20902 : case 0: return gen_1;
1383 121058 : case 1: return utoi(v[1]);
1384 978419 : case 2: return muluu(v[1], v[2]);
1385 : }
1386 741247 : m = n >> 1;
1387 741247 : V = cgetg(m + (odd(n)? 2: 1), t_VEC);
1388 84985339 : for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
1389 741180 : if (odd(n)) gel(V,k) = utoipos(v[n]);
1390 741241 : 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 8811185 : ZV_prod(GEN v)
1412 : {
1413 8811185 : pari_sp av = avma;
1414 8811185 : long i, l = lg(v);
1415 : GEN n;
1416 8811185 : if (l == 1) return gen_1;
1417 8629832 : if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
1418 1287728 : n = gel(v,1);
1419 1287728 : if (l == 2) return icopy(n);
1420 2088155 : for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
1421 839214 : return gerepileuptoint(av, n);
1422 : }
1423 : /* assumes no overflow */
1424 : long
1425 16464 : zv_sum(GEN v)
1426 : {
1427 16464 : long n, i, l = lg(v);
1428 16464 : if (l == 1) return 0;
1429 95924 : n = v[1]; for (i = 2; i < l; i++) n += v[i];
1430 16443 : 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 360456 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
1463 : {
1464 360456 : long i, qq = itos_or_0(q);
1465 360457 : if (!qq)
1466 : {
1467 27955 : for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
1468 6453 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
1469 6453 : return;
1470 : }
1471 354004 : if (qq == 1) {
1472 149995 : for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
1473 101431 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
1474 252573 : } else if (qq == -1) {
1475 154396 : for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
1476 89806 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
1477 : } else {
1478 273721 : for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
1479 162784 : gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
1480 : }
1481 : }
1482 :
1483 : /* update L[k,] */
1484 : static void
1485 1032821 : ZRED(long k, long l, GEN x, GEN L, GEN B)
1486 : {
1487 1032821 : GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
1488 1032860 : if (!signe(q)) return;
1489 360429 : q = negi(q);
1490 360452 : Zupdate_row(k,l,q,L,B);
1491 360443 : 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 1182910 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
1497 : {
1498 : long i, j;
1499 3777913 : for (j=1; j<=k; j++)
1500 : {
1501 2595610 : pari_sp av = avma;
1502 2595610 : GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
1503 5603748 : for (i=1; i<j; i++)
1504 : {
1505 3009416 : u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
1506 3008853 : u = diviiexact(u, gel(B,i));
1507 : }
1508 2594332 : gcoeff(L,k,j) = gerepileuptoint(av, u);
1509 : }
1510 1182303 : gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
1511 1182303 : }
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 110114 : ZC_reducemodmatrix_i(GEN v, GEN y)
1517 : {
1518 110114 : GEN B, L, x = shallowconcat(y, v);
1519 110114 : long k, lx = lg(x), nx = lx-1;
1520 :
1521 110114 : B = scalarcol_shallow(gen_1, lx);
1522 110115 : L = zeromatcopy(nx, nx);
1523 453928 : for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
1524 343812 : for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1525 110101 : return gel(x,nx);
1526 : }
1527 : GEN
1528 110114 : ZC_reducemodmatrix(GEN v, GEN y) {
1529 110114 : pari_sp av = avma;
1530 110114 : 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 226674 : ZM_reducemodmatrix_i(GEN v, GEN y)
1537 : {
1538 : GEN B, L, V;
1539 226674 : long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
1540 :
1541 226674 : V = cgetg(lv, t_MAT);
1542 226676 : B = scalarcol_shallow(gen_1, lx);
1543 226681 : L = zeromatcopy(nx, nx);
1544 600467 : for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
1545 692001 : for (j = 1; j < lg(v); j++)
1546 : {
1547 465347 : GEN x = shallowconcat(y, gel(v,j));
1548 465370 : ZincrementalGS(x, L, B, nx); /* overwrite last */
1549 1264494 : for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
1550 465332 : gel(V,j) = gel(x,nx);
1551 : }
1552 226654 : return V;
1553 : }
1554 : GEN
1555 226671 : ZM_reducemodmatrix(GEN v, GEN y) {
1556 226671 : pari_sp av = avma;
1557 226671 : return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
1558 : }
1559 :
1560 : GEN
1561 98861 : ZC_reducemodlll(GEN x,GEN y)
1562 : {
1563 98861 : pari_sp av = avma;
1564 98861 : GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
1565 98862 : 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 : }
|
