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 1232389 : RgM_is_ZM(GEN x)
20 : {
21 1232389 : long i, j, h, l = lg(x);
22 1232389 : if (l == 1) return 1;
23 1231724 : h = lgcols(x);
24 1231713 : if (h == 1) return 1;
25 4668246 : for (j = l-1; j > 0; j--)
26 22317295 : for (i = h-1; i > 0; i--)
27 18880601 : if (typ(gcoeff(x,i,j)) != t_INT) return 0;
28 1074300 : 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 8295455 : RgMrow_zc_mul_i(GEN x, GEN y, long c, long i)
64 : {
65 8295455 : pari_sp av = avma;
66 8295455 : GEN s = NULL;
67 : long j;
68 829682939 : for (j=1; j<c; j++)
69 : {
70 821390616 : long t = y[j];
71 821390616 : if (!t) continue;
72 28627276 : if (!s) { s = gmulgs(gcoeff(x,i,j),t); continue; }
73 20550495 : switch(t)
74 : {
75 15034372 : case 1: s = gadd(s, gcoeff(x,i,j)); break;
76 239931 : case -1: s = gsub(s, gcoeff(x,i,j)); break;
77 5276192 : default: s = gadd(s, gmulgs(gcoeff(x,i,j), t)); break;
78 : }
79 : }
80 8292323 : return s? gerepileupto(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 239081 : RgM_zc_mul_i(GEN x, GEN y, long c, long l)
87 : {
88 239081 : GEN z = cgetg(l,t_COL);
89 : long i;
90 8464481 : for (i = 1; i < l; i++) gel(z,i) = RgMrow_zc_mul_i(x,y,c,i);
91 239068 : 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 56749 : RgM_zm_mul(GEN x, GEN y)
98 : {
99 56749 : long j, c, l = lg(x), ly = lg(y);
100 56749 : GEN z = cgetg(ly, t_MAT);
101 56749 : if (l == 1) return z;
102 56749 : c = lgcols(x);
103 223849 : for (j = 1; j < ly; j++) gel(z,j) = RgM_zc_mul_i(x, gel(y,j), l,c);
104 56749 : return z;
105 : }
106 :
107 : /* x[i,]*y, l = lg(y) > 1 */
108 : static GEN
109 32324143 : RgMrow_ZC_mul_i(GEN x, GEN y, long i, long l)
110 : {
111 32324143 : pari_sp av = avma;
112 32324143 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
113 : long j;
114 5022131365 : for (j=2; j<l; j++)
115 4990113637 : if (signe(gel(y,j))) t = gadd(t, gmul(gcoeff(x,i,j), gel(y,j)));
116 32017728 : return gerepileupto(av,t);
117 : }
118 :
119 : /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
120 : static GEN
121 1928402 : RgM_ZC_mul_i(GEN x, GEN y, long lx, long l)
122 : {
123 1928402 : GEN z = cgetg(l,t_COL);
124 : long i;
125 34255808 : for (i=1; i<l; i++) gel(z,i) = RgMrow_ZC_mul_i(x,y,i,lx);
126 1925319 : return z;
127 : }
128 :
129 : GEN
130 1928430 : RgM_ZM_mul_worker(GEN y, GEN x)
131 1928430 : { return RgM_ZC_mul_i(x, y, lg(x), lgcols(x)); }
132 :
133 : /* mostly useful when y is sparse */
134 : GEN
135 354834 : RgM_ZM_mul(GEN x, GEN y)
136 : {
137 354834 : pari_sp av = avma;
138 : GEN worker;
139 354834 : if (lg(x) == 1) return cgetg(lg(y), t_MAT);
140 354834 : worker = snm_closure(is_entry("_RgM_ZM_mul_worker"),mkvec(x));
141 354836 : return gerepileupto(av, gen_parapply(worker, y));
142 : }
143 :
144 : static GEN
145 210197 : RgV_zc_mul_i(GEN x, GEN y, long l)
146 : {
147 : long i;
148 210197 : GEN z = gen_0;
149 210197 : pari_sp av = avma;
150 8215426 : for (i = 1; i < l; i++) z = gadd(z, gmulgs(gel(x,i), y[i]));
151 210195 : return gerepileupto(av, z);
152 : }
153 : GEN
154 83384 : RgV_zc_mul(GEN x, GEN y) { return RgV_zc_mul_i(x, y, lg(x)); }
155 :
156 : GEN
157 34274 : RgV_zm_mul(GEN x, GEN y)
158 : {
159 34274 : long j, l = lg(x), ly = lg(y);
160 34274 : GEN z = cgetg(ly, t_VEC);
161 161086 : for (j = 1; j < ly; j++) gel(z,j) = RgV_zc_mul_i(x, gel(y,j), l);
162 34273 : return z;
163 : }
164 :
165 : /* scalar product x.x */
166 : GEN
167 9446143 : RgV_dotsquare(GEN x)
168 : {
169 9446143 : long i, lx = lg(x);
170 9446143 : pari_sp av = avma;
171 : GEN z;
172 9446143 : if (lx == 1) return gen_0;
173 9446143 : z = gsqr(gel(x,1));
174 21161690 : for (i=2; i<lx; i++)
175 : {
176 11715787 : z = gadd(z, gsqr(gel(x,i)));
177 11715721 : if (gc_needed(av,3))
178 : {
179 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotsquare, i = %ld",i);
180 0 : z = gerepileupto(av, z);
181 : }
182 : }
183 9445903 : return gerepileupto(av,z);
184 : }
185 :
186 : /* scalar product x.y, lx = lg(x) = lg(y) */
187 : static GEN
188 13607066 : RgV_dotproduct_i(GEN x, GEN y, long lx)
189 : {
190 13607066 : pari_sp av = avma;
191 : long i;
192 : GEN z;
193 13607066 : if (lx == 1) return gen_0;
194 13606667 : z = gmul(gel(x,1),gel(y,1));
195 218996000 : for (i=2; i<lx; i++)
196 : {
197 205390772 : z = gadd(z, gmul(gel(x,i), gel(y,i)));
198 205391097 : if (gc_needed(av,3))
199 : {
200 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotproduct, i = %ld",i);
201 0 : z = gerepileupto(av, z);
202 : }
203 : }
204 13605228 : return gerepileupto(av,z);
205 : }
206 : GEN
207 1012433 : RgV_dotproduct(GEN x,GEN y)
208 : {
209 1012433 : if (x == y) return RgV_dotsquare(x);
210 1012433 : return RgV_dotproduct_i(x, y, lg(x));
211 : }
212 : /* v[1] + ... + v[lg(v)-1] */
213 : GEN
214 3032796 : RgV_sum(GEN v)
215 : {
216 : GEN p;
217 3032796 : long i, l = lg(v);
218 3032796 : if (l == 1) return gen_0;
219 7381196 : p = gel(v,1); for (i=2; i<l; i++) p = gadd(p, gel(v,i));
220 3032761 : 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) = gerepileupto(av, s);
258 : }
259 216 : return v;
260 : }
261 :
262 : static GEN
263 3335973 : _gmul(void *data, GEN x, GEN y)
264 3335973 : { (void)data; return gmul(x,y); }
265 :
266 : GEN
267 452689 : RgV_prod(GEN x)
268 : {
269 452689 : 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 7580 : RgM_Rg_add(GEN x, GEN y)
276 : {
277 7580 : long l = lg(x), i, j;
278 7580 : GEN z = cgetg(l,t_MAT);
279 :
280 7580 : if (l==1) return z;
281 7580 : if (l != lgcols(x)) pari_err_OP("+", x, y);
282 57215 : for (i=1; i<l; i++)
283 : {
284 49635 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
285 49635 : gel(z,i) = zi;
286 2056800 : for (j=1; j<l; j++)
287 2007165 : gel(zi,j) = i==j? gadd(y,gel(xi,j)): gcopy(gel(xi,j));
288 : }
289 7580 : 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 313933 : RgM_Rg_sub_shallow(GEN x, GEN y)
327 : {
328 313933 : long l = lg(x), i, j;
329 313933 : GEN z = cgetg(l,t_MAT);
330 :
331 313933 : if (l==1) return z;
332 313933 : if (l != lgcols(x)) pari_err_OP( "-", x, y);
333 2027277 : for (i=1; i<l; i++)
334 : {
335 1713363 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
336 1713363 : gel(z,i) = zi;
337 27263497 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
338 1713363 : gel(zi,i) = gsub(gel(zi,i), y);
339 : }
340 313914 : return z;
341 : }
342 :
343 : GEN
344 2388162 : RgC_Rg_add(GEN x, GEN y)
345 : {
346 2388162 : long k, lx = lg(x);
347 2388162 : GEN z = cgetg(lx, t_COL);
348 2388166 : if (lx == 1)
349 : {
350 0 : if (isintzero(y)) return z;
351 0 : pari_err_TYPE2("+",x,y);
352 : }
353 2388166 : gel(z,1) = gadd(y,gel(x,1));
354 6739934 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
355 2388155 : return z;
356 : }
357 : GEN
358 23934 : RgC_Rg_sub(GEN x, GEN y)
359 : {
360 23934 : long k, lx = lg(x);
361 23934 : GEN z = cgetg(lx, t_COL);
362 23934 : if (lx == 1)
363 : {
364 0 : if (isintzero(y)) return z;
365 0 : pari_err_TYPE2("-",x,y);
366 : }
367 23934 : gel(z,1) = gsub(gel(x,1), y);
368 68878 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
369 23934 : return z;
370 : }
371 : /* a - x */
372 : GEN
373 32738 : Rg_RgC_sub(GEN a, GEN x)
374 : {
375 32738 : long k, lx = lg(x);
376 32738 : GEN z = cgetg(lx,t_COL);
377 32738 : if (lx == 1)
378 : {
379 0 : if (isintzero(a)) return z;
380 0 : pari_err_TYPE2("-",a,x);
381 : }
382 32738 : gel(z,1) = gsub(a, gel(x,1));
383 74661 : for (k = 2; k < lx; k++) gel(z,k) = gneg(gel(x,k));
384 32738 : return z;
385 : }
386 :
387 : static GEN
388 18683539 : RgC_add_i(GEN x, GEN y, long lx)
389 : {
390 18683539 : GEN A = cgetg(lx, t_COL);
391 : long i;
392 142288898 : for (i=1; i<lx; i++) gel(A,i) = gadd(gel(x,i), gel(y,i));
393 18683216 : return A;
394 : }
395 : GEN
396 15416861 : RgC_add(GEN x, GEN y) { return RgC_add_i(x, y, lg(x)); }
397 : GEN
398 2547021 : RgV_add(GEN x, GEN y)
399 10436267 : { pari_APPLY_type(t_VEC, gadd(gel(x,i), gel(y,i))) }
400 :
401 : static GEN
402 5748634 : RgC_sub_i(GEN x, GEN y, long lx)
403 : {
404 : long i;
405 5748634 : GEN A = cgetg(lx, t_COL);
406 634253643 : for (i=1; i<lx; i++) gel(A,i) = gsub(gel(x,i), gel(y,i));
407 5747596 : return A;
408 : }
409 : GEN
410 4844679 : RgC_sub(GEN x, GEN y) { return RgC_sub_i(x, y, lg(x)); }
411 : GEN
412 356557 : RgV_sub(GEN x, GEN y)
413 3032182 : { pari_APPLY_type(t_VEC, gsub(gel(x,i), gel(y,i))) }
414 :
415 : GEN
416 690830 : RgM_add(GEN x, GEN y)
417 : {
418 690830 : long lx = lg(x), l, j;
419 : GEN z;
420 690830 : if (lx == 1) return cgetg(1, t_MAT);
421 690830 : z = cgetg(lx, t_MAT); l = lgcols(x);
422 3957506 : for (j = 1; j < lx; j++) gel(z,j) = RgC_add_i(gel(x,j), gel(y,j), l);
423 690830 : return z;
424 : }
425 : GEN
426 172513 : RgM_sub(GEN x, GEN y)
427 : {
428 172513 : long lx = lg(x), l, j;
429 : GEN z;
430 172513 : if (lx == 1) return cgetg(1, t_MAT);
431 172513 : z = cgetg(lx, t_MAT); l = lgcols(x);
432 1076500 : for (j = 1; j < lx; j++) gel(z,j) = RgC_sub_i(gel(x,j), gel(y,j), l);
433 172511 : return z;
434 : }
435 :
436 : static GEN
437 3919547 : RgC_neg_i(GEN x, long lx)
438 : {
439 : long i;
440 3919547 : GEN y = cgetg(lx, t_COL);
441 29923790 : for (i=1; i<lx; i++) gel(y,i) = gneg(gel(x,i));
442 3919544 : return y;
443 : }
444 : GEN
445 932034 : RgC_neg(GEN x) { return RgC_neg_i(x, lg(x)); }
446 : GEN
447 4515 : RgV_neg(GEN x)
448 60473 : { pari_APPLY_type(t_VEC, gneg(gel(x,i))) }
449 : GEN
450 542170 : RgM_neg(GEN x)
451 : {
452 542170 : long i, hx, lx = lg(x);
453 542170 : GEN y = cgetg(lx, t_MAT);
454 542170 : if (lx == 1) return y;
455 542163 : hx = lgcols(x);
456 3529676 : for (i=1; i<lx; i++) gel(y,i) = RgC_neg_i(gel(x,i), hx);
457 542163 : return y;
458 : }
459 :
460 : GEN
461 4725190 : RgV_RgC_mul(GEN x, GEN y)
462 : {
463 4725190 : long lx = lg(x);
464 4725190 : if (lx != lg(y)) pari_err_OP("operation 'RgV_RgC_mul'", x, y);
465 4725106 : 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 115809480 : RgMrow_RgC_mul_i(GEN x, GEN y, long i, long l)
494 : {
495 115809480 : pari_sp av = avma;
496 115809480 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
497 : long j;
498 985385186 : for (j=2; j<l; j++)
499 : {
500 869585737 : GEN c = gcoeff(x,i,j);
501 869585737 : if (!isintzero(c)) t = gadd(t, gmul(c, gel(y,j)));
502 : }
503 115799449 : return gerepileupto(av,t);
504 : }
505 : GEN
506 4655 : RgMrow_RgC_mul(GEN x, GEN y, long i)
507 4655 : { 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 gerepileupto(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 gerepileupto(av, FqC_to_mod(b, T, p));
533 : }
534 :
535 : #define code(t1,t2) ((t1 << 6) | t2)
536 : static GEN
537 31171990 : RgM_RgC_mul_fast(GEN x, GEN y)
538 : {
539 : GEN p, pol;
540 : long pa;
541 31171990 : long t = RgM_RgC_type(x,y, &p,&pol,&pa);
542 31172252 : switch(t)
543 : {
544 8830155 : case t_INT: return ZM_ZC_mul(x,y);
545 126535 : case t_FRAC: return QM_QC_mul(x,y);
546 91 : case t_FFELT: return FFM_FFC_mul(x, y, pol);
547 28 : case t_INTMOD: return RgM_RgC_mul_FpM(x, y, p);
548 24 : case code(t_POLMOD, t_INTMOD):
549 24 : return RgM_RgC_mul_FqM(x, y, pol, p);
550 22215423 : default: return NULL;
551 : }
552 : }
553 : #undef code
554 :
555 : /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
556 : static GEN
557 34004544 : RgM_RgC_mul_i(GEN x, GEN y, long lx, long l)
558 : {
559 34004544 : GEN z = cgetg(l,t_COL);
560 : long i;
561 149804991 : for (i=1; i<l; i++) gel(z,i) = RgMrow_RgC_mul_i(x,y,i,lx);
562 34000721 : return z;
563 : }
564 :
565 : GEN
566 31171998 : RgM_RgC_mul(GEN x, GEN y)
567 : {
568 31171998 : long lx = lg(x);
569 : GEN z;
570 31171998 : if (lx != lg(y)) pari_err_OP("operation 'RgM_RgC_mul'", x,y);
571 31171999 : if (lx == 1) return cgetg(1,t_COL);
572 31171999 : z = RgM_RgC_mul_fast(x, y);
573 31172240 : if (z) return z;
574 22215404 : return RgM_RgC_mul_i(x, y, lx, lgcols(x));
575 : }
576 :
577 : GEN
578 295684 : RgV_RgM_mul(GEN x, GEN y)
579 : {
580 295684 : long i, lx, ly = lg(y);
581 : GEN z;
582 295684 : if (ly == 1) return cgetg(1,t_VEC);
583 295677 : lx = lg(x);
584 295677 : if (lx != lgcols(y)) pari_err_OP("operation 'RgV_RgM_mul'", x,y);
585 295670 : z = cgetg(ly, t_VEC);
586 3423839 : for (i=1; i<ly; i++) gel(z,i) = RgV_dotproduct_i(x, gel(y,i), lx);
587 295560 : return z;
588 : }
589 :
590 : static GEN
591 57148 : RgM_mul_FpM_i(GEN x, GEN y, GEN p)
592 : {
593 57148 : if (lgefint(p) == 3)
594 : {
595 57076 : ulong pp = uel(p, 2);
596 57076 : if (pp==2) return F2m_to_mod(F2m_mul(RgM_to_F2m(x), RgM_to_F2m(y)));
597 57027 : if (pp==3) return F3m_to_mod(F3m_mul(RgM_to_F3m(x), RgM_to_F3m(y)));
598 56873 : return Flm_to_mod(Flm_mul(RgM_to_Flm(x, pp), RgM_to_Flm(y, pp), pp), pp);
599 : }
600 72 : return FpM_to_mod(FpM_mul(RgM_to_FpM(x,p), RgM_to_FpM(y,p), p), p);
601 : }
602 : static GEN
603 57148 : RgM_mul_FpM(GEN x, GEN y, GEN p)
604 57148 : { pari_sp av = avma; return gerepileupto(av, RgM_mul_FpM_i(x, y, p)); }
605 :
606 : static GEN
607 26005 : RgM_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
608 : {
609 26005 : pari_sp av = avma;
610 26005 : GEN b, T = RgX_to_FpX(pol, p);
611 26005 : if (signe(T) == 0) pari_err_OP("*", x, y);
612 26005 : b = FqM_mul(RgM_to_FqM(x, T, p), RgM_to_FqM(y, T, p), T, p);
613 26005 : return gerepileupto(av, FqM_to_mod(b, T, p));
614 : }
615 :
616 : static GEN
617 12446 : RgM_liftred(GEN x, GEN T)
618 12446 : { return RgXQM_red(liftpol_shallow(x), T); }
619 :
620 : static GEN
621 1449 : RgM_mul_ZXQM(GEN x, GEN y, GEN T)
622 : {
623 1449 : pari_sp av = avma;
624 1449 : GEN b = ZXQM_mul(RgM_liftred(x,T), RgM_liftred(y, T), T);
625 1449 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
626 : }
627 :
628 : static GEN
629 133 : RgM_sqr_ZXQM(GEN x, GEN T)
630 : {
631 133 : pari_sp av = avma;
632 133 : GEN b = ZXQM_sqr(RgM_liftred(x, T), T);
633 133 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
634 : }
635 :
636 : static GEN
637 4704 : RgM_mul_QXQM(GEN x, GEN y, GEN T)
638 : {
639 4704 : pari_sp av = avma;
640 4704 : GEN b = QXQM_mul(RgM_liftred(x, T), RgM_liftred(y, T), T);
641 4704 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
642 : }
643 :
644 : static GEN
645 7 : RgM_sqr_QXQM(GEN x, GEN T)
646 : {
647 7 : pari_sp av = avma;
648 7 : GEN b = QXQM_sqr(RgM_liftred(x, T), T);
649 7 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
650 : }
651 :
652 : INLINE int
653 4753 : RgX_is_monic_ZX(GEN pol)
654 4753 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
655 :
656 : #define code(t1,t2) ((t1 << 6) | t2)
657 : static GEN
658 7571993 : RgM_mul_fast(GEN x, GEN y)
659 : {
660 : GEN p, pol;
661 : long pa;
662 7571993 : long t = RgM_type2(x,y, &p,&pol,&pa);
663 7572060 : switch(t)
664 : {
665 2689583 : case t_INT: return ZM_mul(x,y);
666 117138 : case t_FRAC: return QM_mul(x,y);
667 4403 : case t_FFELT: return FFM_mul(x, y, pol);
668 57085 : case t_INTMOD: return RgM_mul_FpM(x, y, p);
669 1456 : case code(t_POLMOD, t_INT):
670 1456 : return ZX_is_monic(pol)? RgM_mul_ZXQM(x, y, pol): NULL;
671 4725 : case code(t_POLMOD, t_FRAC):
672 4725 : return RgX_is_monic_ZX(pol)? RgM_mul_QXQM(x, y, pol): NULL;
673 26004 : case code(t_POLMOD, t_INTMOD):
674 26004 : return RgM_mul_FqM(x, y, pol, p);
675 4671666 : default: return NULL;
676 : }
677 : }
678 :
679 : static GEN
680 1372 : RgM_sqr_fast(GEN x)
681 : {
682 : GEN p, pol;
683 : long pa;
684 1372 : long t = RgM_type(x, &p,&pol,&pa);
685 1372 : switch(t)
686 : {
687 231 : case t_INT: return ZM_sqr(x);
688 700 : case t_FRAC: return QM_sqr(x);
689 112 : case t_FFELT: return FFM_mul(x, x, pol);
690 63 : case t_INTMOD: return RgM_mul_FpM(x, x, p);
691 140 : case code(t_POLMOD, t_INT):
692 140 : return ZX_is_monic(pol)? RgM_sqr_ZXQM(x, pol): NULL;
693 28 : case code(t_POLMOD, t_FRAC):
694 28 : return RgX_is_monic_ZX(pol)? RgM_sqr_QXQM(x, pol): NULL;
695 0 : case code(t_POLMOD, t_INTMOD):
696 0 : return RgM_mul_FqM(x, x, pol, p);
697 98 : default: return NULL;
698 : }
699 : }
700 :
701 : #undef code
702 :
703 : /* lx, ly > 1 */
704 : static GEN
705 4671694 : RgM_mul_i(GEN x, GEN y, long lx, long ly)
706 : {
707 4671694 : GEN z = cgetg(ly, t_MAT);
708 4671689 : long j, l = lgcols(x);
709 16460644 : for (j = 1; j < ly; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(y,j), lx, l);
710 4671658 : return z;
711 : }
712 : GEN
713 7587771 : RgM_mul(GEN x, GEN y)
714 : {
715 7587771 : long lx, ly = lg(y);
716 : GEN z;
717 7587771 : if (ly == 1) return cgetg(1,t_MAT);
718 7572021 : lx = lg(x);
719 7572021 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
720 7572021 : if (lx == 1) return zeromat(0,ly-1);
721 7571993 : z = RgM_mul_fast(x, y);
722 7572056 : if (z) return z;
723 4671689 : return RgM_mul_i(x, y, lx, ly);
724 : }
725 :
726 : GEN
727 1407 : RgM_sqr(GEN x)
728 : {
729 1407 : long j, lx = lg(x);
730 : GEN z;
731 1407 : if (lx == 1) return cgetg(1, t_MAT);
732 1372 : if (lx != lgcols(x)) pari_err_OP("operation 'RgM_mul'", x,x);
733 1372 : z = RgM_sqr_fast(x);
734 1372 : if (z) return z;
735 126 : z = cgetg(lx, t_MAT);
736 497 : for (j=1; j<lx; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(x,j), lx, lx);
737 126 : return z;
738 : }
739 :
740 : /* assume result is symmetric */
741 : GEN
742 0 : RgM_multosym(GEN x, GEN y)
743 : {
744 0 : long j, lx, ly = lg(y);
745 : GEN M;
746 0 : if (ly == 1) return cgetg(1,t_MAT);
747 0 : lx = lg(x);
748 0 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_multosym'", x,y);
749 0 : if (lx == 1) return cgetg(1,t_MAT);
750 0 : if (ly != lgcols(x)) pari_err_OP("operation 'RgM_multosym'", x,y);
751 0 : M = cgetg(ly, t_MAT);
752 0 : for (j=1; j<ly; j++)
753 : {
754 0 : GEN z = cgetg(ly,t_COL), yj = gel(y,j);
755 : long i;
756 0 : for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
757 0 : for (i=j; i<ly; i++)gel(z,i) = RgMrow_RgC_mul_i(x,yj,i,lx);
758 0 : gel(M,j) = z;
759 : }
760 0 : return M;
761 : }
762 : /* x~ * y, assuming result is symmetric */
763 : GEN
764 7663 : RgM_transmultosym(GEN x, GEN y)
765 : {
766 7663 : long i, j, l, ly = lg(y);
767 : GEN M;
768 7663 : if (ly == 1) return cgetg(1,t_MAT);
769 7663 : if (lg(x) != ly) pari_err_OP("operation 'RgM_transmultosym'", x,y);
770 7663 : l = lgcols(y);
771 7663 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmultosym'", x,y);
772 7663 : M = cgetg(ly, t_MAT);
773 30820 : for (i=1; i<ly; i++)
774 : {
775 23157 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
776 23157 : gel(M,i) = c;
777 47490 : for (j=1; j<i; j++)
778 24333 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(y,j),l);
779 23157 : gel(c,i) = RgV_dotproduct_i(xi,gel(y,i),l);
780 : }
781 7663 : return M;
782 : }
783 : /* x~ * y */
784 : GEN
785 0 : RgM_transmul(GEN x, GEN y)
786 : {
787 0 : long i, j, l, lx, ly = lg(y);
788 : GEN M;
789 0 : if (ly == 1) return cgetg(1,t_MAT);
790 0 : lx = lg(x);
791 0 : l = lgcols(y);
792 0 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmul'", x,y);
793 0 : M = cgetg(ly, t_MAT);
794 0 : for (i=1; i<ly; i++)
795 : {
796 0 : GEN yi = gel(y,i), c = cgetg(lx,t_COL);
797 0 : gel(M,i) = c;
798 0 : for (j=1; j<lx; j++) gel(c,j) = RgV_dotproduct_i(yi,gel(x,j),l);
799 : }
800 0 : return M;
801 : }
802 :
803 : GEN
804 4693817 : gram_matrix(GEN x)
805 : {
806 4693817 : long i,j, l, lx = lg(x);
807 : GEN M;
808 4693817 : if (!is_matvec_t(typ(x))) pari_err_TYPE("gram",x);
809 4693817 : if (lx == 1) return cgetg(1,t_MAT);
810 4693803 : l = lgcols(x);
811 4693803 : M = cgetg(lx,t_MAT);
812 14081347 : for (i=1; i<lx; i++)
813 : {
814 9387557 : GEN xi = gel(x,i), c = cgetg(lx,t_COL);
815 9387549 : gel(M,i) = c;
816 14081320 : for (j=1; j<i; j++)
817 4693770 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(x,j),l);
818 9387550 : gel(c,i) = RgV_dotsquare(xi);
819 : }
820 4693790 : return M;
821 : }
822 :
823 : static GEN
824 3801 : _RgM_add(void *E, GEN x, GEN y) { (void)E; return RgM_add(x, y); }
825 :
826 : static GEN
827 0 : _RgM_sub(void *E, GEN x, GEN y) { (void)E; return RgM_sub(x, y); }
828 :
829 : static GEN
830 5684 : _RgM_cmul(void *E, GEN P, long a, GEN x) { (void)E; return RgM_Rg_mul(x,gel(P,a+2)); }
831 :
832 : static GEN
833 238 : _RgM_sqr(void *E, GEN x) { (void) E; return RgM_sqr(x); }
834 :
835 : static GEN
836 749 : _RgM_mul(void *E, GEN x, GEN y) { (void) E; return RgM_mul(x, y); }
837 :
838 : static GEN
839 4046 : _RgM_one(void *E) { long *n = (long*) E; return matid(*n); }
840 :
841 : static GEN
842 0 : _RgM_zero(void *E) { long *n = (long*) E; return zeromat(*n,*n); }
843 :
844 : static GEN
845 2828 : _RgM_red(void *E, GEN x) { (void)E; return x; }
846 :
847 : static struct bb_algebra RgM_algebra = { _RgM_red, _RgM_add, _RgM_sub,
848 : _RgM_mul, _RgM_sqr, _RgM_one, _RgM_zero };
849 :
850 : /* generates the list of powers of x of degree 0,1,2,...,l*/
851 : GEN
852 168 : RgM_powers(GEN x, long l)
853 : {
854 168 : long n = lg(x)-1;
855 168 : return gen_powers(x,l,1,(void *) &n, &_RgM_sqr, &_RgM_mul, &_RgM_one);
856 : }
857 :
858 : GEN
859 490 : RgX_RgMV_eval(GEN Q, GEN x)
860 : {
861 490 : long n = lg(x)>1 ? lg(gel(x,1))-1:0;
862 490 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&n,&RgM_algebra,&_RgM_cmul);
863 : }
864 :
865 : GEN
866 1393 : RgX_RgM_eval(GEN Q, GEN x)
867 : {
868 1393 : long n = lg(x)-1;
869 1393 : return gen_bkeval(Q,degpol(Q),x,1,(void*)&n,&RgM_algebra,&_RgM_cmul);
870 : }
871 :
872 : GEN
873 4406687 : RgC_Rg_div(GEN x, GEN y)
874 19018979 : { pari_APPLY_type(t_COL, gdiv(gel(x,i),y)) }
875 :
876 : GEN
877 11842681 : RgC_Rg_mul(GEN x, GEN y)
878 75457831 : { pari_APPLY_type(t_COL, gmul(gel(x,i),y)) }
879 :
880 : GEN
881 31157 : RgV_Rg_mul(GEN x, GEN y)
882 1200766 : { pari_APPLY_type(t_VEC, gmul(gel(x,i),y)) }
883 :
884 : GEN
885 603281 : RgM_Rg_div(GEN X, GEN c) {
886 603281 : long i, j, h, l = lg(X);
887 603281 : GEN A = cgetg(l, t_MAT);
888 603281 : if (l == 1) return A;
889 603281 : h = lgcols(X);
890 2781781 : for (j=1; j<l; j++)
891 : {
892 2178529 : GEN a = cgetg(h, t_COL), x = gel(X, j);
893 17366494 : for (i = 1; i < h; i++) gel(a,i) = gdiv(gel(x,i), c);
894 2178500 : gel(A,j) = a;
895 : }
896 603252 : return A;
897 : }
898 : GEN
899 218726 : RgM_Rg_mul(GEN X, GEN c) {
900 218726 : long i, j, h, l = lg(X);
901 218726 : GEN A = cgetg(l, t_MAT);
902 218725 : if (l == 1) return A;
903 218725 : h = lgcols(X);
904 997152 : for (j=1; j<l; j++)
905 : {
906 778440 : GEN a = cgetg(h, t_COL), x = gel(X, j);
907 6866809 : for (i = 1; i < h; i++) gel(a,i) = gmul(gel(x,i), c);
908 778426 : gel(A,j) = a;
909 : }
910 218712 : return A;
911 : }
912 :
913 : /********************************************************************/
914 : /* */
915 : /* SCALAR TO MATRIX/VECTOR */
916 : /* */
917 : /********************************************************************/
918 : /* fill the square nxn matrix equal to t*Id */
919 : static void
920 12703064 : fill_scalmat(GEN y, GEN t, long n)
921 : {
922 : long i;
923 47869083 : for (i = 1; i <= n; i++)
924 : {
925 35166312 : gel(y,i) = zerocol(n);
926 35166019 : gcoeff(y,i,i) = t;
927 : }
928 12702771 : }
929 :
930 : GEN
931 867575 : scalarmat(GEN x, long n) {
932 867575 : GEN y = cgetg(n+1, t_MAT);
933 867551 : if (!n) return y;
934 867551 : fill_scalmat(y, gcopy(x), n); return y;
935 : }
936 : GEN
937 3153978 : scalarmat_shallow(GEN x, long n) {
938 3153978 : GEN y = cgetg(n+1, t_MAT);
939 3154019 : fill_scalmat(y, x, n); return y;
940 : }
941 : GEN
942 217 : scalarmat_s(long x, long n) {
943 217 : GEN y = cgetg(n+1, t_MAT);
944 217 : if (!n) return y;
945 217 : fill_scalmat(y, stoi(x), n); return y;
946 : }
947 : GEN
948 8681965 : matid(long n) {
949 : GEN y;
950 8681965 : if (n < 0) pari_err_DOMAIN("matid", "size", "<", gen_0, stoi(n));
951 8681958 : y = cgetg(n+1, t_MAT);
952 8682011 : fill_scalmat(y, gen_1, n); return y;
953 : }
954 :
955 : INLINE GEN
956 2257695 : scalarcol_i(GEN x, long n, long c)
957 : {
958 : long i;
959 2257695 : GEN y = cgetg(n+1,t_COL);
960 2257682 : if (!n) return y;
961 2257682 : gel(y,1) = c? gcopy(x): x;
962 6737746 : for (i=2; i<=n; i++) gel(y,i) = gen_0;
963 2257683 : return y;
964 : }
965 :
966 : GEN
967 779540 : scalarcol(GEN x, long n) { return scalarcol_i(x,n,1); }
968 :
969 : GEN
970 1478155 : scalarcol_shallow(GEN x, long n) { return scalarcol_i(x,n,0); }
971 :
972 : int
973 17138 : RgM_isscalar(GEN x, GEN s)
974 : {
975 17138 : long i, j, lx = lg(x);
976 :
977 17138 : if (lx == 1) return 1;
978 17138 : if (lx != lgcols(x)) return 0;
979 17138 : if (!s) s = gcoeff(x,1,1);
980 :
981 48578 : for (j=1; j<lx; j++)
982 : {
983 39607 : GEN c = gel(x,j);
984 152571 : for (i=1; i<j; )
985 119773 : if (!gequal0(gel(c,i++))) return 0;
986 : /* i = j */
987 32798 : if (!gequal(gel(c,i++),s)) return 0;
988 174345 : for ( ; i<lx; )
989 142905 : if (!gequal0(gel(c,i++))) return 0;
990 : }
991 8971 : return 1;
992 : }
993 :
994 : int
995 3576 : RgM_isidentity(GEN x)
996 : {
997 3576 : long i,j, lx = lg(x);
998 :
999 3576 : if (lx == 1) return 1;
1000 3576 : if (lx != lgcols(x)) return 0;
1001 8561 : for (j=1; j<lx; j++)
1002 : {
1003 8316 : GEN c = gel(x,j);
1004 39661 : for (i=1; i<j; )
1005 33178 : if (!gequal0(gel(c,i++))) return 0;
1006 : /* i = j */
1007 6483 : if (!gequal1(gel(c,i++))) return 0;
1008 49268 : for ( ; i<lx; )
1009 44283 : if (!gequal0(gel(c,i++))) return 0;
1010 : }
1011 245 : return 1;
1012 : }
1013 :
1014 : long
1015 581 : RgC_is_ei(GEN x)
1016 : {
1017 581 : long i, j = 0, l = lg(x);
1018 3157 : for (i = 1; i < l; i++)
1019 : {
1020 2576 : GEN c = gel(x,i);
1021 2576 : if (gequal0(c)) continue;
1022 581 : if (!gequal1(c) || j) return 0;
1023 581 : j = i;
1024 : }
1025 581 : return j;
1026 : }
1027 :
1028 : int
1029 336 : RgM_isdiagonal(GEN x)
1030 : {
1031 336 : long i,j, lx = lg(x);
1032 336 : if (lx == 1) return 1;
1033 336 : if (lx != lgcols(x)) return 0;
1034 :
1035 3220 : for (j=1; j<lx; j++)
1036 : {
1037 2891 : GEN c = gel(x,j);
1038 18452 : for (i=1; i<j; i++)
1039 15561 : if (!gequal0(gel(c,i))) return 0;
1040 18452 : for (i++; i<lx; i++)
1041 15568 : if (!gequal0(gel(c,i))) return 0;
1042 : }
1043 329 : return 1;
1044 : }
1045 : int
1046 315 : isdiagonal(GEN x) { return (typ(x)==t_MAT) && RgM_isdiagonal(x); }
1047 :
1048 : GEN
1049 21981 : RgM_det_triangular(GEN mat)
1050 : {
1051 21981 : long i,l = lg(mat);
1052 : pari_sp av;
1053 : GEN s;
1054 :
1055 21981 : if (l<3) return l<2? gen_1: gcopy(gcoeff(mat,1,1));
1056 20343 : av = avma; s = gcoeff(mat,1,1);
1057 121091 : for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1058 20343 : return av==avma? gcopy(s): gerepileupto(av,s);
1059 : }
1060 :
1061 : GEN
1062 478032 : RgV_kill0(GEN v)
1063 : {
1064 : long i, l;
1065 478032 : GEN w = cgetg_copy(v, &l);
1066 129984645 : for (i = 1; i < l; i++)
1067 : {
1068 129507377 : GEN a = gel(v,i);
1069 129507377 : gel(w,i) = gequal0(a) ? NULL: a;
1070 : }
1071 477268 : return w;
1072 : }
|
