Line data Source code
1 : /* Copyright (C) 2012 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 : /********************************************************************/
16 : /** **/
17 : /** LINEAR ALGEBRA **/
18 : /** (third part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_mat
25 :
26 : /*******************************************************************/
27 : /* */
28 : /* SUM */
29 : /* */
30 : /*******************************************************************/
31 :
32 : GEN
33 110127 : vecsum(GEN v)
34 : {
35 110127 : pari_sp av = avma;
36 : long i, l;
37 : GEN p;
38 110127 : if (!is_vec_t(typ(v)))
39 7 : pari_err_TYPE("vecsum", v);
40 110120 : l = lg(v);
41 110120 : if (l == 1) return gen_0;
42 110113 : p = gel(v,1);
43 110113 : if (l == 2) return gcopy(p);
44 188714 : for (i=2; i<l; i++)
45 : {
46 134986 : p = gadd(p, gel(v,i));
47 134986 : if (gc_needed(av, 2))
48 : {
49 0 : if (DEBUGMEM>1) pari_warn(warnmem,"sum");
50 0 : p = gerepileupto(av, p);
51 : }
52 : }
53 53728 : return gerepileupto(av, p);
54 : }
55 :
56 : /*******************************************************************/
57 : /* */
58 : /* TRANSPOSE */
59 : /* */
60 : /*******************************************************************/
61 : /* A[x0,]~ */
62 : static GEN
63 16220532 : row_transpose(GEN A, long x0)
64 : {
65 16220532 : long i, lB = lg(A);
66 16220532 : GEN B = cgetg(lB, t_COL);
67 148100374 : for (i=1; i<lB; i++) gel(B, i) = gcoeff(A, x0, i);
68 16220584 : return B;
69 : }
70 : static GEN
71 18589 : row_transposecopy(GEN A, long x0)
72 : {
73 18589 : long i, lB = lg(A);
74 18589 : GEN B = cgetg(lB, t_COL);
75 161132 : for (i=1; i<lB; i++) gel(B, i) = gcopy(gcoeff(A, x0, i));
76 18589 : return B;
77 : }
78 :
79 : /* No copy*/
80 : GEN
81 4435030 : shallowtrans(GEN x)
82 : {
83 : long i, dx, lx;
84 : GEN y;
85 4435030 : switch(typ(x))
86 : {
87 1750 : case t_VEC: y = leafcopy(x); settyp(y,t_COL); break;
88 82216 : case t_COL: y = leafcopy(x); settyp(y,t_VEC); break;
89 4351125 : case t_MAT:
90 4351125 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
91 4350096 : dx = lgcols(x); y = cgetg(dx,t_MAT);
92 20571243 : for (i = 1; i < dx; i++) gel(y,i) = row_transpose(x,i);
93 4350755 : break;
94 0 : default: pari_err_TYPE("shallowtrans",x);
95 : return NULL;/*LCOV_EXCL_LINE*/
96 : }
97 4434721 : return y;
98 : }
99 :
100 : GEN
101 43520 : gtrans(GEN x)
102 : {
103 : long i, dx, lx;
104 : GEN y;
105 43520 : switch(typ(x))
106 : {
107 36694 : case t_VEC: y = gcopy(x); settyp(y,t_COL); break;
108 4872 : case t_COL: y = gcopy(x); settyp(y,t_VEC); break;
109 1947 : case t_MAT:
110 1947 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
111 1940 : dx = lgcols(x); y = cgetg(dx,t_MAT);
112 20529 : for (i = 1; i < dx; i++) gel(y,i) = row_transposecopy(x,i);
113 1940 : break;
114 7 : default: pari_err_TYPE("gtrans",x);
115 : return NULL;/*LCOV_EXCL_LINE*/
116 : }
117 43506 : return y;
118 : }
119 :
120 : /*******************************************************************/
121 : /* */
122 : /* EXTRACTION */
123 : /* */
124 : /*******************************************************************/
125 :
126 : static long
127 182 : str_to_long(char *s, char **pt)
128 : {
129 182 : long a = atol(s);
130 182 : while (isspace((int)*s)) s++;
131 182 : if (*s == '-' || *s == '+') s++;
132 385 : while (isdigit((int)*s) || isspace((int)*s)) s++;
133 182 : *pt = s; return a;
134 : }
135 :
136 : static int
137 112 : get_range(char *s, long *a, long *b, long *cmpl, long lx)
138 : {
139 112 : long max = lx - 1;
140 :
141 112 : *a = 1; *b = max;
142 112 : if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0;
143 112 : if (!*s) return 0;
144 112 : if (*s != '.')
145 : {
146 105 : *a = str_to_long(s, &s);
147 105 : if (*a < 0) *a += lx;
148 105 : if (*a<1 || *a>max) return 0;
149 : }
150 112 : if (*s == '.')
151 : {
152 105 : s++; if (*s != '.') return 0;
153 105 : do s++; while (isspace((int)*s));
154 105 : if (*s)
155 : {
156 77 : *b = str_to_long(s, &s);
157 77 : if (*b < 0) *b += lx;
158 77 : if (*b<1 || *b>max || *s) return 0;
159 : }
160 98 : return 1;
161 : }
162 7 : if (*s) return 0;
163 7 : *b = *a; return 1;
164 : }
165 :
166 : static int
167 35 : extract_selector_ok(long lx, GEN L)
168 : {
169 : long i, l;
170 35 : switch (typ(L))
171 : {
172 7 : case t_INT: {
173 : long maxj;
174 7 : if (!signe(L)) return 1;
175 7 : l = lgefint(L)-1;
176 7 : maxj = BITS_IN_LONG - bfffo(*int_MSW(L));
177 7 : return ((l-2) * BITS_IN_LONG + maxj < lx);
178 : }
179 7 : case t_STR: {
180 : long first, last, cmpl;
181 7 : return get_range(GSTR(L), &first, &last, &cmpl, lx);
182 : }
183 14 : case t_VEC: case t_COL:
184 14 : l = lg(L);
185 28 : for (i=1; i<l; i++)
186 : {
187 21 : long j = itos(gel(L,i));
188 21 : if (j>=lx || j<=0) return 0;
189 : }
190 7 : return 1;
191 7 : case t_VECSMALL:
192 7 : l = lg(L);
193 21 : for (i=1; i<l; i++)
194 : {
195 14 : long j = L[i];
196 14 : if (j>=lx || j<=0) return 0;
197 : }
198 7 : return 1;
199 : }
200 0 : return 0;
201 : }
202 :
203 : GEN
204 10234 : shallowmatextract(GEN x, GEN l1, GEN l2)
205 : {
206 10234 : long i, j, n1 = lg(l1), n2 = lg(l2);
207 10234 : GEN M = cgetg(n2, t_MAT);
208 64988 : for(i=1; i < n2; i++)
209 : {
210 54754 : long ii = l2[i];
211 54754 : GEN C = cgetg(n1, t_COL);
212 955696 : for (j=1; j < n1; j++)
213 : {
214 900942 : long jj = l1[j];
215 900942 : gel(C, j) = gmael(x, ii, jj);
216 : }
217 54754 : gel(M, i) = C;
218 : }
219 10234 : return M;
220 : }
221 :
222 : GEN
223 39470 : shallowextract(GEN x, GEN L)
224 : {
225 39470 : long i,j, tl = typ(L), tx = typ(x), lx = lg(x);
226 : GEN y;
227 :
228 39470 : switch(tx)
229 : {
230 39463 : case t_VEC:
231 : case t_COL:
232 : case t_MAT:
233 39463 : case t_VECSMALL: break;
234 7 : default: pari_err_TYPE("extract",x);
235 :
236 : }
237 39463 : if (tl==t_INT)
238 : { /* extract components of x as per the bits of mask L */
239 : long k, l, ix, iy, maxj;
240 : GEN Ld;
241 3311 : if (!signe(L)) return cgetg(1,tx);
242 3304 : y = new_chunk(lx);
243 3304 : l = lgefint(L)-1; ix = iy = 1;
244 3304 : maxj = BITS_IN_LONG - bfffo(*int_MSW(L));
245 3304 : if ((l-2) * BITS_IN_LONG + maxj >= lx)
246 7 : pari_err_TYPE("vecextract [mask too large]", L);
247 3672 : for (k = 2, Ld = int_LSW(L); k < l; k++, Ld = int_nextW(Ld))
248 : {
249 375 : ulong B = *Ld;
250 20439 : for (j = 0; j < BITS_IN_LONG; j++, B >>= 1, ix++)
251 20064 : if (B & 1) y[iy++] = x[ix];
252 : }
253 : { /* k = l */
254 3297 : ulong B = *Ld;
255 28551 : for (j = 0; j < maxj; j++, B >>= 1, ix++)
256 25254 : if (B & 1) y[iy++] = x[ix];
257 : }
258 3297 : y[0] = evaltyp(tx) | evallg(iy);
259 3297 : return y;
260 : }
261 36152 : if (tl==t_STR)
262 : {
263 105 : char *s = GSTR(L);
264 : long first, last, cmpl, d;
265 105 : if (! get_range(s, &first, &last, &cmpl, lx))
266 7 : pari_err_TYPE("vecextract [incorrect range]", L);
267 98 : if (lx == 1) return cgetg(1,tx);
268 98 : d = last - first;
269 98 : if (cmpl)
270 : {
271 21 : if (d >= 0)
272 : {
273 14 : y = cgetg(lx - (1+d),tx);
274 469 : for (j=1; j<first; j++) gel(y,j) = gel(x,j);
275 266 : for (i=last+1; i<lx; i++,j++) gel(y,j) = gel(x,i);
276 : }
277 : else
278 : {
279 7 : y = cgetg(lx - (1-d),tx);
280 14 : for (j=1,i=lx-1; i>first; i--,j++) gel(y,j) = gel(x,i);
281 14 : for (i=last-1; i>0; i--,j++) gel(y,j) = gel(x,i);
282 : }
283 : }
284 : else
285 : {
286 77 : if (d >= 0)
287 : {
288 35 : y = cgetg(d+2,tx);
289 112 : for (i=first,j=1; i<=last; i++,j++) gel(y,j) = gel(x,i);
290 : }
291 : else
292 : {
293 42 : y = cgetg(2-d,tx);
294 203 : for (i=first,j=1; i>=last; i--,j++) gel(y,j) = gel(x,i);
295 : }
296 : }
297 98 : return y;
298 : }
299 :
300 36047 : if (is_vec_t(tl))
301 : {
302 77 : long ll=lg(L); y=cgetg(ll,tx);
303 196 : for (i=1; i<ll; i++)
304 : {
305 133 : j = itos(gel(L,i));
306 133 : if (j<=0) pari_err_COMPONENT("vecextract","<=",gen_0,stoi(j));
307 126 : if (j>=lx) pari_err_COMPONENT("vecextract",">=",stoi(lx),stoi(j));
308 119 : gel(y,i) = gel(x,j);
309 : }
310 63 : return y;
311 : }
312 35970 : if (tl == t_VECSMALL)
313 : {
314 35963 : long ll=lg(L); y=cgetg(ll,tx);
315 161124 : for (i=1; i<ll; i++)
316 : {
317 125161 : j = L[i];
318 125161 : if (j<=0) pari_err_COMPONENT("vecextract","<=",gen_0,stoi(j));
319 125161 : if (j>=lx) pari_err_COMPONENT("vecextract",">=",stoi(lx),stoi(j));
320 125161 : gel(y,i) = gel(x,j);
321 : }
322 35963 : return y;
323 : }
324 7 : pari_err_TYPE("vecextract [mask]", L);
325 : return NULL; /* LCOV_EXCL_LINE */
326 : }
327 :
328 : /* does the component selector l select 0 component ? */
329 : static int
330 78 : select_0(GEN l)
331 : {
332 78 : switch(typ(l))
333 : {
334 14 : case t_INT:
335 14 : return (!signe(l));
336 43 : case t_VEC: case t_COL: case t_VECSMALL:
337 43 : return (lg(l) == 1);
338 : }
339 21 : return 0;
340 : }
341 :
342 : GEN
343 28674 : extract0(GEN x, GEN l1, GEN l2)
344 : {
345 28674 : pari_sp av = avma, av2;
346 : GEN y;
347 28674 : if (! l2)
348 : {
349 28596 : y = shallowextract(x, l1);
350 28554 : if (lg(y) == 1 || typ(y) == t_VECSMALL) return y;
351 28547 : av2 = avma;
352 28547 : y = gcopy(y);
353 : }
354 : else
355 : {
356 78 : if (typ(x) != t_MAT) pari_err_TYPE("extract",x);
357 78 : y = shallowextract(x,l2);
358 78 : if (select_0(l1)) { set_avma(av); return zeromat(0, lg(y)-1); }
359 64 : if (lg(y) == 1 && lg(x) > 1)
360 : {
361 35 : if (!extract_selector_ok(lgcols(x), l1))
362 7 : pari_err_TYPE("vecextract [incorrect mask]", l1);
363 28 : set_avma(av); return cgetg(1, t_MAT);
364 : }
365 29 : y = shallowextract(shallowtrans(y), l1);
366 29 : av2 = avma;
367 29 : y = gtrans(y);
368 : }
369 28576 : stackdummy(av, av2);
370 28576 : return y;
371 : }
372 :
373 : static long
374 5876 : vecslice_parse_arg(long lA, long *y1, long *y2, long *skip)
375 : {
376 5876 : *skip=0;
377 5876 : if (*y1==LONG_MAX)
378 : {
379 231 : if (*y2!=LONG_MAX)
380 : {
381 140 : if (*y2<0) *y2 += lA;
382 140 : if (*y2<0 || *y2==LONG_MAX || *y2>=lA)
383 0 : pari_err_DIM("_[..]");
384 140 : *skip=*y2;
385 : }
386 231 : *y1 = 1; *y2 = lA-1;
387 : }
388 5645 : else if (*y2==LONG_MAX) *y2 = *y1;
389 5876 : if (*y1<=0) *y1 += lA;
390 5876 : if (*y2<0) *y2 += lA;
391 5876 : if (*y1<=0 || *y1>*y2+1 || *y2>=lA) pari_err_DIM("_[..]");
392 5862 : return *y2 - *y1 + 2 - !!*skip;
393 : }
394 :
395 : static GEN
396 6359 : vecslice_i(GEN A, long t, long lB, long y1, long skip)
397 : {
398 6359 : GEN B = cgetg(lB, t);
399 : long i;
400 149884 : for (i=1; i<lB; i++, y1++)
401 : {
402 143525 : if (y1 == skip) { i--; continue; }
403 143385 : gel(B,i) = gcopy(gel(A,y1));
404 : }
405 6359 : return B;
406 : }
407 :
408 : static GEN
409 14 : rowslice_i(GEN A, long lB, long x1, long y1, long skip)
410 : {
411 14 : GEN B = cgetg(lB, t_VEC);
412 : long i;
413 77 : for (i=1; i<lB; i++, y1++)
414 : {
415 63 : if (y1 == skip) { i--; continue; }
416 56 : gel(B,i) = gcopy(gcoeff(A,x1,y1));
417 : }
418 14 : return B;
419 : }
420 :
421 : static GEN
422 0 : rowsmallslice_i(GEN A, long lB, long x1, long y1, long skip)
423 : {
424 0 : GEN B = cgetg(lB, t_VECSMALL);
425 : long i;
426 0 : for (i=1; i<lB; i++, y1++)
427 : {
428 0 : if (y1 == skip) { i--; continue; }
429 0 : B[i] = coeff(A,x1,y1);
430 : }
431 0 : return B;
432 : }
433 :
434 : static GEN
435 28 : vecsmallslice_i(GEN A, long t, long lB, long y1, long skip)
436 : {
437 28 : GEN B = cgetg(lB, t);
438 : long i;
439 126 : for (i=1; i<lB; i++, y1++)
440 : {
441 98 : if (y1 == skip) { i--; continue; }
442 91 : B[i] = A[y1];
443 : }
444 28 : return B;
445 : }
446 : GEN
447 5533 : vecslice0(GEN A, long y1, long y2)
448 : {
449 5533 : long skip, lB, t = typ(A);
450 5533 : switch(t)
451 : {
452 5435 : case t_VEC: case t_COL:
453 5435 : lB = vecslice_parse_arg(lg(A), &y1, &y2, &skip);
454 5421 : return vecslice_i(A, t,lB,y1,skip);
455 28 : case t_VECSMALL:
456 28 : lB = vecslice_parse_arg(lg(A), &y1, &y2, &skip);
457 28 : return vecsmallslice_i(A, t,lB,y1,skip);
458 63 : case t_LIST:
459 63 : if (list_typ(A) == t_LIST_RAW)
460 : {
461 63 : GEN y, z = list_data(A);
462 63 : long l = z? lg(z): 1;
463 63 : lB = vecslice_parse_arg(l, &y1, &y2, &skip);
464 63 : y = mklist(); if (!z) return y;
465 63 : list_data(y) = vecslice_i(z, t_VEC,lB,y1,skip);
466 63 : return y;
467 : }
468 : default:
469 7 : pari_err_TYPE("_[_.._]",A);
470 : return NULL;/*LCOV_EXCL_LINE*/
471 : }
472 : }
473 :
474 : GEN
475 182 : matslice0(GEN A, long x1, long x2, long y1, long y2)
476 : {
477 : GEN B;
478 182 : long i, lB, lA = lg(A), rA, t, skip, rskip, rlB;
479 182 : long is_col = y1!=LONG_MAX && y2==LONG_MAX;
480 182 : long is_row = x1!=LONG_MAX && x2==LONG_MAX;
481 : GEN (*slice)(GEN A, long t, long lB, long y1, long skip);
482 182 : if (typ(A)!=t_MAT) pari_err_TYPE("_[_.._,_.._]",A);
483 182 : lB = vecslice_parse_arg(lA, &y1, &y2, &skip);
484 182 : if (is_col) return vecslice0(gel(A, y1), x1, x2);
485 168 : rA = lg(A)==1 ? 1: lgcols(A);
486 168 : rlB = vecslice_parse_arg(rA, &x1, &x2, &rskip);
487 168 : t = lg(A)==1 ? t_COL: typ(gel(A,1));
488 168 : if (is_row) return t == t_COL ? rowslice_i(A, lB, x1, y1, skip):
489 0 : rowsmallslice_i(A, lB, x1, y1, skip);
490 154 : slice = t == t_COL? &vecslice_i: &vecsmallslice_i;
491 :
492 154 : B = cgetg(lB, t_MAT);
493 1043 : for (i=1; i<lB; i++, y1++)
494 : {
495 889 : if (y1 == skip) { i--; continue; }
496 875 : gel(B,i) = slice(gel(A,y1),t,rlB, x1, rskip);
497 : }
498 154 : return B;
499 : }
500 :
501 : GEN
502 10472 : vecrange(GEN a, GEN b)
503 : {
504 : GEN y;
505 : long i, l;
506 10472 : if (typ(a)!=t_INT) pari_err_TYPE("[_.._]",a);
507 10465 : if (typ(b)!=t_INT) pari_err_TYPE("[_.._]",b);
508 10458 : if (cmpii(a,b)>0) return cgetg(1,t_VEC);
509 10451 : l = itos(subii(b,a))+1;
510 10451 : a = setloop(a);
511 10451 : y = cgetg(l+1, t_VEC);
512 26377185 : for (i=1; i<=l; a = incloop(a), i++) gel(y,i) = icopy(a);
513 10451 : return y;
514 : }
515 :
516 : GEN
517 0 : vecrangess(long a, long b)
518 : {
519 : GEN y;
520 : long i, l;
521 0 : if (a>b) return cgetg(1,t_VEC);
522 0 : l = b-a+1;
523 0 : y = cgetg(l+1, t_VEC);
524 0 : for (i=1; i<=l; a++, i++) gel(y,i) = stoi(a);
525 0 : return y;
526 : }
527 :
528 : GEN
529 96 : genindexselect(void *E, long (*f)(void* E, GEN x), GEN A)
530 : {
531 : long l, i, lv;
532 : GEN v, z;
533 : pari_sp av;
534 96 : clone_lock(A);
535 96 : switch(typ(A))
536 : {
537 7 : case t_LIST:
538 7 : z = list_data(A);
539 7 : l = z? lg(z): 1;
540 7 : break;
541 82 : case t_VEC: case t_COL: case t_MAT:
542 82 : l = lg(A);
543 82 : z = A;
544 82 : break;
545 7 : default:
546 7 : pari_err_TYPE("select",A);
547 : return NULL;/*LCOV_EXCL_LINE*/
548 : }
549 89 : v = cgetg(l, t_VECSMALL);
550 89 : av = avma;
551 12801 : for (i = lv = 1; i < l; i++) {
552 12712 : if (f(E, gel(z,i))) v[lv++] = i;
553 12712 : set_avma(av);
554 : }
555 89 : clone_unlock_deep(A); fixlg(v, lv); return v;
556 : }
557 : static GEN
558 88 : extract_copy(GEN A, GEN v)
559 : {
560 88 : long i, l = lg(v);
561 88 : GEN B = cgetg(l, typ(A));
562 387 : for (i = 1; i < l; i++) gel(B,i) = gcopy(gel(A,v[i]));
563 88 : return B;
564 : }
565 : /* as genselect, but treat A [ t_VEC,t_COL, or t_MAT] as a t_VEC */
566 : GEN
567 0 : vecselect(void *E, long (*f)(void* E, GEN x), GEN A)
568 : {
569 : GEN v;
570 0 : clone_lock(A);
571 0 : v = genindexselect(E, f, A);
572 0 : A = extract_copy(A, v); settyp(A, t_VEC);
573 0 : clone_unlock_deep(A); return A;
574 : }
575 : GEN
576 83 : genselect(void *E, long (*f)(void* E, GEN x), GEN A)
577 : {
578 : GEN y, z, v;/* v left on stack for efficiency */
579 83 : clone_lock(A);
580 83 : switch(typ(A))
581 : {
582 14 : case t_LIST:
583 14 : z = list_data(A);
584 14 : if (!z) y = mklist();
585 : else
586 : {
587 : GEN B;
588 14 : y = cgetg(3, t_LIST);
589 14 : v = genindexselect(E, f, z);
590 14 : B = extract_copy(z, v);
591 14 : y[1] = lg(B)-1;
592 14 : list_data(y) = B;
593 : }
594 14 : break;
595 62 : case t_VEC: case t_COL: case t_MAT:
596 62 : v = genindexselect(E, f, A);
597 62 : y = extract_copy(A, v);
598 62 : break;
599 7 : default:
600 7 : pari_err_TYPE("select",A);
601 : return NULL;/*LCOV_EXCL_LINE*/
602 : }
603 76 : clone_unlock_deep(A); return y;
604 : }
605 :
606 : static void
607 54683 : check_callgen1(GEN f, const char *s)
608 : {
609 54683 : if (typ(f) != t_CLOSURE || closure_is_variadic(f) || closure_arity(f) < 1)
610 0 : pari_err_TYPE(s, f);
611 54683 : }
612 :
613 : GEN
614 103 : select0(GEN f, GEN x, long flag)
615 : {
616 103 : check_callgen1(f, "select");
617 103 : switch(flag)
618 : {
619 83 : case 0: return genselect((void *) f, gp_callbool, x);
620 20 : case 1: return genindexselect((void *) f, gp_callbool, x);
621 0 : default: pari_err_FLAG("select");
622 : return NULL;/*LCOV_EXCL_LINE*/
623 : }
624 : }
625 :
626 : GEN
627 23920 : parselect_worker(GEN d, GEN C)
628 : {
629 23920 : return gequal0(closure_callgen1(C, d))? gen_0: gen_1;
630 : }
631 :
632 : GEN
633 24 : parselect(GEN C, GEN D, long flag)
634 : {
635 : pari_sp av;
636 24 : long lv, l = lg(D), i;
637 : GEN V, W, worker;
638 24 : check_callgen1(C, "parselect");
639 24 : if (!is_vec_t(typ(D))) pari_err_TYPE("parselect",D);
640 24 : W = cgetg(l, t_VECSMALL); av = avma;
641 24 : worker = snm_closure(is_entry("_parselect_worker"), mkvec(C));
642 24 : V = gen_parapply(worker, D);
643 24048 : for (lv=1, i=1; i<l; i++)
644 24024 : if (signe(gel(V,i))) W[lv++] = i;
645 24 : fixlg(W, lv);
646 24 : set_avma(av);
647 24 : return flag? W: extract_copy(D, W);
648 : }
649 :
650 : GEN
651 0 : veccatapply(void *E, GEN (*f)(void* E, GEN x), GEN x)
652 : {
653 0 : pari_sp av = avma;
654 0 : GEN v = vecapply(E, f, x);
655 0 : return lg(v) == 1? v: gerepilecopy(av, shallowconcat1(v));
656 : }
657 :
658 : static GEN
659 14 : vecapply2(void *E, GEN (*f)(void* E, GEN x), GEN x)
660 : {
661 : long i, lx;
662 14 : GEN y = cgetg_copy(x, &lx); y[1] = x[1];
663 56 : for (i=2; i<lx; i++) gel(y,i) = f(E, gel(x,i));
664 14 : return y;
665 : }
666 : static GEN
667 167777 : vecapply1(void *E, GEN (*f)(void* E, GEN x), GEN x)
668 : {
669 : long i, lx;
670 167777 : GEN y = cgetg_copy(x, &lx);
671 861259 : for (i=1; i<lx; i++) gel(y,i) = f(E, gel(x,i));
672 167777 : return y;
673 : }
674 : static GEN
675 7 : mapapply1(void *E, GEN (*f)(void* E, GEN x), GEN x)
676 : {
677 : long i, lx;
678 7 : GEN y = cgetg_copy(x, &lx);
679 28 : for (i=1; i<lx; i++)
680 : {
681 21 : GEN xi = gel(x, i);
682 21 : gel(y,i) = mkvec2(mkvec2(gcopy(gmael(xi, 1, 1)), f(E,gmael(xi, 1, 2))),
683 21 : gcopy(gel(xi, 2)));
684 : }
685 7 : return y;
686 : }
687 : /* as genapply, but treat A [ t_VEC,t_COL, or t_MAT] as a t_VEC */
688 : GEN
689 115070 : vecapply(void *E, GEN (*f)(void* E, GEN x), GEN x)
690 : {
691 : GEN y;
692 115070 : clone_lock(x); y = vecapply1(E,f,x);
693 115070 : clone_unlock_deep(x); settyp(y, t_VEC); return y;
694 : }
695 : GEN
696 52693 : genapply(void *E, GEN (*f)(void* E, GEN x), GEN x)
697 : {
698 52693 : long i, lx, tx = typ(x);
699 : GEN y, z;
700 52693 : if (is_scalar_t(tx)) return f(E, x);
701 52693 : clone_lock(x);
702 52693 : switch(tx) {
703 7 : case t_POL: y = normalizepol(vecapply2(E,f,x)); break;
704 7 : case t_SER:
705 7 : y = ser_isexactzero(x)? gcopy(x): normalizeser(vecapply2(E,f,x));
706 7 : break;
707 28 : case t_LIST:
708 : {
709 28 : long t = list_typ(x);
710 28 : z = list_data(x);
711 28 : if (!z)
712 7 : y = mklist_typ(t);
713 : else
714 : {
715 21 : y = cgetg(3, t_LIST);
716 21 : y[1] = evaltyp(t)|evallg(lg(z)-1);
717 : switch(t)
718 : {
719 14 : case t_LIST_RAW:
720 14 : list_data(y) = vecapply1(E,f,z);
721 14 : break;
722 7 : case t_LIST_MAP:
723 7 : list_data(y) = mapapply1(E,f,z);
724 7 : break;
725 : }
726 28 : }
727 : }
728 28 : break;
729 42 : case t_MAT:
730 42 : y = cgetg_copy(x, &lx);
731 126 : for (i = 1; i < lx; i++) gel(y,i) = vecapply1(E,f,gel(x,i));
732 42 : break;
733 :
734 52609 : case t_VEC: case t_COL: y = vecapply1(E,f,x); break;
735 0 : default:
736 0 : pari_err_TYPE("apply",x);
737 : return NULL;/*LCOV_EXCL_LINE*/
738 : }
739 52693 : clone_unlock_deep(x); return y;
740 : }
741 :
742 : GEN
743 52693 : apply0(GEN f, GEN x)
744 : {
745 52693 : check_callgen1(f, "apply");
746 52693 : return genapply((void *) f, gp_call, x);
747 : }
748 :
749 : GEN
750 448 : vecselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
751 : GEN (*fun)(void* E, GEN x), GEN A)
752 : {
753 : GEN y;
754 448 : long i, l = lg(A), nb=1;
755 448 : clone_lock(A); y = cgetg(l, t_VEC);
756 26165475 : for (i=1; i<l; i++)
757 26165027 : if (pred(Epred, gel(A,i))) gel(y,nb++) = fun(Efun, gel(A,i));
758 448 : fixlg(y,nb); clone_unlock_deep(A); return y;
759 : }
760 :
761 : GEN
762 0 : veccatselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
763 : GEN (*fun)(void* E, GEN x), GEN A)
764 : {
765 0 : pari_sp av = avma;
766 0 : GEN v = vecselapply(Epred, pred, Efun, fun, A);
767 0 : return lg(v) == 1? v: gerepilecopy(av, shallowconcat1(v));
768 : }
769 :
770 : /* suitable for gerepileupto */
771 : GEN
772 857 : parapply_slice_worker(GEN D, GEN worker)
773 : {
774 : long l, i;
775 857 : GEN w = cgetg_copy(D, &l);
776 73374 : for (i = 1; i < l; i++) gel(w,i) = closure_callgen1(worker, gel(D,i));
777 775 : return w;
778 : }
779 :
780 : /* B <- {A[i] : i = r (mod m)}, 1 <= r <= m */
781 : static void
782 72575 : arithprogset(GEN B, GEN A, long r, long m)
783 : {
784 72575 : long i, k, l = lg(A);
785 165178 : for (k = 1, i = r; i < l; i += m, k++) gel(B, k) = gel(A, i);
786 72575 : setlg(B, k);
787 72575 : }
788 : GEN
789 6736 : gen_parapply_slice(GEN worker, GEN D, long mmin)
790 : {
791 6736 : long l, r, n = lg(D)-1, m = minss(mmin, n), pending = 0;
792 6736 : GEN L = cgetg(n / m + 2, t_VEC), va = mkvec(L), V = cgetg_copy(D, &l);
793 : struct pari_mt pt;
794 6736 : mt_queue_start_lim(&pt, worker, m);
795 145150 : for (r = 1; r <= m || pending; r++)
796 : {
797 : long workid;
798 : GEN done;
799 138414 : if (r <= m) arithprogset(L, D, r, m);
800 138414 : mt_queue_submit(&pt, r, r <= m? va: NULL);
801 138414 : done = mt_queue_get(&pt, &workid, &pending);
802 138414 : if (done)
803 : {
804 72575 : long j, k, J = lg(done)-1;
805 165178 : for (j = 1, k = workid; j <= J; j++, k +=m) gel(V, k) = gel(done, j);
806 : }
807 : }
808 6736 : mt_queue_end(&pt); return V;
809 : }
810 :
811 : GEN
812 20231 : gen_parapply(GEN worker, GEN D)
813 : {
814 20231 : long l = lg(D), i, pending = 0;
815 : GEN W, V;
816 : struct pari_mt pt;
817 :
818 20231 : if (l == 1) return cgetg(1, typ(D));
819 20203 : W = cgetg(2, t_VEC);
820 20203 : V = cgetg(l, typ(D));
821 20203 : mt_queue_start_lim(&pt, worker, l-1);
822 374851 : for (i = 1; i < l || pending; i++)
823 : {
824 : long workid;
825 : GEN done;
826 354649 : if (i < l) gel(W,1) = gel(D,i);
827 354649 : mt_queue_submit(&pt, i, i<l? W: NULL);
828 354648 : done = mt_queue_get(&pt, &workid, &pending);
829 354648 : if (done) gel(V,workid) = done;
830 : }
831 20202 : mt_queue_end(&pt); return V;
832 : }
833 :
834 : GEN
835 1863 : parapply(GEN C, GEN D)
836 : {
837 1863 : pari_sp av = avma;
838 1863 : check_callgen1(C, "parapply");
839 1863 : if (!is_vec_t(typ(D))) pari_err_TYPE("parapply",D);
840 1863 : return gerepileupto(av, gen_parapply(C, D));
841 : }
842 :
843 : GEN
844 28 : genfold(void *E, GEN (*f)(void* E, GEN x, GEN y), GEN x)
845 : {
846 28 : pari_sp av = avma;
847 : GEN z;
848 28 : long i, l = lg(x);
849 28 : if (!is_vec_t(typ(x))|| l==1 ) pari_err_TYPE("fold",x);
850 28 : clone_lock(x);
851 28 : z = gel(x,1);
852 119 : for (i=2; i<l; i++)
853 : {
854 91 : z = f(E,z,gel(x,i));
855 91 : if (gc_needed(av, 2))
856 : {
857 0 : if (DEBUGMEM>1) pari_warn(warnmem,"fold");
858 0 : z = gerepilecopy(av, z);
859 : }
860 : }
861 28 : clone_unlock_deep(x);
862 28 : return gerepilecopy(av, z);
863 : }
864 :
865 : GEN
866 28 : fold0(GEN f, GEN x)
867 : {
868 28 : if (typ(f) != t_CLOSURE || closure_arity(f) < 2) pari_err_TYPE("fold",f);
869 28 : return genfold((void *) f, gp_call2, x);
870 : }
871 : /*******************************************************************/
872 : /* */
873 : /* SCALAR-MATRIX OPERATIONS */
874 : /* */
875 : /*******************************************************************/
876 : GEN
877 374926 : gtomat(GEN x)
878 : {
879 : long lx, i;
880 : GEN y;
881 :
882 374926 : if (!x) return cgetg(1, t_MAT);
883 374905 : switch(typ(x))
884 : {
885 28 : case t_LIST:
886 28 : if (list_typ(x)==t_LIST_MAP)
887 14 : return maptomat(x);
888 14 : x = list_data(x);
889 14 : if (!x) return cgetg(1, t_MAT);
890 : /* fall through */
891 : case t_VEC: {
892 2734 : lx=lg(x); y=cgetg(lx,t_MAT);
893 2737 : if (lx == 1) break;
894 2737 : if (typ(gel(x,1)) == t_COL) {
895 2394 : long h = lgcols(x);
896 5649 : for (i=2; i<lx; i++) {
897 3255 : if (typ(gel(x,i)) != t_COL || lg(gel(x,i)) != h) break;
898 : }
899 2394 : if (i == lx) { /* matrix with h-1 rows */
900 2394 : y = cgetg(lx, t_MAT);
901 8043 : for (i=1 ; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
902 2394 : return y;
903 : }
904 : }
905 1029 : for (i=1; i<lx; i++) gel(y,i) = mkcolcopy(gel(x,i));
906 343 : break;
907 : }
908 117969 : case t_COL:
909 117969 : lx = lg(x);
910 117969 : if (lx == 1) return cgetg(1, t_MAT);
911 117955 : if (typ(gel(x,1)) == t_VEC) {
912 7 : long j, h = lg(gel(x,1));
913 14 : for (i=2; i<lx; i++) {
914 7 : if (typ(gel(x,i)) != t_VEC || lg(gel(x,i)) != h) break;
915 : }
916 7 : if (i == lx) { /* matrix with h cols */
917 7 : y = cgetg(h, t_MAT);
918 28 : for (j=1 ; j<h; j++) {
919 21 : gel(y,j) = cgetg(lx, t_COL);
920 63 : for (i=1; i<lx; i++) gcoeff(y,i,j) = gcopy(gmael(x,i,j));
921 : }
922 7 : return y;
923 : }
924 : }
925 117948 : y = mkmatcopy(x); break;
926 253355 : case t_MAT:
927 253355 : y = gcopy(x); break;
928 7 : case t_QFB: {
929 : GEN b;
930 7 : y = cgetg(3,t_MAT); b = gmul2n(gel(x,2),-1);
931 7 : gel(y,1) = mkcol2(icopy(gel(x,1)), b);
932 7 : gel(y,2) = mkcol2(b, icopy(gel(x,3)));
933 7 : break;
934 : }
935 819 : default:
936 819 : y = cgetg(2,t_MAT); gel(y,1) = mkcolcopy(x);
937 819 : break;
938 : }
939 372478 : return y;
940 : }
941 :
942 : /* create the diagonal matrix, whose diagonal is given by x */
943 : GEN
944 798 : diagonal(GEN x)
945 : {
946 798 : long j, lx, tx = typ(x);
947 : GEN y;
948 :
949 798 : if (! is_matvec_t(tx)) return scalarmat(x,1);
950 791 : if (tx==t_MAT)
951 : {
952 14 : if (RgM_isdiagonal(x)) return gcopy(x);
953 7 : pari_err_TYPE("diagonal",x);
954 : }
955 777 : lx=lg(x); y=cgetg(lx,t_MAT);
956 2296 : for (j=1; j<lx; j++)
957 : {
958 1519 : gel(y,j) = zerocol(lx-1);
959 1519 : gcoeff(y,j,j) = gcopy(gel(x,j));
960 : }
961 777 : return y;
962 : }
963 : /* same, assuming x is a t_VEC/t_COL. Not memory clean. */
964 : GEN
965 348937 : diagonal_shallow(GEN x)
966 : {
967 348937 : long j, lx = lg(x);
968 348937 : GEN y = cgetg(lx,t_MAT);
969 :
970 959708 : for (j=1; j<lx; j++)
971 : {
972 610746 : gel(y,j) = zerocol(lx-1);
973 610751 : gcoeff(y,j,j) = gel(x,j);
974 : }
975 348962 : return y;
976 : }
977 :
978 : GEN
979 448 : zv_diagonal(GEN x)
980 : {
981 448 : long j, l = lg(x), n = l-1;
982 448 : GEN y = cgetg(l,t_MAT);
983 :
984 1428 : for (j = 1; j < l; j++)
985 : {
986 980 : gel(y,j) = zero_Flv(n);
987 980 : ucoeff(y,j,j) = uel(x,j);
988 : }
989 448 : return y;
990 : }
991 :
992 : /* compute m*diagonal(d) */
993 : GEN
994 70 : matmuldiagonal(GEN m, GEN d)
995 : {
996 : long j, lx;
997 70 : GEN y = cgetg_copy(m, &lx);
998 :
999 70 : if (typ(m)!=t_MAT) pari_err_TYPE("matmuldiagonal",m);
1000 70 : if (! is_vec_t(typ(d))) pari_err_TYPE("matmuldiagonal",d);
1001 70 : if (lg(d) != lx) pari_err_OP("operation 'matmuldiagonal'", m,d);
1002 350 : for (j=1; j<lx; j++) gel(y,j) = RgC_Rg_mul(gel(m,j), gel(d,j));
1003 70 : return y;
1004 : }
1005 :
1006 : /* compute A*B assuming the result is a diagonal matrix */
1007 : GEN
1008 7 : matmultodiagonal(GEN A, GEN B)
1009 : {
1010 7 : long i, j, hA, hB, lA = lg(A), lB = lg(B);
1011 7 : GEN y = matid(lB-1);
1012 :
1013 7 : if (typ(A) != t_MAT) pari_err_TYPE("matmultodiagonal",A);
1014 7 : if (typ(B) != t_MAT) pari_err_TYPE("matmultodiagonal",B);
1015 7 : hA = (lA == 1)? lB: lgcols(A);
1016 7 : hB = (lB == 1)? lA: lgcols(B);
1017 7 : if (lA != hB || lB != hA) pari_err_OP("operation 'matmultodiagonal'", A,B);
1018 56 : for (i=1; i<lB; i++)
1019 : {
1020 49 : GEN z = gen_0;
1021 392 : for (j=1; j<lA; j++) z = gadd(z, gmul(gcoeff(A,i,j),gcoeff(B,j,i)));
1022 49 : gcoeff(y,i,i) = z;
1023 : }
1024 7 : return y;
1025 : }
1026 :
1027 : /* [m[1,1], ..., m[l,l]], internal */
1028 : GEN
1029 854066 : RgM_diagonal_shallow(GEN m)
1030 : {
1031 854066 : long i, lx = lg(m);
1032 854066 : GEN y = cgetg(lx,t_VEC);
1033 2814993 : for (i=1; i<lx; i++) gel(y, i) = gcoeff(m,i,i);
1034 854080 : return y;
1035 : }
1036 :
1037 : /* same, public function */
1038 : GEN
1039 0 : RgM_diagonal(GEN m)
1040 : {
1041 0 : long i, lx = lg(m);
1042 0 : GEN y = cgetg(lx,t_VEC);
1043 0 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gcoeff(m,i,i));
1044 0 : return y;
1045 : }
|
