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 1362759 : RgM_is_ZM(GEN x)
20 : {
21 1362759 : long i, j, h, l = lg(x);
22 1362759 : if (l == 1) return 1;
23 1362080 : h = lgcols(x);
24 1362080 : if (h == 1) return 1;
25 5105074 : for (j = l-1; j > 0; j--)
26 23593040 : for (i = h-1; i > 0; i--)
27 19849885 : if (typ(gcoeff(x,i,j)) != t_INT) return 0;
28 1197001 : return 1;
29 : }
30 :
31 : int
32 203 : RgM_is_QM(GEN x)
33 : {
34 203 : long i, j, h, l = lg(x);
35 203 : if (l == 1) return 1;
36 189 : h = lgcols(x);
37 189 : if (h == 1) return 1;
38 1218 : for (j = l-1; j > 0; j--)
39 15057 : for (i = h-1; i > 0; i--)
40 14028 : if (!is_rational_t(typ(gcoeff(x,i,j)))) return 0;
41 175 : return 1;
42 : }
43 :
44 : int
45 49 : RgV_is_ZMV(GEN V)
46 : {
47 49 : long i, l = lg(V);
48 399 : for (i=1; i<l; i++)
49 350 : if (typ(gel(V,i))!=t_MAT || !RgM_is_ZM(gel(V,i)))
50 0 : return 0;
51 49 : 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 8373768 : RgMrow_zc_mul_i(GEN x, GEN y, long c, long i)
64 : {
65 8373768 : pari_sp av = avma;
66 8373768 : GEN s = NULL;
67 : long j;
68 829158209 : for (j=1; j<c; j++)
69 : {
70 820784441 : long t = y[j];
71 820784441 : if (!t) continue;
72 28419149 : if (!s) { s = gmulgs(gcoeff(x,i,j),t); continue; }
73 20260893 : switch(t)
74 : {
75 14798162 : case 1: s = gadd(s, gcoeff(x,i,j)); break;
76 239730 : case -1: s = gsub(s, gcoeff(x,i,j)); break;
77 5223001 : default: s = gadd(s, gmulgs(gcoeff(x,i,j), t)); break;
78 : }
79 : }
80 8373768 : return s? gc_upto(av, s): gc_const(av, gen_0);
81 : }
82 : GEN
83 70028 : 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 242996 : RgM_zc_mul_i(GEN x, GEN y, long c, long l)
87 : {
88 242996 : GEN z = cgetg(l,t_COL);
89 : long i;
90 8546736 : for (i = 1; i < l; i++) gel(z,i) = RgMrow_zc_mul_i(x,y,c,i);
91 242996 : return z;
92 : }
93 : GEN
94 71981 : 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 57925 : RgM_zm_mul(GEN x, GEN y)
98 : {
99 57925 : long j, c, l = lg(x), ly = lg(y);
100 57925 : GEN z = cgetg(ly, t_MAT);
101 57925 : if (l == 1) return z;
102 57925 : c = lgcols(x);
103 228940 : for (j = 1; j < ly; j++) gel(z,j) = RgM_zc_mul_i(x, gel(y,j), l,c);
104 57925 : return z;
105 : }
106 :
107 : /* x[i,]*y, l = lg(y) > 1 */
108 : static GEN
109 32961182 : RgMrow_ZC_mul_i(GEN x, GEN y, long i, long l)
110 : {
111 32961182 : pari_sp av = avma;
112 32961182 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
113 : long j;
114 5809762156 : for (j=2; j<l; j++)
115 5776800974 : if (signe(gel(y,j))) t = gadd(t, gmul(gcoeff(x,i,j), gel(y,j)));
116 32961182 : 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 1942830 : RgM_ZC_mul_i(GEN x, GEN y, long lx, long l)
122 : {
123 1942830 : GEN z = cgetg(l,t_COL);
124 : long i;
125 34904012 : for (i=1; i<l; i++) gel(z,i) = RgMrow_ZC_mul_i(x,y,i,lx);
126 1942830 : return z;
127 : }
128 :
129 : GEN
130 1942830 : RgM_ZM_mul_worker(GEN y, GEN x)
131 1942830 : { return RgM_ZC_mul_i(x, y, lg(x), lgcols(x)); }
132 :
133 : /* mostly useful when y is sparse */
134 : GEN
135 356708 : RgM_ZM_mul(GEN x, GEN y)
136 : {
137 356708 : pari_sp av = avma;
138 : GEN worker;
139 356708 : if (lg(x) == 1) return cgetg(lg(y), t_MAT);
140 356708 : worker = snm_closure(is_entry("_RgM_ZM_mul_worker"),mkvec(x));
141 356708 : return gc_upto(av, gen_parapply(worker, y));
142 : }
143 :
144 : static GEN
145 293289 : RgV_zc_mul_i(GEN x, GEN y, long l)
146 : {
147 : long i;
148 293289 : GEN z = gen_0;
149 293289 : pari_sp av = avma;
150 9624631 : for (i = 1; i < l; i++) z = gadd(z, gmulgs(gel(x,i), y[i]));
151 293289 : return gc_upto(av, z);
152 : }
153 : GEN
154 83454 : RgV_zc_mul(GEN x, GEN y) { return RgV_zc_mul_i(x, y, lg(x)); }
155 :
156 : GEN
157 58931 : RgV_zm_mul(GEN x, GEN y)
158 : {
159 58931 : long j, l = lg(x), ly = lg(y);
160 58931 : GEN z = cgetg(ly, t_VEC);
161 268766 : for (j = 1; j < ly; j++) gel(z,j) = RgV_zc_mul_i(x, gel(y,j), l);
162 58931 : return z;
163 : }
164 :
165 : /* scalar product x.x */
166 : GEN
167 9756425 : RgV_dotsquare(GEN x)
168 : {
169 9756425 : long i, lx = lg(x);
170 9756425 : pari_sp av = avma;
171 : GEN z;
172 9756425 : if (lx == 1) return gen_0;
173 9756425 : z = gsqr(gel(x,1));
174 21796204 : for (i=2; i<lx; i++)
175 : {
176 12039779 : z = gadd(z, gsqr(gel(x,i)));
177 12039779 : 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 9756425 : return gc_upto(av,z);
184 : }
185 :
186 : /* scalar product x.y, lx = lg(x) = lg(y) */
187 : static GEN
188 15980052 : RgV_dotproduct_i(GEN x, GEN y, long lx)
189 : {
190 15980052 : pari_sp av = avma;
191 : long i;
192 : GEN z;
193 15980052 : if (lx == 1) return gen_0;
194 15979653 : z = gmul(gel(x,1),gel(y,1));
195 253551363 : for (i=2; i<lx; i++)
196 : {
197 237571710 : z = gadd(z, gmul(gel(x,i), gel(y,i)));
198 237571710 : 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 15979653 : return gc_upto(av,z);
205 : }
206 : GEN
207 1031617 : RgV_dotproduct(GEN x,GEN y)
208 : {
209 1031617 : if (x == y) return RgV_dotsquare(x);
210 1031617 : return RgV_dotproduct_i(x, y, lg(x));
211 : }
212 : /* v[1] + ... + v[lg(v)-1] */
213 : GEN
214 3192968 : RgV_sum(GEN v)
215 : {
216 : GEN p;
217 3192968 : long i, l = lg(v);
218 3192968 : if (l == 1) return gen_0;
219 8062742 : p = gel(v,1); for (i=2; i<l; i++) p = gadd(p, gel(v,i));
220 3192968 : 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 400 : RgM_sumcol(GEN A)
244 : {
245 400 : long i,j,m,l = lg(A);
246 : GEN v;
247 :
248 400 : if (l == 1) return cgetg(1,t_MAT);
249 400 : if (l == 2) return gcopy(gel(A,1));
250 216 : m = lgcols(A);
251 216 : v = cgetg(m, t_COL);
252 718 : for (i = 1; i < m; i++)
253 : {
254 502 : pari_sp av = avma;
255 502 : GEN s = gcoeff(A,i,1);
256 1214 : for (j = 2; j < l; j++) s = gadd(s, gcoeff(A,i,j));
257 502 : gel(v, i) = gc_upto(av, s);
258 : }
259 216 : return v;
260 : }
261 :
262 : static GEN
263 3395912 : _gmul(void *data, GEN x, GEN y)
264 3395912 : { (void)data; return gmul(x,y); }
265 :
266 : GEN
267 458545 : RgV_prod(GEN x)
268 : {
269 458545 : 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 7878 : RgM_Rg_add(GEN x, GEN y)
276 : {
277 7878 : long l = lg(x), i, j;
278 7878 : GEN z = cgetg(l,t_MAT);
279 :
280 7878 : if (l==1) return z;
281 7878 : if (l != lgcols(x)) pari_err_OP("+", x, y);
282 57833 : for (i=1; i<l; i++)
283 : {
284 49955 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
285 49955 : gel(z,i) = zi;
286 2057330 : for (j=1; j<l; j++)
287 2007375 : gel(zi,j) = i==j? gadd(y,gel(xi,j)): gcopy(gel(xi,j));
288 : }
289 7878 : 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 739 : RgM_Rg_add_shallow(GEN x, GEN y)
310 : {
311 739 : long l = lg(x), i, j;
312 739 : GEN z = cgetg(l,t_MAT);
313 :
314 739 : if (l==1) return z;
315 739 : if (l != lgcols(x)) pari_err_OP( "+", x, y);
316 2384 : for (i=1; i<l; i++)
317 : {
318 1645 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
319 1645 : gel(z,i) = zi;
320 9456 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
321 1645 : gel(zi,i) = gadd(gel(zi,i), y);
322 : }
323 739 : return z;
324 : }
325 : GEN
326 316623 : RgM_Rg_sub_shallow(GEN x, GEN y)
327 : {
328 316623 : long l = lg(x), i, j;
329 316623 : GEN z = cgetg(l,t_MAT);
330 :
331 316623 : if (l==1) return z;
332 316623 : if (l != lgcols(x)) pari_err_OP( "-", x, y);
333 2043211 : for (i=1; i<l; i++)
334 : {
335 1726588 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
336 1726588 : gel(z,i) = zi;
337 27515378 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
338 1726588 : gel(zi,i) = gsub(gel(zi,i), y);
339 : }
340 316623 : return z;
341 : }
342 :
343 : GEN
344 2519686 : RgC_Rg_add(GEN x, GEN y)
345 : {
346 2519686 : long k, lx = lg(x);
347 2519686 : GEN z = cgetg(lx, t_COL);
348 2519686 : if (lx == 1)
349 : {
350 0 : if (isintzero(y)) return z;
351 0 : pari_err_TYPE2("+",x,y);
352 : }
353 2519686 : gel(z,1) = gadd(y,gel(x,1));
354 7050551 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
355 2519686 : return z;
356 : }
357 : GEN
358 24155 : RgC_Rg_sub(GEN x, GEN y)
359 : {
360 24155 : long k, lx = lg(x);
361 24155 : GEN z = cgetg(lx, t_COL);
362 24155 : if (lx == 1)
363 : {
364 0 : if (isintzero(y)) return z;
365 0 : pari_err_TYPE2("-",x,y);
366 : }
367 24155 : gel(z,1) = gsub(gel(x,1), y);
368 69449 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
369 24155 : return z;
370 : }
371 : /* a - x */
372 : GEN
373 32766 : Rg_RgC_sub(GEN a, GEN x)
374 : {
375 32766 : long k, lx = lg(x);
376 32766 : GEN z = cgetg(lx,t_COL);
377 32766 : if (lx == 1)
378 : {
379 0 : if (isintzero(a)) return z;
380 0 : pari_err_TYPE2("-",a,x);
381 : }
382 32766 : gel(z,1) = gsub(a, gel(x,1));
383 74766 : for (k = 2; k < lx; k++) gel(z,k) = gneg(gel(x,k));
384 32766 : return z;
385 : }
386 :
387 : static GEN
388 19037458 : RgC_add_i(GEN x, GEN y, long lx)
389 : {
390 19037458 : GEN A = cgetg(lx, t_COL);
391 : long i;
392 143806327 : for (i=1; i<lx; i++) gel(A,i) = gadd(gel(x,i), gel(y,i));
393 19037458 : return A;
394 : }
395 : GEN
396 15773701 : RgC_add(GEN x, GEN y) { return RgC_add_i(x, y, lg(x)); }
397 : GEN
398 2606578 : RgV_add(GEN x, GEN y)
399 10799093 : { pari_APPLY_type(t_VEC, gadd(gel(x,i), gel(y,i))) }
400 :
401 : static GEN
402 5918758 : RgC_sub_i(GEN x, GEN y, long lx)
403 : {
404 : long i;
405 5918758 : GEN A = cgetg(lx, t_COL);
406 641069511 : for (i=1; i<lx; i++) gel(A,i) = gsub(gel(x,i), gel(y,i));
407 5918758 : return A;
408 : }
409 : GEN
410 5004477 : RgC_sub(GEN x, GEN y) { return RgC_sub_i(x, y, lg(x)); }
411 : GEN
412 381067 : RgV_sub(GEN x, GEN y)
413 3139640 : { pari_APPLY_type(t_VEC, gsub(gel(x,i), gel(y,i))) }
414 :
415 : GEN
416 690886 : RgM_add(GEN x, GEN y)
417 : {
418 690886 : long lx = lg(x), l, j;
419 : GEN z;
420 690886 : if (lx == 1) return cgetg(1, t_MAT);
421 690886 : z = cgetg(lx, t_MAT); l = lgcols(x);
422 3954643 : for (j = 1; j < lx; j++) gel(z,j) = RgC_add_i(gel(x,j), gel(y,j), l);
423 690886 : return z;
424 : }
425 : GEN
426 174241 : RgM_sub(GEN x, GEN y)
427 : {
428 174241 : long lx = lg(x), l, j;
429 : GEN z;
430 174241 : if (lx == 1) return cgetg(1, t_MAT);
431 174241 : z = cgetg(lx, t_MAT); l = lgcols(x);
432 1088522 : for (j = 1; j < lx; j++) gel(z,j) = RgC_sub_i(gel(x,j), gel(y,j), l);
433 174241 : return z;
434 : }
435 :
436 : static GEN
437 4040633 : RgC_neg_i(GEN x, long lx)
438 : {
439 : long i;
440 4040633 : GEN y = cgetg(lx, t_COL);
441 30590643 : for (i=1; i<lx; i++) gel(y,i) = gneg(gel(x,i));
442 4040633 : return y;
443 : }
444 : GEN
445 1052504 : RgC_neg(GEN x) { return RgC_neg_i(x, lg(x)); }
446 : GEN
447 5271 : RgV_neg(GEN x)
448 63357 : { pari_APPLY_type(t_VEC, gneg(gel(x,i))) }
449 : GEN
450 542310 : RgM_neg(GEN x)
451 : {
452 542310 : long i, hx, lx = lg(x);
453 542310 : GEN y = cgetg(lx, t_MAT);
454 542310 : if (lx == 1) return y;
455 542303 : hx = lgcols(x);
456 3530432 : for (i=1; i<lx; i++) gel(y,i) = RgC_neg_i(gel(x,i), hx);
457 542303 : return y;
458 : }
459 :
460 : GEN
461 5127157 : RgV_RgC_mul(GEN x, GEN y)
462 : {
463 5127157 : long lx = lg(x);
464 5127157 : if (lx != lg(y)) pari_err_OP("operation 'RgV_RgC_mul'", x, y);
465 5127073 : return RgV_dotproduct_i(x, y, lx);
466 : }
467 : GEN
468 1687 : RgC_RgV_mul(GEN x, GEN y)
469 : {
470 1687 : long i, ly = lg(y);
471 1687 : GEN z = cgetg(ly,t_MAT);
472 5061 : for (i=1; i<ly; i++) gel(z,i) = RgC_Rg_mul(x, gel(y,i));
473 1687 : 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 118855310 : RgMrow_RgC_mul_i(GEN x, GEN y, long i, long l)
494 : {
495 118855310 : pari_sp av = avma;
496 118855310 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
497 : long j;
498 1004505590 : for (j=2; j<l; j++)
499 : {
500 885650280 : GEN c = gcoeff(x,i,j);
501 885650280 : if (!isintzero(c)) t = gadd(t, gmul(c, gel(y,j)));
502 : }
503 118855310 : return gc_upto(av,t);
504 : }
505 : GEN
506 4795 : RgMrow_RgC_mul(GEN x, GEN y, long i)
507 4795 : { return RgMrow_RgC_mul_i(x, y, i, lg(x)); }
508 :
509 : static GEN
510 28 : RgM_RgC_mul_FpM(GEN x, GEN y, GEN p)
511 : {
512 28 : pari_sp av = avma;
513 : GEN r;
514 28 : if (lgefint(p) == 3)
515 : {
516 14 : ulong pp = uel(p, 2);
517 14 : r = Flm_Flc_mul(RgM_to_Flm(x, pp), RgV_to_Flv(y, pp), pp);
518 14 : r = Flc_to_ZC_inplace(r);
519 : }
520 : else
521 14 : r = FpM_FpC_mul(RgM_to_FpM(x, p), RgC_to_FpC(y, p), p);
522 28 : return gc_upto(av, FpC_to_mod(r, p));
523 : }
524 :
525 : static GEN
526 28 : RgM_RgC_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
527 : {
528 28 : pari_sp av = avma;
529 28 : GEN b, T = RgX_to_FpX(pol, p);
530 28 : if (signe(T) == 0) pari_err_OP("*", x, y);
531 28 : b = FqM_FqC_mul(RgM_to_FqM(x, T, p), RgC_to_FqC(y, T, p), T, p);
532 28 : return gc_upto(av, FqC_to_mod(b, T, p));
533 : }
534 :
535 : static GEN
536 32245254 : RgM_RgC_mul_fast(GEN x, GEN y)
537 : {
538 : GEN p, pol;
539 : long pa;
540 32245254 : long t = RgM_RgC_type(x,y, &p,&pol,&pa);
541 32245254 : switch(t)
542 : {
543 8852351 : case t_INT: return ZM_ZC_mul(x,y);
544 127070 : case t_FRAC: return QM_QC_mul(x,y);
545 91 : case t_FFELT: return FFM_FFC_mul(x, y, pol);
546 28 : case t_INTMOD: return RgM_RgC_mul_FpM(x, y, p);
547 28 : case RgX_type_code(t_POLMOD, t_INTMOD):
548 28 : return RgM_RgC_mul_FqM(x, y, pol, p);
549 23265686 : 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 35195472 : RgM_RgC_mul_i(GEN x, GEN y, long lx, long l)
556 : {
557 35195472 : GEN z = cgetg(l,t_COL);
558 : long i;
559 154045987 : for (i=1; i<l; i++) gel(z,i) = RgMrow_RgC_mul_i(x,y,i,lx);
560 35195472 : return z;
561 : }
562 :
563 : GEN
564 32245254 : RgM_RgC_mul(GEN x, GEN y)
565 : {
566 32245254 : long lx = lg(x);
567 : GEN z;
568 32245254 : if (lx != lg(y)) pari_err_OP("operation 'RgM_RgC_mul'", x,y);
569 32245254 : if (lx == 1) return cgetg(1,t_COL);
570 32245254 : z = RgM_RgC_mul_fast(x, y);
571 32245254 : if (z) return z;
572 23265686 : return RgM_RgC_mul_i(x, y, lx, lgcols(x));
573 : }
574 :
575 : GEN
576 426340 : RgV_RgM_mul(GEN x, GEN y)
577 : {
578 426340 : long i, lx, ly = lg(y);
579 : GEN z;
580 426340 : if (ly == 1) return cgetg(1,t_VEC);
581 426333 : lx = lg(x);
582 426333 : if (lx != lgcols(y)) pari_err_OP("operation 'RgV_RgM_mul'", x,y);
583 426326 : z = cgetg(ly, t_VEC);
584 5356283 : for (i=1; i<ly; i++) gel(z,i) = RgV_dotproduct_i(x, gel(y,i), lx);
585 426326 : return z;
586 : }
587 :
588 : static GEN
589 57148 : RgM_mul_FpM_i(GEN x, GEN y, GEN p)
590 : {
591 57148 : if (lgefint(p) == 3)
592 : {
593 57076 : ulong pp = uel(p, 2);
594 57076 : if (pp==2) return F2m_to_mod(F2m_mul(RgM_to_F2m(x), RgM_to_F2m(y)));
595 57027 : if (pp==3) return F3m_to_mod(F3m_mul(RgM_to_F3m(x), RgM_to_F3m(y)));
596 56873 : return Flm_to_mod(Flm_mul(RgM_to_Flm(x, pp), RgM_to_Flm(y, pp), pp), pp);
597 : }
598 72 : return FpM_to_mod(FpM_mul(RgM_to_FpM(x,p), RgM_to_FpM(y,p), p), p);
599 : }
600 : static GEN
601 57148 : RgM_mul_FpM(GEN x, GEN y, GEN p)
602 57148 : { pari_sp av = avma; return gc_upto(av, RgM_mul_FpM_i(x, y, p)); }
603 :
604 : static GEN
605 26005 : RgM_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
606 : {
607 26005 : pari_sp av = avma;
608 26005 : GEN b, T = RgX_to_FpX(pol, p);
609 26005 : if (signe(T) == 0) pari_err_OP("*", x, y);
610 26005 : b = FqM_mul(RgM_to_FqM(x, T, p), RgM_to_FqM(y, T, p), T, p);
611 26005 : return gc_upto(av, FqM_to_mod(b, T, p));
612 : }
613 :
614 : static GEN
615 13398 : RgM_liftred(GEN x, GEN T)
616 13398 : { return RgXQM_red(liftpol_shallow(x), T); }
617 :
618 : static GEN
619 1505 : RgM_mul_ZXQM(GEN x, GEN y, GEN T)
620 : {
621 1505 : pari_sp av = avma;
622 1505 : GEN b = ZXQM_mul(RgM_liftred(x,T), RgM_liftred(y, T), T);
623 1505 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
624 : }
625 :
626 : static GEN
627 133 : RgM_sqr_ZXQM(GEN x, GEN T)
628 : {
629 133 : pari_sp av = avma;
630 133 : GEN b = ZXQM_sqr(RgM_liftred(x, T), T);
631 133 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
632 : }
633 :
634 : static GEN
635 5124 : RgM_mul_QXQM(GEN x, GEN y, GEN T)
636 : {
637 5124 : pari_sp av = avma;
638 5124 : GEN b = QXQM_mul(RgM_liftred(x, T), RgM_liftred(y, T), T);
639 5124 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
640 : }
641 :
642 : static GEN
643 7 : RgM_sqr_QXQM(GEN x, GEN T)
644 : {
645 7 : pari_sp av = avma;
646 7 : GEN b = QXQM_sqr(RgM_liftred(x, T), T);
647 7 : return gc_GEN(av, QXQM_to_mod_shallow(b,T));
648 : }
649 :
650 : INLINE int
651 5173 : RgX_is_monic_ZX(GEN pol)
652 5173 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
653 :
654 : static GEN
655 7628306 : RgM_mul_fast(GEN x, GEN y)
656 : {
657 : GEN p, pol;
658 : long pa;
659 7628306 : long t = RgM_type2(x,y, &p,&pol,&pa);
660 7628306 : switch(t)
661 : {
662 2690240 : case t_INT: return ZM_mul(x,y);
663 117684 : case t_FRAC: return QM_mul(x,y);
664 4389 : case t_FFELT: return FFM_mul(x, y, pol);
665 57085 : case t_INTMOD: return RgM_mul_FpM(x, y, p);
666 1512 : case RgX_type_code(t_POLMOD, t_INT):
667 1512 : return ZX_is_monic(pol)? RgM_mul_ZXQM(x, y, pol): NULL;
668 5145 : case RgX_type_code(t_POLMOD, t_FRAC):
669 5145 : return RgX_is_monic_ZX(pol)? RgM_mul_QXQM(x, y, pol): NULL;
670 26005 : case RgX_type_code(t_POLMOD, t_INTMOD):
671 26005 : return RgM_mul_FqM(x, y, pol, p);
672 4726246 : default: return NULL;
673 : }
674 : }
675 :
676 : static GEN
677 1372 : RgM_sqr_fast(GEN x)
678 : {
679 : GEN p, pol;
680 : long pa;
681 1372 : long t = RgM_type(x, &p,&pol,&pa);
682 1372 : switch(t)
683 : {
684 238 : case t_INT: return ZM_sqr(x);
685 700 : case t_FRAC: return QM_sqr(x);
686 105 : case t_FFELT: return FFM_mul(x, x, pol);
687 63 : case t_INTMOD: return RgM_mul_FpM(x, x, p);
688 140 : case RgX_type_code(t_POLMOD, t_INT):
689 140 : return ZX_is_monic(pol)? RgM_sqr_ZXQM(x, pol): NULL;
690 28 : case RgX_type_code(t_POLMOD, t_FRAC):
691 28 : 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 98 : default: return NULL;
695 : }
696 : }
697 :
698 : /* lx, ly > 1 */
699 : static GEN
700 4726274 : RgM_mul_i(GEN x, GEN y, long lx, long ly)
701 : {
702 4726274 : GEN z = cgetg(ly, t_MAT);
703 4726274 : long j, l = lgcols(x);
704 16655689 : for (j = 1; j < ly; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(y,j), lx, l);
705 4726274 : return z;
706 : }
707 : GEN
708 7644084 : RgM_mul(GEN x, GEN y)
709 : {
710 7644084 : long lx, ly = lg(y);
711 : GEN z;
712 7644084 : if (ly == 1) return cgetg(1,t_MAT);
713 7628334 : lx = lg(x);
714 7628334 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
715 7628334 : if (lx == 1) return zeromat(0,ly-1);
716 7628306 : z = RgM_mul_fast(x, y);
717 7628306 : if (z) return z;
718 4726274 : return RgM_mul_i(x, y, lx, ly);
719 : }
720 :
721 : GEN
722 1407 : RgM_sqr(GEN x)
723 : {
724 1407 : long j, lx = lg(x);
725 : GEN z;
726 1407 : if (lx == 1) return cgetg(1, t_MAT);
727 1372 : if (lx != lgcols(x)) pari_err_OP("operation 'RgM_mul'", x,x);
728 1372 : z = RgM_sqr_fast(x);
729 1372 : if (z) return z;
730 126 : z = cgetg(lx, t_MAT);
731 497 : for (j=1; j<lx; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(x,j), lx, lx);
732 126 : 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 7663 : RgM_transmultosym(GEN x, GEN y)
760 : {
761 7663 : long i, j, l, ly = lg(y);
762 : GEN M;
763 7663 : if (ly == 1) return cgetg(1,t_MAT);
764 7663 : if (lg(x) != ly) pari_err_OP("operation 'RgM_transmultosym'", x,y);
765 7663 : l = lgcols(y);
766 7663 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmultosym'", x,y);
767 7663 : M = cgetg(ly, t_MAT);
768 30820 : for (i=1; i<ly; i++)
769 : {
770 23157 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
771 23157 : gel(M,i) = c;
772 47490 : for (j=1; j<i; j++)
773 24333 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(y,j),l);
774 23157 : gel(c,i) = RgV_dotproduct_i(xi,gel(y,i),l);
775 : }
776 7663 : 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 4843943 : gram_matrix(GEN x)
800 : {
801 4843943 : long i,j, l, lx = lg(x);
802 : GEN M;
803 4843943 : if (!is_matvec_t(typ(x))) pari_err_TYPE("gram",x);
804 4843943 : if (lx == 1) return cgetg(1,t_MAT);
805 4843929 : l = lgcols(x);
806 4843929 : M = cgetg(lx,t_MAT);
807 14531766 : for (i=1; i<lx; i++)
808 : {
809 9687837 : GEN xi = gel(x,i), c = cgetg(lx,t_COL);
810 9687837 : gel(M,i) = c;
811 14531752 : for (j=1; j<i; j++)
812 4843915 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(x,j),l);
813 9687837 : gel(c,i) = RgV_dotsquare(xi);
814 : }
815 4843929 : return M;
816 : }
817 :
818 : static GEN
819 3836 : _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 5796 : _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 238 : _RgM_sqr(void *E, GEN x) { (void) E; return RgM_sqr(x); }
829 :
830 : static GEN
831 763 : _RgM_mul(void *E, GEN x, GEN y) { (void) E; return RgM_mul(x, y); }
832 :
833 : static GEN
834 4221 : _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 2954 : _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 168 : RgM_powers(GEN x, long l)
848 : {
849 168 : long n = lg(x)-1;
850 168 : return gen_powers(x,l,1,(void *) &n, &_RgM_sqr, &_RgM_mul, &_RgM_one);
851 : }
852 :
853 : GEN
854 490 : RgX_RgMV_eval(GEN Q, GEN x)
855 : {
856 490 : long n = lg(x)>1 ? lg(gel(x,1))-1:0;
857 490 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&n,&RgM_algebra,&_RgM_cmul);
858 : }
859 :
860 : GEN
861 1470 : RgX_RgM_eval(GEN Q, GEN x)
862 : {
863 1470 : long n = lg(x)-1;
864 1470 : return gen_bkeval(Q,degpol(Q),x,1,(void*)&n,&RgM_algebra,&_RgM_cmul);
865 : }
866 :
867 : GEN
868 4492894 : RgC_Rg_div(GEN x, GEN y)
869 19313119 : { pari_APPLY_type(t_COL, gdiv(gel(x,i),y)) }
870 :
871 : GEN
872 12048979 : RgC_Rg_mul(GEN x, GEN y)
873 76275651 : { pari_APPLY_type(t_COL, gmul(gel(x,i),y)) }
874 :
875 : GEN
876 38990 : RgV_Rg_mul(GEN x, GEN y)
877 1465968 : { pari_APPLY_type(t_VEC, gmul(gel(x,i),y)) }
878 :
879 : GEN
880 606163 : RgM_Rg_div(GEN X, GEN c) {
881 606163 : long i, j, h, l = lg(X);
882 606163 : GEN A = cgetg(l, t_MAT);
883 606163 : if (l == 1) return A;
884 606163 : h = lgcols(X);
885 2808481 : for (j=1; j<l; j++)
886 : {
887 2202318 : GEN a = cgetg(h, t_COL), x = gel(X, j);
888 18099305 : for (i = 1; i < h; i++) gel(a,i) = gdiv(gel(x,i), c);
889 2202318 : gel(A,j) = a;
890 : }
891 606163 : return A;
892 : }
893 : GEN
894 1377899 : RgM_Rg_mul(GEN X, GEN c) {
895 1377899 : long i, j, h, l = lg(X);
896 1377899 : GEN A = cgetg(l, t_MAT);
897 1377899 : if (l == 1) return A;
898 1377899 : h = lgcols(X);
899 4471636 : for (j=1; j<l; j++)
900 : {
901 3093737 : GEN a = cgetg(h, t_COL), x = gel(X, j);
902 13718936 : for (i = 1; i < h; i++) gel(a,i) = gmul(gel(x,i), c);
903 3093737 : gel(A,j) = a;
904 : }
905 1377899 : 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 13435281 : fill_scalmat(GEN y, GEN t, long n)
916 : {
917 : long i;
918 51949593 : for (i = 1; i <= n; i++)
919 : {
920 38514312 : gel(y,i) = zerocol(n);
921 38514312 : gcoeff(y,i,i) = t;
922 : }
923 13435281 : }
924 :
925 : GEN
926 873816 : scalarmat(GEN x, long n) {
927 873816 : GEN y = cgetg(n+1, t_MAT);
928 873816 : if (!n) return y;
929 873816 : fill_scalmat(y, gcopy(x), n); return y;
930 : }
931 : GEN
932 3376953 : scalarmat_shallow(GEN x, long n) {
933 3376953 : GEN y = cgetg(n+1, t_MAT);
934 3376953 : fill_scalmat(y, x, n); return y;
935 : }
936 : GEN
937 252 : scalarmat_s(long x, long n) {
938 252 : GEN y = cgetg(n+1, t_MAT);
939 252 : if (!n) return y;
940 252 : fill_scalmat(y, stoi(x), n); return y;
941 : }
942 : GEN
943 9184267 : matid(long n) {
944 : GEN y;
945 9184267 : if (n < 0) pari_err_DOMAIN("matid", "size", "<", gen_0, stoi(n));
946 9184260 : y = cgetg(n+1, t_MAT);
947 9184260 : fill_scalmat(y, gen_1, n); return y;
948 : }
949 :
950 : INLINE GEN
951 2285907 : scalarcol_i(GEN x, long n, long c)
952 : {
953 : long i;
954 2285907 : GEN y = cgetg(n+1,t_COL);
955 2285907 : if (!n) return y;
956 2285907 : gel(y,1) = c? gcopy(x): x;
957 6832861 : for (i=2; i<=n; i++) gel(y,i) = gen_0;
958 2285907 : return y;
959 : }
960 :
961 : GEN
962 788982 : scalarcol(GEN x, long n) { return scalarcol_i(x,n,1); }
963 :
964 : GEN
965 1496925 : scalarcol_shallow(GEN x, long n) { return scalarcol_i(x,n,0); }
966 :
967 : int
968 18417 : RgM_isscalar(GEN x, GEN s)
969 : {
970 18417 : long i, j, lx = lg(x);
971 :
972 18417 : if (lx == 1) return 1;
973 18417 : if (lx != lgcols(x)) return 0;
974 18417 : if (!s) s = gcoeff(x,1,1);
975 :
976 52228 : for (j=1; j<lx; j++)
977 : {
978 42137 : GEN c = gel(x,j);
979 156536 : for (i=1; i<j; )
980 121311 : if (!gequal0(gel(c,i++))) return 0;
981 : /* i = j */
982 35225 : if (!gequal(gel(c,i++),s)) return 0;
983 178461 : for ( ; i<lx; )
984 144650 : if (!gequal0(gel(c,i++))) return 0;
985 : }
986 10091 : return 1;
987 : }
988 :
989 : int
990 3664 : RgM_isidentity(GEN x)
991 : {
992 3664 : long i,j, lx = lg(x);
993 :
994 3664 : if (lx == 1) return 1;
995 3664 : if (lx != lgcols(x)) return 0;
996 8805 : for (j=1; j<lx; j++)
997 : {
998 8525 : GEN c = gel(x,j);
999 40066 : for (i=1; i<j; )
1000 33380 : if (!gequal0(gel(c,i++))) return 0;
1001 : /* i = j */
1002 6686 : if (!gequal1(gel(c,i++))) return 0;
1003 49771 : for ( ; i<lx; )
1004 44630 : if (!gequal0(gel(c,i++))) return 0;
1005 : }
1006 280 : return 1;
1007 : }
1008 :
1009 : long
1010 644 : RgC_is_ei(GEN x)
1011 : {
1012 644 : long i, j = 0, l = lg(x);
1013 3472 : for (i = 1; i < l; i++)
1014 : {
1015 2828 : GEN c = gel(x,i);
1016 2828 : if (gequal0(c)) continue;
1017 644 : if (!gequal1(c) || j) return 0;
1018 644 : j = i;
1019 : }
1020 644 : return j;
1021 : }
1022 :
1023 : int
1024 336 : RgM_isdiagonal(GEN x)
1025 : {
1026 336 : long i,j, lx = lg(x);
1027 336 : if (lx == 1) return 1;
1028 336 : if (lx != lgcols(x)) return 0;
1029 :
1030 3220 : for (j=1; j<lx; j++)
1031 : {
1032 2891 : GEN c = gel(x,j);
1033 18452 : for (i=1; i<j; i++)
1034 15561 : if (!gequal0(gel(c,i))) return 0;
1035 18452 : for (i++; i<lx; i++)
1036 15568 : if (!gequal0(gel(c,i))) return 0;
1037 : }
1038 329 : return 1;
1039 : }
1040 : int
1041 315 : isdiagonal(GEN x) { return (typ(x)==t_MAT) && RgM_isdiagonal(x); }
1042 :
1043 : GEN
1044 22855 : RgM_det_triangular(GEN mat)
1045 : {
1046 22855 : long i,l = lg(mat);
1047 : pari_sp av;
1048 : GEN s;
1049 :
1050 22855 : if (l<3) return l<2? gen_1: gcopy(gcoeff(mat,1,1));
1051 20944 : av = avma; s = gcoeff(mat,1,1);
1052 122360 : for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1053 20944 : return av==avma? gcopy(s): gc_upto(av,s);
1054 : }
1055 :
1056 : GEN
1057 484910 : RgV_kill0(GEN v)
1058 : {
1059 : long i, l;
1060 484910 : GEN w = cgetg_copy(v, &l);
1061 133689059 : for (i = 1; i < l; i++)
1062 : {
1063 133204149 : GEN a = gel(v,i);
1064 133204149 : gel(w,i) = gequal0(a) ? NULL: a;
1065 : }
1066 484910 : return w;
1067 : }
|
