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 : int
19 1110488 : RgM_is_ZM(GEN x)
20 : {
21 1110488 : long i, j, h, l = lg(x);
22 1110488 : if (l == 1) return 1;
23 1109914 : h = lgcols(x);
24 1109914 : if (h == 1) return 1;
25 4193679 : for (j = l-1; j > 0; j--)
26 19718062 : for (i = h-1; i > 0; i--)
27 16634159 : if (typ(gcoeff(x,i,j)) != t_INT) return 0;
28 969602 : return 1;
29 : }
30 :
31 : int
32 174 : RgM_is_QM(GEN x)
33 : {
34 174 : long i, j, h, l = lg(x);
35 174 : if (l == 1) return 1;
36 162 : h = lgcols(x);
37 162 : if (h == 1) return 1;
38 1044 : for (j = l-1; j > 0; j--)
39 12906 : for (i = h-1; i > 0; i--)
40 12024 : if (!is_rational_t(typ(gcoeff(x,i,j)))) return 0;
41 150 : return 1;
42 : }
43 :
44 : int
45 35 : RgV_is_ZMV(GEN V)
46 : {
47 35 : long i, l = lg(V);
48 285 : for (i=1; i<l; i++)
49 250 : if (typ(gel(V,i))!=t_MAT || !RgM_is_ZM(gel(V,i)))
50 0 : return 0;
51 35 : return 1;
52 : }
53 :
54 : /********************************************************************/
55 : /** **/
56 : /** GENERIC LINEAR ALGEBRA **/
57 : /** **/
58 : /********************************************************************/
59 : /* GENERIC MULTIPLICATION involving zc/zm */
60 :
61 : /* x[i,] * y */
62 : static GEN
63 7043170 : RgMrow_zc_mul_i(GEN x, GEN y, long c, long i)
64 : {
65 7043170 : pari_sp av = avma;
66 7043170 : GEN s = NULL;
67 : long j;
68 696400541 : for (j=1; j<c; j++)
69 : {
70 689357371 : long t = y[j];
71 689357371 : if (!t) continue;
72 24013723 : if (!s) { s = gmulgs(gcoeff(x,i,j),t); continue; }
73 17155761 : switch(t)
74 : {
75 12459573 : case 1: s = gadd(s, gcoeff(x,i,j)); break;
76 202738 : case -1: s = gsub(s, gcoeff(x,i,j)); break;
77 4493450 : default: s = gadd(s, gmulgs(gcoeff(x,i,j), t)); break;
78 : }
79 : }
80 7043170 : return s? gc_upto(av, s): gc_const(av, gen_0);
81 : }
82 : GEN
83 60024 : RgMrow_zc_mul(GEN x, GEN y, long i) { return RgMrow_zc_mul_i(x,y,lg(y),i); }
84 : /* x nonempty t_MAT, y a compatible zc (dimension > 0). */
85 : static GEN
86 197638 : RgM_zc_mul_i(GEN x, GEN y, long c, long l)
87 : {
88 197638 : GEN z = cgetg(l,t_COL);
89 : long i;
90 7180784 : for (i = 1; i < l; i++) gel(z,i) = RgMrow_zc_mul_i(x,y,c,i);
91 197638 : return z;
92 : }
93 : GEN
94 51632 : RgM_zc_mul(GEN x, GEN y) { return RgM_zc_mul_i(x,y, lg(x), lgcols(x)); }
95 : /* x t_MAT, y a compatible zm (dimension > 0). */
96 : GEN
97 49688 : RgM_zm_mul(GEN x, GEN y)
98 : {
99 49688 : long j, c, l = lg(x), ly = lg(y);
100 49688 : GEN z = cgetg(ly, t_MAT);
101 49688 : if (l == 1) return z;
102 49688 : c = lgcols(x);
103 195694 : for (j = 1; j < ly; j++) gel(z,j) = RgM_zc_mul_i(x, gel(y,j), l,c);
104 49688 : return z;
105 : }
106 :
107 : /* x[i,]*y, l = lg(y) > 1 */
108 : static GEN
109 27852196 : RgMrow_ZC_mul_i(GEN x, GEN y, long i, long l)
110 : {
111 27852196 : pari_sp av = avma;
112 27852196 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
113 : long j;
114 4877023098 : for (j=2; j<l; j++)
115 4849170902 : if (signe(gel(y,j))) t = gadd(t, gmul(gcoeff(x,i,j), gel(y,j)));
116 27852196 : return gc_upto(av,t);
117 : }
118 :
119 : /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
120 : static GEN
121 1654678 : RgM_ZC_mul_i(GEN x, GEN y, long lx, long l)
122 : {
123 1654678 : GEN z = cgetg(l,t_COL);
124 : long i;
125 29506874 : for (i=1; i<l; i++) gel(z,i) = RgMrow_ZC_mul_i(x,y,i,lx);
126 1654678 : return z;
127 : }
128 :
129 : GEN
130 1654678 : RgM_ZM_mul_worker(GEN y, GEN x)
131 1654678 : { return RgM_ZC_mul_i(x, y, lg(x), lgcols(x)); }
132 :
133 : /* mostly useful when y is sparse */
134 : GEN
135 304230 : RgM_ZM_mul(GEN x, GEN y)
136 : {
137 304230 : pari_sp av = avma;
138 : GEN worker;
139 304230 : if (lg(x) == 1) return cgetg(lg(y), t_MAT);
140 304230 : worker = snm_closure(is_entry("_RgM_ZM_mul_worker"),mkvec(x));
141 304230 : return gc_upto(av, gen_parapply(worker, y));
142 : }
143 :
144 : static GEN
145 161766 : RgV_zc_mul_i(GEN x, GEN y, long l)
146 : {
147 : long i;
148 161766 : GEN z = gen_0;
149 161766 : pari_sp av = avma;
150 6489021 : for (i = 1; i < l; i++) z = gadd(z, gmulgs(gel(x,i), y[i]));
151 161766 : return gc_upto(av, z);
152 : }
153 : GEN
154 71532 : RgV_zc_mul(GEN x, GEN y) { return RgV_zc_mul_i(x, y, lg(x)); }
155 :
156 : GEN
157 24389 : RgV_zm_mul(GEN x, GEN y)
158 : {
159 24389 : long j, l = lg(x), ly = lg(y);
160 24389 : GEN z = cgetg(ly, t_VEC);
161 114623 : for (j = 1; j < ly; j++) gel(z,j) = RgV_zc_mul_i(x, gel(y,j), l);
162 24389 : return z;
163 : }
164 :
165 : /* scalar product x.x */
166 : GEN
167 8141684 : RgV_dotsquare(GEN x)
168 : {
169 8141684 : long i, lx = lg(x);
170 8141684 : pari_sp av = avma;
171 : GEN z;
172 8141684 : if (lx == 1) return gen_0;
173 8141684 : z = gsqr(gel(x,1));
174 18127908 : for (i=2; i<lx; i++)
175 : {
176 9986224 : z = gadd(z, gsqr(gel(x,i)));
177 9986224 : if (gc_needed(av,3))
178 : {
179 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotsquare, i = %ld",i);
180 0 : z = gc_upto(av, z);
181 : }
182 : }
183 8141684 : return gc_upto(av,z);
184 : }
185 :
186 : /* scalar product x.y, lx = lg(x) = lg(y) */
187 : static GEN
188 11540814 : RgV_dotproduct_i(GEN x, GEN y, long lx)
189 : {
190 11540814 : pari_sp av = avma;
191 : long i;
192 : GEN z;
193 11540814 : if (lx == 1) return gen_0;
194 11540472 : z = gmul(gel(x,1),gel(y,1));
195 177157636 : for (i=2; i<lx; i++)
196 : {
197 165617164 : z = gadd(z, gmul(gel(x,i), gel(y,i)));
198 165617164 : if (gc_needed(av,3))
199 : {
200 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotproduct, i = %ld",i);
201 0 : z = gc_upto(av, z);
202 : }
203 : }
204 11540472 : return gc_upto(av,z);
205 : }
206 : GEN
207 869364 : RgV_dotproduct(GEN x,GEN y)
208 : {
209 869364 : if (x == y) return RgV_dotsquare(x);
210 869364 : return RgV_dotproduct_i(x, y, lg(x));
211 : }
212 : /* v[1] + ... + v[lg(v)-1] */
213 : GEN
214 2682587 : RgV_sum(GEN v)
215 : {
216 : GEN p;
217 2682587 : long i, l = lg(v);
218 2682587 : if (l == 1) return gen_0;
219 6697532 : p = gel(v,1); for (i=2; i<l; i++) p = gadd(p, gel(v,i));
220 2682587 : return p;
221 : }
222 : /* v[1] + ... + v[n]. Assume lg(v) > n. */
223 : GEN
224 0 : RgV_sumpart(GEN v, long n)
225 : {
226 : GEN p;
227 : long i;
228 0 : if (!n) return gen_0;
229 0 : p = gel(v,1); for (i=2; i<=n; i++) p = gadd(p, gel(v,i));
230 0 : return p;
231 : }
232 : /* v[m] + ... + v[n]. Assume lg(v) > n, m > 0. */
233 : GEN
234 0 : RgV_sumpart2(GEN v, long m, long n)
235 : {
236 : GEN p;
237 : long i;
238 0 : if (n < m) return gen_0;
239 0 : p = gel(v,m); for (i=m+1; i<=n; i++) p = gadd(p, gel(v,i));
240 0 : return p;
241 : }
242 : GEN
243 288 : RgM_sumcol(GEN A)
244 : {
245 288 : long i,j,m,l = lg(A);
246 : GEN v;
247 :
248 288 : if (l == 1) return cgetg(1,t_MAT);
249 288 : if (l == 2) return gcopy(gel(A,1));
250 156 : m = lgcols(A);
251 156 : v = cgetg(m, t_COL);
252 518 : for (i = 1; i < m; i++)
253 : {
254 362 : pari_sp av = avma;
255 362 : GEN s = gcoeff(A,i,1);
256 874 : for (j = 2; j < l; j++) s = gadd(s, gcoeff(A,i,j));
257 362 : gel(v, i) = gc_upto(av, s);
258 : }
259 156 : return v;
260 : }
261 :
262 : static GEN
263 2819855 : _gmul(void *data, GEN x, GEN y)
264 2819855 : { (void)data; return gmul(x,y); }
265 :
266 : GEN
267 386537 : RgV_prod(GEN x)
268 : {
269 386537 : return gen_product(x, NULL, _gmul);
270 : }
271 :
272 : /* ADDITION SCALAR + MATRIX */
273 : /* x square matrix, y scalar; create the square matrix x + y*Id */
274 : GEN
275 6533 : RgM_Rg_add(GEN x, GEN y)
276 : {
277 6533 : long l = lg(x), i, j;
278 6533 : GEN z = cgetg(l,t_MAT);
279 :
280 6533 : if (l==1) return z;
281 6533 : if (l != lgcols(x)) pari_err_OP("+", x, y);
282 49177 : for (i=1; i<l; i++)
283 : {
284 42644 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
285 42644 : gel(z,i) = zi;
286 1763496 : for (j=1; j<l; j++)
287 1720852 : gel(zi,j) = i==j? gadd(y,gel(xi,j)): gcopy(gel(xi,j));
288 : }
289 6533 : return z;
290 : }
291 : GEN
292 0 : RgM_Rg_sub(GEN x, GEN y)
293 : {
294 0 : long l = lg(x), i, j;
295 0 : GEN z = cgetg(l,t_MAT);
296 :
297 0 : if (l==1) return z;
298 0 : if (l != lgcols(x)) pari_err_OP("-", x, y);
299 0 : for (i=1; i<l; i++)
300 : {
301 0 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
302 0 : gel(z,i) = zi;
303 0 : for (j=1; j<l; j++)
304 0 : gel(zi,j) = i==j? gsub(gel(xi,j), y): gcopy(gel(xi,j));
305 : }
306 0 : return z;
307 : }
308 : GEN
309 604 : RgM_Rg_add_shallow(GEN x, GEN y)
310 : {
311 604 : long l = lg(x), i, j;
312 604 : GEN z = cgetg(l,t_MAT);
313 :
314 604 : if (l==1) return z;
315 604 : if (l != lgcols(x)) pari_err_OP( "+", x, y);
316 1965 : for (i=1; i<l; i++)
317 : {
318 1361 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
319 1361 : gel(z,i) = zi;
320 7958 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
321 1361 : gel(zi,i) = gadd(gel(zi,i), y);
322 : }
323 604 : return z;
324 : }
325 : GEN
326 269864 : RgM_Rg_sub_shallow(GEN x, GEN y)
327 : {
328 269864 : long l = lg(x), i, j;
329 269864 : GEN z = cgetg(l,t_MAT);
330 :
331 269864 : if (l==1) return z;
332 269864 : if (l != lgcols(x)) pari_err_OP( "-", x, y);
333 1728331 : for (i=1; i<l; i++)
334 : {
335 1458467 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
336 1458467 : gel(z,i) = zi;
337 23098976 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
338 1458467 : gel(zi,i) = gsub(gel(zi,i), y);
339 : }
340 269864 : return z;
341 : }
342 :
343 : GEN
344 2059425 : RgC_Rg_add(GEN x, GEN y)
345 : {
346 2059425 : long k, lx = lg(x);
347 2059425 : GEN z = cgetg(lx, t_COL);
348 2059425 : if (lx == 1)
349 : {
350 0 : if (isintzero(y)) return z;
351 0 : pari_err_TYPE2("+",x,y);
352 : }
353 2059425 : gel(z,1) = gadd(y,gel(x,1));
354 5794403 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
355 2059425 : return z;
356 : }
357 : GEN
358 19808 : RgC_Rg_sub(GEN x, GEN y)
359 : {
360 19808 : long k, lx = lg(x);
361 19808 : GEN z = cgetg(lx, t_COL);
362 19808 : if (lx == 1)
363 : {
364 0 : if (isintzero(y)) return z;
365 0 : pari_err_TYPE2("-",x,y);
366 : }
367 19808 : gel(z,1) = gsub(gel(x,1), y);
368 56774 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
369 19808 : return z;
370 : }
371 : /* a - x */
372 : GEN
373 26734 : Rg_RgC_sub(GEN a, GEN x)
374 : {
375 26734 : long k, lx = lg(x);
376 26734 : GEN z = cgetg(lx,t_COL);
377 26734 : if (lx == 1)
378 : {
379 0 : if (isintzero(a)) return z;
380 0 : pari_err_TYPE2("-",a,x);
381 : }
382 26734 : gel(z,1) = gsub(a, gel(x,1));
383 60640 : for (k = 2; k < lx; k++) gel(z,k) = gneg(gel(x,k));
384 26734 : return z;
385 : }
386 :
387 : static GEN
388 16120481 : RgC_add_i(GEN x, GEN y, long lx)
389 : {
390 16120481 : GEN A = cgetg(lx, t_COL);
391 : long i;
392 122321470 : for (i=1; i<lx; i++) gel(A,i) = gadd(gel(x,i), gel(y,i));
393 16120481 : return A;
394 : }
395 : GEN
396 13322979 : RgC_add(GEN x, GEN y) { return RgC_add_i(x, y, lg(x)); }
397 : GEN
398 2210696 : RgV_add(GEN x, GEN y)
399 8906703 : { pari_APPLY_type(t_VEC, gadd(gel(x,i), gel(y,i))) }
400 :
401 : static GEN
402 4962179 : RgC_sub_i(GEN x, GEN y, long lx)
403 : {
404 : long i;
405 4962179 : GEN A = cgetg(lx, t_COL);
406 536838622 : for (i=1; i<lx; i++) gel(A,i) = gsub(gel(x,i), gel(y,i));
407 4962179 : return A;
408 : }
409 : GEN
410 4181090 : RgC_sub(GEN x, GEN y) { return RgC_sub_i(x, y, lg(x)); }
411 : GEN
412 292150 : RgV_sub(GEN x, GEN y)
413 2515367 : { pari_APPLY_type(t_VEC, gsub(gel(x,i), gel(y,i))) }
414 :
415 : GEN
416 592187 : RgM_add(GEN x, GEN y)
417 : {
418 592187 : long lx = lg(x), l, j;
419 : GEN z;
420 592187 : if (lx == 1) return cgetg(1, t_MAT);
421 592187 : z = cgetg(lx, t_MAT); l = lgcols(x);
422 3389689 : for (j = 1; j < lx; j++) gel(z,j) = RgC_add_i(gel(x,j), gel(y,j), l);
423 592187 : return z;
424 : }
425 : GEN
426 148906 : RgM_sub(GEN x, GEN y)
427 : {
428 148906 : long lx = lg(x), l, j;
429 : GEN z;
430 148906 : if (lx == 1) return cgetg(1, t_MAT);
431 148906 : z = cgetg(lx, t_MAT); l = lgcols(x);
432 929995 : for (j = 1; j < lx; j++) gel(z,j) = RgC_sub_i(gel(x,j), gel(y,j), l);
433 148906 : return z;
434 : }
435 :
436 : static GEN
437 3435861 : RgC_neg_i(GEN x, long lx)
438 : {
439 : long i;
440 3435861 : GEN y = cgetg(lx, t_COL);
441 26113607 : for (i=1; i<lx; i++) gel(y,i) = gneg(gel(x,i));
442 3435861 : return y;
443 : }
444 : GEN
445 880526 : RgC_neg(GEN x) { return RgC_neg_i(x, lg(x)); }
446 : GEN
447 4009 : RgV_neg(GEN x)
448 52505 : { pari_APPLY_type(t_VEC, gneg(gel(x,i))) }
449 : GEN
450 463350 : RgM_neg(GEN x)
451 : {
452 463350 : long i, hx, lx = lg(x);
453 463350 : GEN y = cgetg(lx, t_MAT);
454 463350 : if (lx == 1) return y;
455 463344 : hx = lgcols(x);
456 3018679 : for (i=1; i<lx; i++) gel(y,i) = RgC_neg_i(gel(x,i), hx);
457 463344 : return y;
458 : }
459 :
460 : GEN
461 4285184 : RgV_RgC_mul(GEN x, GEN y)
462 : {
463 4285184 : long lx = lg(x);
464 4285184 : if (lx != lg(y)) pari_err_OP("operation 'RgV_RgC_mul'", x, y);
465 4285124 : return RgV_dotproduct_i(x, y, lx);
466 : }
467 : GEN
468 1446 : RgC_RgV_mul(GEN x, GEN y)
469 : {
470 1446 : long i, ly = lg(y);
471 1446 : GEN z = cgetg(ly,t_MAT);
472 4338 : for (i=1; i<ly; i++) gel(z,i) = RgC_Rg_mul(x, gel(y,i));
473 1446 : return z;
474 : }
475 : GEN
476 0 : RgC_RgM_mul(GEN x, GEN y)
477 : {
478 0 : long i, ly = lg(y);
479 0 : GEN z = cgetg(ly,t_MAT);
480 0 : if (ly != 1 && lgcols(y) != 2) pari_err_OP("operation 'RgC_RgM_mul'",x,y);
481 0 : for (i=1; i<ly; i++) gel(z,i) = RgC_Rg_mul(x, gcoeff(y,1,i));
482 0 : return z;
483 : }
484 : GEN
485 0 : RgM_RgV_mul(GEN x, GEN y)
486 : {
487 0 : if (lg(x) != 2) pari_err_OP("operation 'RgM_RgV_mul'", x,y);
488 0 : return RgC_RgV_mul(gel(x,1), y);
489 : }
490 :
491 : /* x[i,]*y, l = lg(y) > 1 */
492 : static GEN
493 99132113 : RgMrow_RgC_mul_i(GEN x, GEN y, long i, long l)
494 : {
495 99132113 : pari_sp av = avma;
496 99132113 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
497 : long j;
498 834473876 : for (j=2; j<l; j++)
499 : {
500 735341763 : GEN c = gcoeff(x,i,j);
501 735341763 : if (!isintzero(c)) t = gadd(t, gmul(c, gel(y,j)));
502 : }
503 99132113 : return gc_upto(av,t);
504 : }
505 : GEN
506 4110 : RgMrow_RgC_mul(GEN x, GEN y, long i)
507 4110 : { return RgMrow_RgC_mul_i(x, y, i, lg(x)); }
508 :
509 : static GEN
510 24 : RgM_RgC_mul_FpM(GEN x, GEN y, GEN p)
511 : {
512 24 : pari_sp av = avma;
513 : GEN r;
514 24 : if (lgefint(p) == 3)
515 : {
516 12 : ulong pp = uel(p, 2);
517 12 : r = Flm_Flc_mul(RgM_to_Flm(x, pp), RgV_to_Flv(y, pp), pp);
518 12 : r = Flc_to_ZC_inplace(r);
519 : }
520 : else
521 12 : r = FpM_FpC_mul(RgM_to_FpM(x, p), RgC_to_FpC(y, p), p);
522 24 : return gc_upto(av, FpC_to_mod(r, p));
523 : }
524 :
525 : static GEN
526 24 : RgM_RgC_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
527 : {
528 24 : pari_sp av = avma;
529 24 : GEN b, T = RgX_to_FpX(pol, p);
530 24 : if (signe(T) == 0) pari_err_OP("*", x, y);
531 24 : b = FqM_FqC_mul(RgM_to_FqM(x, T, p), RgC_to_FqC(y, T, p), T, p);
532 24 : return gc_upto(av, FqC_to_mod(b, T, p));
533 : }
534 :
535 : static GEN
536 27202337 : RgM_RgC_mul_fast(GEN x, GEN y)
537 : {
538 : GEN p, pol;
539 : long pa;
540 27202337 : long t = RgM_RgC_type(x,y, &p,&pol,&pa);
541 27202337 : switch(t)
542 : {
543 7579574 : case t_INT: return ZM_ZC_mul(x,y);
544 108316 : case t_FRAC: return QM_QC_mul(x,y);
545 78 : case t_FFELT: return FFM_FFC_mul(x, y, pol);
546 24 : case t_INTMOD: return RgM_RgC_mul_FpM(x, y, p);
547 24 : case RgX_type_code(t_POLMOD, t_INTMOD):
548 24 : return RgM_RgC_mul_FqM(x, y, pol, p);
549 19514321 : default: return NULL;
550 : }
551 : }
552 :
553 : /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
554 : static GEN
555 29329347 : RgM_RgC_mul_i(GEN x, GEN y, long lx, long l)
556 : {
557 29329347 : GEN z = cgetg(l,t_COL);
558 : long i;
559 128457350 : for (i=1; i<l; i++) gel(z,i) = RgMrow_RgC_mul_i(x,y,i,lx);
560 29329347 : return z;
561 : }
562 :
563 : GEN
564 27202337 : RgM_RgC_mul(GEN x, GEN y)
565 : {
566 27202337 : long lx = lg(x);
567 : GEN z;
568 27202337 : if (lx != lg(y)) pari_err_OP("operation 'RgM_RgC_mul'", x,y);
569 27202337 : if (lx == 1) return cgetg(1,t_COL);
570 27202337 : z = RgM_RgC_mul_fast(x, y);
571 27202337 : if (z) return z;
572 19514321 : return RgM_RgC_mul_i(x, y, lx, lgcols(x));
573 : }
574 :
575 : GEN
576 228914 : RgV_RgM_mul(GEN x, GEN y)
577 : {
578 228914 : long i, lx, ly = lg(y);
579 : GEN z;
580 228914 : if (ly == 1) return cgetg(1,t_VEC);
581 228909 : lx = lg(x);
582 228909 : if (lx != lgcols(y)) pari_err_OP("operation 'RgV_RgM_mul'", x,y);
583 228904 : z = cgetg(ly, t_VEC);
584 2532673 : for (i=1; i<ly; i++) gel(z,i) = RgV_dotproduct_i(x, gel(y,i), lx);
585 228904 : return z;
586 : }
587 :
588 : static GEN
589 48984 : RgM_mul_FpM_i(GEN x, GEN y, GEN p)
590 : {
591 48984 : if (lgefint(p) == 3)
592 : {
593 48922 : ulong pp = uel(p, 2);
594 48922 : if (pp==2) return F2m_to_mod(F2m_mul(RgM_to_F2m(x), RgM_to_F2m(y)));
595 48880 : if (pp==3) return F3m_to_mod(F3m_mul(RgM_to_F3m(x), RgM_to_F3m(y)));
596 48748 : return Flm_to_mod(Flm_mul(RgM_to_Flm(x, pp), RgM_to_Flm(y, pp), pp), pp);
597 : }
598 62 : return FpM_to_mod(FpM_mul(RgM_to_FpM(x,p), RgM_to_FpM(y,p), p), p);
599 : }
600 : static GEN
601 48984 : RgM_mul_FpM(GEN x, GEN y, GEN p)
602 48984 : { pari_sp av = avma; return gc_upto(av, RgM_mul_FpM_i(x, y, p)); }
603 :
604 : static GEN
605 22290 : RgM_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
606 : {
607 22290 : pari_sp av = avma;
608 22290 : GEN b, T = RgX_to_FpX(pol, p);
609 22290 : if (signe(T) == 0) pari_err_OP("*", x, y);
610 22290 : b = FqM_mul(RgM_to_FqM(x, T, p), RgM_to_FqM(y, T, p), T, p);
611 22290 : return gc_upto(av, FqM_to_mod(b, T, p));
612 : }
613 :
614 : static GEN
615 10888 : RgM_liftred(GEN x, GEN T)
616 10888 : { return RgXQM_red(liftpol_shallow(x), T); }
617 :
618 : static GEN
619 1276 : RgM_mul_ZXQM(GEN x, GEN y, GEN T)
620 : {
621 1276 : pari_sp av = avma;
622 1276 : GEN b = ZXQM_mul(RgM_liftred(x,T), RgM_liftred(y, T), T);
623 1276 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
624 : }
625 :
626 : static GEN
627 114 : RgM_sqr_ZXQM(GEN x, GEN T)
628 : {
629 114 : pari_sp av = avma;
630 114 : GEN b = ZXQM_sqr(RgM_liftred(x, T), T);
631 114 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
632 : }
633 :
634 : static GEN
635 4108 : RgM_mul_QXQM(GEN x, GEN y, GEN T)
636 : {
637 4108 : pari_sp av = avma;
638 4108 : GEN b = QXQM_mul(RgM_liftred(x, T), RgM_liftred(y, T), T);
639 4108 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
640 : }
641 :
642 : static GEN
643 6 : RgM_sqr_QXQM(GEN x, GEN T)
644 : {
645 6 : pari_sp av = avma;
646 6 : GEN b = QXQM_sqr(RgM_liftred(x, T), T);
647 6 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
648 : }
649 :
650 : INLINE int
651 4150 : RgX_is_monic_ZX(GEN pol)
652 4150 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
653 :
654 : static GEN
655 6367454 : RgM_mul_fast(GEN x, GEN y)
656 : {
657 : GEN p, pol;
658 : long pa;
659 6367454 : long t = RgM_type2(x,y, &p,&pol,&pa);
660 6367454 : switch(t)
661 : {
662 2303976 : case t_INT: return ZM_mul(x,y);
663 100322 : case t_FRAC: return QM_mul(x,y);
664 3762 : case t_FFELT: return FFM_mul(x, y, pol);
665 48930 : case t_INTMOD: return RgM_mul_FpM(x, y, p);
666 1282 : case RgX_type_code(t_POLMOD, t_INT):
667 1282 : return ZX_is_monic(pol)? RgM_mul_ZXQM(x, y, pol): NULL;
668 4126 : case RgX_type_code(t_POLMOD, t_FRAC):
669 4126 : return RgX_is_monic_ZX(pol)? RgM_mul_QXQM(x, y, pol): NULL;
670 22290 : case RgX_type_code(t_POLMOD, t_INTMOD):
671 22290 : return RgM_mul_FqM(x, y, pol, p);
672 3882766 : default: return NULL;
673 : }
674 : }
675 :
676 : static GEN
677 1176 : RgM_sqr_fast(GEN x)
678 : {
679 : GEN p, pol;
680 : long pa;
681 1176 : long t = RgM_type(x, &p,&pol,&pa);
682 1176 : switch(t)
683 : {
684 204 : case t_INT: return ZM_sqr(x);
685 600 : case t_FRAC: return QM_sqr(x);
686 90 : case t_FFELT: return FFM_mul(x, x, pol);
687 54 : case t_INTMOD: return RgM_mul_FpM(x, x, p);
688 120 : case RgX_type_code(t_POLMOD, t_INT):
689 120 : return ZX_is_monic(pol)? RgM_sqr_ZXQM(x, pol): NULL;
690 24 : case RgX_type_code(t_POLMOD, t_FRAC):
691 24 : return RgX_is_monic_ZX(pol)? RgM_sqr_QXQM(x, pol): NULL;
692 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
693 0 : return RgM_mul_FqM(x, x, pol, p);
694 84 : default: return NULL;
695 : }
696 : }
697 :
698 : /* lx, ly > 1 */
699 : static GEN
700 3882790 : RgM_mul_i(GEN x, GEN y, long lx, long ly)
701 : {
702 3882790 : GEN z = cgetg(ly, t_MAT);
703 3882790 : long j, l = lgcols(x);
704 13697498 : for (j = 1; j < ly; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(y,j), lx, l);
705 3882790 : return z;
706 : }
707 : GEN
708 6380978 : RgM_mul(GEN x, GEN y)
709 : {
710 6380978 : long lx, ly = lg(y);
711 : GEN z;
712 6380978 : if (ly == 1) return cgetg(1,t_MAT);
713 6367478 : lx = lg(x);
714 6367478 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
715 6367478 : if (lx == 1) return zeromat(0,ly-1);
716 6367454 : z = RgM_mul_fast(x, y);
717 6367454 : if (z) return z;
718 3882790 : return RgM_mul_i(x, y, lx, ly);
719 : }
720 :
721 : GEN
722 1206 : RgM_sqr(GEN x)
723 : {
724 1206 : long j, lx = lg(x);
725 : GEN z;
726 1206 : if (lx == 1) return cgetg(1, t_MAT);
727 1176 : if (lx != lgcols(x)) pari_err_OP("operation 'RgM_mul'", x,x);
728 1176 : z = RgM_sqr_fast(x);
729 1176 : if (z) return z;
730 108 : z = cgetg(lx, t_MAT);
731 426 : for (j=1; j<lx; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(x,j), lx, lx);
732 108 : return z;
733 : }
734 :
735 : /* assume result is symmetric */
736 : GEN
737 0 : RgM_multosym(GEN x, GEN y)
738 : {
739 0 : long j, lx, ly = lg(y);
740 : GEN M;
741 0 : if (ly == 1) return cgetg(1,t_MAT);
742 0 : lx = lg(x);
743 0 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_multosym'", x,y);
744 0 : if (lx == 1) return cgetg(1,t_MAT);
745 0 : if (ly != lgcols(x)) pari_err_OP("operation 'RgM_multosym'", x,y);
746 0 : M = cgetg(ly, t_MAT);
747 0 : for (j=1; j<ly; j++)
748 : {
749 0 : GEN z = cgetg(ly,t_COL), yj = gel(y,j);
750 : long i;
751 0 : for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
752 0 : for (i=j; i<ly; i++)gel(z,i) = RgMrow_RgC_mul_i(x,yj,i,lx);
753 0 : gel(M,j) = z;
754 : }
755 0 : return M;
756 : }
757 : /* x~ * y, assuming result is symmetric */
758 : GEN
759 6507 : RgM_transmultosym(GEN x, GEN y)
760 : {
761 6507 : long i, j, l, ly = lg(y);
762 : GEN M;
763 6507 : if (ly == 1) return cgetg(1,t_MAT);
764 6507 : if (lg(x) != ly) pari_err_OP("operation 'RgM_transmultosym'", x,y);
765 6507 : l = lgcols(y);
766 6507 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmultosym'", x,y);
767 6507 : M = cgetg(ly, t_MAT);
768 26161 : for (i=1; i<ly; i++)
769 : {
770 19654 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
771 19654 : gel(M,i) = c;
772 40252 : for (j=1; j<i; j++)
773 20598 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(y,j),l);
774 19654 : gel(c,i) = RgV_dotproduct_i(xi,gel(y,i),l);
775 : }
776 6507 : return M;
777 : }
778 : /* x~ * y */
779 : GEN
780 0 : RgM_transmul(GEN x, GEN y)
781 : {
782 0 : long i, j, l, lx, ly = lg(y);
783 : GEN M;
784 0 : if (ly == 1) return cgetg(1,t_MAT);
785 0 : lx = lg(x);
786 0 : l = lgcols(y);
787 0 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmul'", x,y);
788 0 : M = cgetg(ly, t_MAT);
789 0 : for (i=1; i<ly; i++)
790 : {
791 0 : GEN yi = gel(y,i), c = cgetg(lx,t_COL);
792 0 : gel(M,i) = c;
793 0 : for (j=1; j<lx; j++) gel(c,j) = RgV_dotproduct_i(yi,gel(x,j),l);
794 : }
795 0 : return M;
796 : }
797 :
798 : GEN
799 4042331 : gram_matrix(GEN x)
800 : {
801 4042331 : long i,j, l, lx = lg(x);
802 : GEN M;
803 4042331 : if (!is_matvec_t(typ(x))) pari_err_TYPE("gram",x);
804 4042331 : if (lx == 1) return cgetg(1,t_MAT);
805 4042319 : l = lgcols(x);
806 4042319 : M = cgetg(lx,t_MAT);
807 12126938 : for (i=1; i<lx; i++)
808 : {
809 8084619 : GEN xi = gel(x,i), c = cgetg(lx,t_COL);
810 8084619 : gel(M,i) = c;
811 12126924 : for (j=1; j<i; j++)
812 4042305 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(x,j),l);
813 8084619 : gel(c,i) = RgV_dotsquare(xi);
814 : }
815 4042319 : return M;
816 : }
817 :
818 : static GEN
819 3288 : _RgM_add(void *E, GEN x, GEN y) { (void)E; return RgM_add(x, y); }
820 :
821 : static GEN
822 0 : _RgM_sub(void *E, GEN x, GEN y) { (void)E; return RgM_sub(x, y); }
823 :
824 : static GEN
825 4968 : _RgM_cmul(void *E, GEN P, long a, GEN x) { (void)E; return RgM_Rg_mul(x,gel(P,a+2)); }
826 :
827 : static GEN
828 204 : _RgM_sqr(void *E, GEN x) { (void) E; return RgM_sqr(x); }
829 :
830 : static GEN
831 654 : _RgM_mul(void *E, GEN x, GEN y) { (void) E; return RgM_mul(x, y); }
832 :
833 : static GEN
834 3618 : _RgM_one(void *E) { long *n = (long*) E; return matid(*n); }
835 :
836 : static GEN
837 0 : _RgM_zero(void *E) { long *n = (long*) E; return zeromat(*n,*n); }
838 :
839 : static GEN
840 2532 : _RgM_red(void *E, GEN x) { (void)E; return x; }
841 :
842 : static struct bb_algebra RgM_algebra = { _RgM_red, _RgM_add, _RgM_sub,
843 : _RgM_mul, _RgM_sqr, _RgM_one, _RgM_zero };
844 :
845 : /* generates the list of powers of x of degree 0,1,2,...,l*/
846 : GEN
847 144 : RgM_powers(GEN x, long l)
848 : {
849 144 : long n = lg(x)-1;
850 144 : return gen_powers(x,l,1,(void *) &n, &_RgM_sqr, &_RgM_mul, &_RgM_one);
851 : }
852 :
853 : GEN
854 420 : RgX_RgMV_eval(GEN Q, GEN x)
855 : {
856 420 : long n = lg(x)>1 ? lg(gel(x,1))-1:0;
857 420 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&n,&RgM_algebra,&_RgM_cmul);
858 : }
859 :
860 : GEN
861 1260 : RgX_RgM_eval(GEN Q, GEN x)
862 : {
863 1260 : long n = lg(x)-1;
864 1260 : return gen_bkeval(Q,degpol(Q),x,1,(void*)&n,&RgM_algebra,&_RgM_cmul);
865 : }
866 :
867 : GEN
868 3821308 : RgC_Rg_div(GEN x, GEN y)
869 16413265 : { pari_APPLY_type(t_COL, gdiv(gel(x,i),y)) }
870 :
871 : GEN
872 10154344 : RgC_Rg_mul(GEN x, GEN y)
873 64586940 : { pari_APPLY_type(t_COL, gmul(gel(x,i),y)) }
874 :
875 : GEN
876 33307 : RgV_Rg_mul(GEN x, GEN y)
877 1256228 : { pari_APPLY_type(t_VEC, gmul(gel(x,i),y)) }
878 :
879 : GEN
880 510295 : RgM_Rg_div(GEN X, GEN c) {
881 510295 : long i, j, h, l = lg(X);
882 510295 : GEN A = cgetg(l, t_MAT);
883 510295 : if (l == 1) return A;
884 510295 : h = lgcols(X);
885 2336177 : for (j=1; j<l; j++)
886 : {
887 1825882 : GEN a = cgetg(h, t_COL), x = gel(X, j);
888 14318754 : for (i = 1; i < h; i++) gel(a,i) = gdiv(gel(x,i), c);
889 1825882 : gel(A,j) = a;
890 : }
891 510295 : return A;
892 : }
893 : GEN
894 1168949 : RgM_Rg_mul(GEN X, GEN c) {
895 1168949 : long i, j, h, l = lg(X);
896 1168949 : GEN A = cgetg(l, t_MAT);
897 1168949 : if (l == 1) return A;
898 1168949 : h = lgcols(X);
899 3795131 : for (j=1; j<l; j++)
900 : {
901 2626182 : GEN a = cgetg(h, t_COL), x = gel(X, j);
902 11640718 : for (i = 1; i < h; i++) gel(a,i) = gmul(gel(x,i), c);
903 2626182 : gel(A,j) = a;
904 : }
905 1168949 : return A;
906 : }
907 :
908 : /********************************************************************/
909 : /* */
910 : /* SCALAR TO MATRIX/VECTOR */
911 : /* */
912 : /********************************************************************/
913 : /* fill the square nxn matrix equal to t*Id */
914 : static void
915 11009905 : fill_scalmat(GEN y, GEN t, long n)
916 : {
917 : long i;
918 41175485 : for (i = 1; i <= n; i++)
919 : {
920 30165580 : gel(y,i) = zerocol(n);
921 30165580 : gcoeff(y,i,i) = t;
922 : }
923 11009905 : }
924 :
925 : GEN
926 740916 : scalarmat(GEN x, long n) {
927 740916 : GEN y = cgetg(n+1, t_MAT);
928 740916 : if (!n) return y;
929 740916 : fill_scalmat(y, gcopy(x), n); return y;
930 : }
931 : GEN
932 2866723 : scalarmat_shallow(GEN x, long n) {
933 2866723 : GEN y = cgetg(n+1, t_MAT);
934 2866723 : fill_scalmat(y, x, n); return y;
935 : }
936 : GEN
937 180 : scalarmat_s(long x, long n) {
938 180 : GEN y = cgetg(n+1, t_MAT);
939 180 : if (!n) return y;
940 180 : fill_scalmat(y, stoi(x), n); return y;
941 : }
942 : GEN
943 7402092 : matid(long n) {
944 : GEN y;
945 7402092 : if (n < 0) pari_err_DOMAIN("matid", "size", "<", gen_0, stoi(n));
946 7402086 : y = cgetg(n+1, t_MAT);
947 7402086 : fill_scalmat(y, gen_1, n); return y;
948 : }
949 :
950 : INLINE GEN
951 1946510 : scalarcol_i(GEN x, long n, long c)
952 : {
953 : long i;
954 1946510 : GEN y = cgetg(n+1,t_COL);
955 1946510 : if (!n) return y;
956 1946510 : gel(y,1) = c? gcopy(x): x;
957 5786209 : for (i=2; i<=n; i++) gel(y,i) = gen_0;
958 1946510 : return y;
959 : }
960 :
961 : GEN
962 669788 : scalarcol(GEN x, long n) { return scalarcol_i(x,n,1); }
963 :
964 : GEN
965 1276722 : scalarcol_shallow(GEN x, long n) { return scalarcol_i(x,n,0); }
966 :
967 : int
968 14782 : RgM_isscalar(GEN x, GEN s)
969 : {
970 14782 : long i, j, lx = lg(x);
971 :
972 14782 : if (lx == 1) return 1;
973 14782 : if (lx != lgcols(x)) return 0;
974 14782 : if (!s) s = gcoeff(x,1,1);
975 :
976 41640 : for (j=1; j<lx; j++)
977 : {
978 33339 : GEN c = gel(x,j);
979 116481 : for (i=1; i<j; )
980 88543 : if (!gequal0(gel(c,i++))) return 0;
981 : /* i = j */
982 27938 : if (!gequal(gel(c,i++),s)) return 0;
983 132256 : for ( ; i<lx; )
984 105398 : if (!gequal0(gel(c,i++))) return 0;
985 : }
986 8301 : return 1;
987 : }
988 :
989 : int
990 2932 : RgM_isidentity(GEN x)
991 : {
992 2932 : long i,j, lx = lg(x);
993 :
994 2932 : if (lx == 1) return 1;
995 2932 : if (lx != lgcols(x)) return 0;
996 6971 : for (j=1; j<lx; j++)
997 : {
998 6732 : GEN c = gel(x,j);
999 29568 : for (i=1; i<j; )
1000 24298 : if (!gequal0(gel(c,i++))) return 0;
1001 : /* i = j */
1002 5270 : if (!gequal1(gel(c,i++))) return 0;
1003 36446 : for ( ; i<lx; )
1004 32407 : if (!gequal0(gel(c,i++))) return 0;
1005 : }
1006 239 : return 1;
1007 : }
1008 :
1009 : long
1010 552 : RgC_is_ei(GEN x)
1011 : {
1012 552 : long i, j = 0, l = lg(x);
1013 2976 : for (i = 1; i < l; i++)
1014 : {
1015 2424 : GEN c = gel(x,i);
1016 2424 : if (gequal0(c)) continue;
1017 552 : if (!gequal1(c) || j) return 0;
1018 552 : j = i;
1019 : }
1020 552 : return j;
1021 : }
1022 :
1023 : int
1024 287 : RgM_isdiagonal(GEN x)
1025 : {
1026 287 : long i,j, lx = lg(x);
1027 287 : if (lx == 1) return 1;
1028 287 : if (lx != lgcols(x)) return 0;
1029 :
1030 2757 : for (j=1; j<lx; j++)
1031 : {
1032 2476 : GEN c = gel(x,j);
1033 15813 : for (i=1; i<j; i++)
1034 13337 : if (!gequal0(gel(c,i))) return 0;
1035 15813 : for (i++; i<lx; i++)
1036 13343 : if (!gequal0(gel(c,i))) return 0;
1037 : }
1038 281 : return 1;
1039 : }
1040 : int
1041 270 : isdiagonal(GEN x) { return (typ(x)==t_MAT) && RgM_isdiagonal(x); }
1042 :
1043 : GEN
1044 19365 : RgM_det_triangular(GEN mat)
1045 : {
1046 19365 : long i,l = lg(mat);
1047 : pari_sp av;
1048 : GEN s;
1049 :
1050 19365 : if (l<3) return l<2? gen_1: gcopy(gcoeff(mat,1,1));
1051 17727 : av = avma; s = gcoeff(mat,1,1);
1052 103861 : for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1053 17727 : return av==avma? gcopy(s): gc_upto(av,s);
1054 : }
1055 :
1056 : GEN
1057 405172 : RgV_kill0(GEN v)
1058 : {
1059 : long i, l;
1060 405172 : GEN w = cgetg_copy(v, &l);
1061 108965060 : for (i = 1; i < l; i++)
1062 : {
1063 108559888 : GEN a = gel(v,i);
1064 108559888 : gel(w,i) = gequal0(a) ? NULL: a;
1065 : }
1066 405172 : return w;
1067 : }
|
