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 127698 : vecsum(GEN v)
34 : {
35 127698 : pari_sp av = avma;
36 : long i, l;
37 : GEN p;
38 127698 : if (!is_vec_t(typ(v)))
39 6 : pari_err_TYPE("vecsum", v);
40 127692 : l = lg(v);
41 127692 : if (l == 1) return gen_0;
42 127686 : p = gel(v,1);
43 127686 : if (l == 2) return gcopy(p);
44 212438 : for (i=2; i<l; i++)
45 : {
46 143591 : p = gadd(p, gel(v,i));
47 143591 : if (gc_needed(av, 2))
48 : {
49 0 : if (DEBUGMEM>1) pari_warn(warnmem,"sum");
50 0 : p = gc_upto(av, p);
51 : }
52 : }
53 68847 : return gc_upto(av, p);
54 : }
55 :
56 : /*******************************************************************/
57 : /* */
58 : /* TRANSPOSE */
59 : /* */
60 : /*******************************************************************/
61 : /* A[x0,]~ */
62 : static GEN
63 22004840 : row_transpose(GEN A, long x0)
64 : {
65 22004840 : long i, lB = lg(A);
66 22004840 : GEN B = cgetg(lB, t_COL);
67 169633187 : for (i=1; i<lB; i++) gel(B, i) = gcoeff(A, x0, i);
68 22004840 : return B;
69 : }
70 : static GEN
71 13701 : row_transposecopy(GEN A, long x0)
72 : {
73 13701 : long i, lB = lg(A);
74 13701 : GEN B = cgetg(lB, t_COL);
75 120306 : for (i=1; i<lB; i++) gel(B, i) = gcopy(gcoeff(A, x0, i));
76 13701 : return B;
77 : }
78 :
79 : /* No copy*/
80 : GEN
81 5999275 : shallowtrans(GEN x)
82 : {
83 : long i, dx, lx;
84 : GEN y;
85 5999275 : switch(typ(x))
86 : {
87 1079 : case t_VEC: y = leafcopy(x); settyp(y,t_COL); break;
88 135 : case t_COL: y = leafcopy(x); settyp(y,t_VEC); break;
89 5998061 : case t_MAT:
90 5998061 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
91 5996939 : dx = lgcols(x); y = cgetg(dx,t_MAT);
92 28001779 : for (i = 1; i < dx; i++) gel(y,i) = row_transpose(x,i);
93 5996939 : break;
94 0 : default: pari_err_TYPE("shallowtrans",x);
95 : return NULL;/*LCOV_EXCL_LINE*/
96 : }
97 5998153 : return y;
98 : }
99 :
100 : GEN
101 36625 : gtrans(GEN x)
102 : {
103 : long i, dx, lx;
104 : GEN y;
105 36625 : switch(typ(x))
106 : {
107 31614 : case t_VEC: y = gcopy(x); settyp(y,t_COL); break;
108 3573 : case t_COL: y = gcopy(x); settyp(y,t_VEC); break;
109 1432 : case t_MAT:
110 1432 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
111 1426 : dx = lgcols(x); y = cgetg(dx,t_MAT);
112 15127 : for (i = 1; i < dx; i++) gel(y,i) = row_transposecopy(x,i);
113 1426 : break;
114 6 : default: pari_err_TYPE("gtrans",x);
115 : return NULL;/*LCOV_EXCL_LINE*/
116 : }
117 36613 : return y;
118 : }
119 :
120 : /*******************************************************************/
121 : /* */
122 : /* EXTRACTION */
123 : /* */
124 : /*******************************************************************/
125 :
126 : static long
127 156 : str_to_long(char *s, char **pt)
128 : {
129 156 : long a = atol(s);
130 156 : while (isspace((unsigned char)*s)) s++;
131 156 : if (*s == '-' || *s == '+') s++;
132 330 : while (isdigit((unsigned char)*s) || isspace((unsigned char)*s)) s++;
133 156 : *pt = s; return a;
134 : }
135 :
136 : static int
137 96 : get_range(char *s, long *a, long *b, long *cmpl, long lx)
138 : {
139 96 : long max = lx - 1;
140 :
141 96 : *a = 1; *b = max;
142 96 : if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0;
143 96 : if (!*s) return 0;
144 96 : if (*s != '.')
145 : {
146 90 : *a = str_to_long(s, &s);
147 90 : if (*a < 0) *a += lx;
148 90 : if (*a<1 || *a>max) return 0;
149 : }
150 96 : if (*s == '.')
151 : {
152 90 : s++; if (*s != '.') return 0;
153 90 : do s++; while (isspace((unsigned char)*s));
154 90 : if (*s)
155 : {
156 66 : *b = str_to_long(s, &s);
157 66 : if (*b < 0) *b += lx;
158 66 : if (*b<1 || *b>max || *s) return 0;
159 : }
160 84 : return 1;
161 : }
162 6 : if (*s) return 0;
163 6 : *b = *a; return 1;
164 : }
165 :
166 : static int
167 30 : extract_selector_ok(long lx, GEN L)
168 : {
169 : long i, l;
170 30 : switch (typ(L))
171 : {
172 6 : case t_INT: {
173 : long maxj;
174 6 : if (!signe(L)) return 1;
175 6 : l = lgefint(L)-1;
176 6 : maxj = BITS_IN_LONG - bfffo(*int_MSW(L));
177 6 : return ((l-2) * BITS_IN_LONG + maxj < lx);
178 : }
179 6 : case t_STR: {
180 : long first, last, cmpl;
181 6 : return get_range(GSTR(L), &first, &last, &cmpl, lx);
182 : }
183 12 : case t_VEC: case t_COL:
184 12 : l = lg(L);
185 24 : for (i=1; i<l; i++)
186 : {
187 18 : long j = itos(gel(L,i));
188 18 : if (j>=lx || j<=0) return 0;
189 : }
190 6 : return 1;
191 6 : case t_VECSMALL:
192 6 : l = lg(L);
193 18 : for (i=1; i<l; i++)
194 : {
195 12 : long j = L[i];
196 12 : if (j>=lx || j<=0) return 0;
197 : }
198 6 : return 1;
199 : }
200 0 : return 0;
201 : }
202 :
203 : GEN
204 11327 : shallowmatextract(GEN x, GEN l1, GEN l2)
205 : {
206 11327 : long i, j, n1 = lg(l1), n2 = lg(l2);
207 11327 : GEN M = cgetg(n2, t_MAT);
208 65860 : for(i=1; i < n2; i++)
209 : {
210 54533 : long ii = l2[i];
211 54533 : GEN C = cgetg(n1, t_COL);
212 880150 : for (j=1; j < n1; j++)
213 : {
214 825617 : long jj = l1[j];
215 825617 : gel(C, j) = gmael(x, ii, jj);
216 : }
217 54533 : gel(M, i) = C;
218 : }
219 11327 : return M;
220 : }
221 :
222 : GEN
223 45057 : shallowextract(GEN x, GEN L)
224 : {
225 45057 : long i,j, tl = typ(L), tx = typ(x), lx = lg(x);
226 : GEN y;
227 :
228 45057 : switch(tx)
229 : {
230 45051 : case t_VEC:
231 : case t_COL:
232 : case t_MAT:
233 45051 : case t_VECSMALL: break;
234 6 : default: pari_err_TYPE("extract",x);
235 :
236 : }
237 45051 : 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 2862 : if (!signe(L)) return cgetg(1,tx);
242 2856 : y = new_chunk(lx);
243 2856 : l = lgefint(L)-1; ix = iy = 1;
244 2856 : maxj = BITS_IN_LONG - bfffo(*int_MSW(L));
245 2856 : if ((l-2) * BITS_IN_LONG + maxj >= lx)
246 6 : pari_err_TYPE("vecextract [mask too large]", L);
247 3207 : for (k = 2, Ld = int_LSW(L); k < l; k++, Ld = int_nextW(Ld))
248 : {
249 357 : ulong B = *Ld;
250 18981 : for (j = 0; j < BITS_IN_LONG; j++, B >>= 1, ix++)
251 18624 : if (B & 1) y[iy++] = x[ix];
252 : }
253 : { /* k = l */
254 2850 : ulong B = *Ld;
255 25008 : for (j = 0; j < maxj; j++, B >>= 1, ix++)
256 22158 : if (B & 1) y[iy++] = x[ix];
257 : }
258 2850 : y[0] = evaltyp(tx) | evallg(iy);
259 2850 : return y;
260 : }
261 42189 : if (tl==t_STR)
262 : {
263 90 : char *s = GSTR(L);
264 : long first, last, cmpl, d;
265 90 : if (! get_range(s, &first, &last, &cmpl, lx))
266 6 : pari_err_TYPE("vecextract [incorrect range]", L);
267 84 : if (lx == 1) return cgetg(1,tx);
268 84 : d = last - first;
269 84 : if (cmpl)
270 : {
271 18 : if (d >= 0)
272 : {
273 12 : y = cgetg(lx - (1+d),tx);
274 402 : for (j=1; j<first; j++) gel(y,j) = gel(x,j);
275 228 : for (i=last+1; i<lx; i++,j++) gel(y,j) = gel(x,i);
276 : }
277 : else
278 : {
279 6 : y = cgetg(lx - (1-d),tx);
280 12 : for (j=1,i=lx-1; i>first; i--,j++) gel(y,j) = gel(x,i);
281 12 : for (i=last-1; i>0; i--,j++) gel(y,j) = gel(x,i);
282 : }
283 : }
284 : else
285 : {
286 66 : if (d >= 0)
287 : {
288 30 : y = cgetg(d+2,tx);
289 96 : for (i=first,j=1; i<=last; i++,j++) gel(y,j) = gel(x,i);
290 : }
291 : else
292 : {
293 36 : y = cgetg(2-d,tx);
294 174 : for (i=first,j=1; i>=last; i--,j++) gel(y,j) = gel(x,i);
295 : }
296 : }
297 84 : return y;
298 : }
299 :
300 42099 : if (is_vec_t(tl))
301 : {
302 66 : long ll=lg(L); y=cgetg(ll,tx);
303 168 : for (i=1; i<ll; i++)
304 : {
305 114 : j = itos(gel(L,i));
306 114 : if (j<=0) pari_err_COMPONENT("vecextract","<=",gen_0,stoi(j));
307 108 : if (j>=lx) pari_err_COMPONENT("vecextract",">=",stoi(lx),stoi(j));
308 102 : gel(y,i) = gel(x,j);
309 : }
310 54 : return y;
311 : }
312 42033 : if (tl == t_VECSMALL)
313 : {
314 42027 : long ll=lg(L); y=cgetg(ll,tx);
315 178595 : for (i=1; i<ll; i++)
316 : {
317 136568 : j = L[i];
318 136568 : if (j<=0) pari_err_COMPONENT("vecextract","<=",gen_0,stoi(j));
319 136568 : if (j>=lx) pari_err_COMPONENT("vecextract",">=",stoi(lx),stoi(j));
320 136568 : gel(y,i) = gel(x,j);
321 : }
322 42027 : return y;
323 : }
324 6 : 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 73 : select_0(GEN l)
331 : {
332 73 : switch(typ(l))
333 : {
334 12 : case t_INT:
335 12 : return (!signe(l));
336 43 : case t_VEC: case t_COL: case t_VECSMALL:
337 43 : return (lg(l) == 1);
338 : }
339 18 : return 0;
340 : }
341 :
342 : GEN
343 34146 : extract0(GEN x, GEN l1, GEN l2)
344 : {
345 34146 : pari_sp av = avma, av2;
346 : GEN y;
347 34146 : if (! l2)
348 : {
349 34073 : y = shallowextract(x, l1);
350 34037 : if (lg(y) == 1 || typ(y) == t_VECSMALL) return y;
351 34031 : av2 = avma;
352 34031 : y = gcopy(y);
353 : }
354 : else
355 : {
356 73 : if (typ(x) != t_MAT) pari_err_TYPE("extract",x);
357 73 : y = shallowextract(x,l2);
358 73 : if (select_0(l1)) { set_avma(av); return zeromat(0, lg(y)-1); }
359 61 : if (lg(y) == 1 && lg(x) > 1)
360 : {
361 30 : if (!extract_selector_ok(lgcols(x), l1))
362 6 : pari_err_TYPE("vecextract [incorrect mask]", l1);
363 24 : retgc_const(av, cgetg(1, t_MAT));
364 : }
365 31 : y = shallowextract(shallowtrans(y), l1);
366 31 : av2 = avma;
367 31 : y = gtrans(y);
368 : }
369 34062 : stackdummy(av, av2);
370 34062 : return y;
371 : }
372 :
373 : static long
374 1691 : vecslice_parse_arg(long lA, long *y1, long *y2, long *skip)
375 : {
376 1691 : *skip=0;
377 1691 : if (*y1==LONG_MAX)
378 : {
379 199 : if (*y2!=LONG_MAX)
380 : {
381 108 : if (*y2<0) *y2 += lA;
382 108 : if (*y2<0 || *y2==LONG_MAX || *y2>=lA)
383 0 : pari_err_DIM("_[..]");
384 108 : *skip=*y2;
385 : }
386 199 : *y1 = 1; *y2 = lA-1;
387 : }
388 1492 : else if (*y2==LONG_MAX) *y2 = *y1;
389 1691 : if (*y1<=0) *y1 += lA;
390 1691 : if (*y2<0) *y2 += lA;
391 1691 : if (*y1<=0 || *y1>*y2+1 || *y2>=lA) pari_err_DIM("_[..]");
392 1679 : return *y2 - *y1 + 2 - !!*skip;
393 : }
394 :
395 : static GEN
396 2229 : vecslice_i(GEN A, long t, long lB, long y1, long skip)
397 : {
398 2229 : GEN B = cgetg(lB, t);
399 : long i;
400 23220 : for (i=1; i<lB; i++, y1++)
401 : {
402 20991 : if (y1 == skip) { i--; continue; }
403 20883 : gel(B,i) = gcopy(gel(A,y1));
404 : }
405 2229 : return B;
406 : }
407 :
408 : static GEN
409 10 : rowslice_i(GEN A, long lB, long x1, long y1, long skip)
410 : {
411 10 : GEN B = cgetg(lB, t_VEC);
412 : long i;
413 55 : for (i=1; i<lB; i++, y1++)
414 : {
415 45 : if (y1 == skip) { i--; continue; }
416 40 : gel(B,i) = gcopy(gcoeff(A,x1,y1));
417 : }
418 10 : 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 20 : vecsmallslice_i(GEN A, long t, long lB, long y1, long skip)
436 : {
437 20 : GEN B = cgetg(lB, t);
438 : long i;
439 90 : for (i=1; i<lB; i++, y1++)
440 : {
441 70 : if (y1 == skip) { i--; continue; }
442 65 : B[i] = A[y1];
443 : }
444 20 : return B;
445 : }
446 : GEN
447 1383 : vecslice0(GEN A, long y1, long y2)
448 : {
449 1383 : long skip, lB, t = typ(A);
450 1383 : switch(t)
451 : {
452 1312 : case t_VEC: case t_COL:
453 1312 : lB = vecslice_parse_arg(lg(A), &y1, &y2, &skip);
454 1300 : return vecslice_i(A, t,lB,y1,skip);
455 20 : case t_VECSMALL:
456 20 : lB = vecslice_parse_arg(lg(A), &y1, &y2, &skip);
457 20 : return vecsmallslice_i(A, t,lB,y1,skip);
458 45 : case t_LIST:
459 45 : if (list_typ(A) == t_LIST_RAW)
460 : {
461 45 : GEN y, z = list_data(A);
462 45 : long l = z? lg(z): 1;
463 45 : lB = vecslice_parse_arg(l, &y1, &y2, &skip);
464 45 : y = mklist(); if (!z) return y;
465 45 : list_data(y) = vecslice_i(z, t_VEC,lB,y1,skip);
466 45 : return y;
467 : }
468 : default:
469 6 : pari_err_TYPE("_[_.._]",A);
470 : return NULL;/*LCOV_EXCL_LINE*/
471 : }
472 : }
473 :
474 : GEN
475 162 : matslice0(GEN A, long x1, long x2, long y1, long y2)
476 : {
477 : GEN B;
478 162 : long i, lB, lA = lg(A), rA, t, skip, rskip, rlB;
479 162 : long is_col = y1!=LONG_MAX && y2==LONG_MAX;
480 162 : long is_row = x1!=LONG_MAX && x2==LONG_MAX;
481 : GEN (*slice)(GEN A, long t, long lB, long y1, long skip);
482 162 : if (typ(A)!=t_MAT) pari_err_TYPE("_[_.._,_.._]",A);
483 162 : lB = vecslice_parse_arg(lA, &y1, &y2, &skip);
484 162 : if (is_col) return vecslice0(gel(A, y1), x1, x2);
485 152 : rA = lg(A)==1 ? 1: lgcols(A);
486 152 : rlB = vecslice_parse_arg(rA, &x1, &x2, &rskip);
487 152 : t = lg(A)==1 ? t_COL: typ(gel(A,1));
488 152 : if (is_row) return t == t_COL ? rowslice_i(A, lB, x1, y1, skip):
489 0 : rowsmallslice_i(A, lB, x1, y1, skip);
490 142 : slice = t == t_COL? &vecslice_i: &vecsmallslice_i;
491 :
492 142 : B = cgetg(lB, t_MAT);
493 1036 : for (i=1; i<lB; i++, y1++)
494 : {
495 894 : if (y1 == skip) { i--; continue; }
496 884 : gel(B,i) = slice(gel(A,y1),t,rlB, x1, rskip);
497 : }
498 142 : return B;
499 : }
500 :
501 : GEN
502 9119 : vecrange(GEN a, GEN b)
503 : {
504 : GEN y;
505 : long i, l;
506 9119 : if (typ(a)!=t_INT) pari_err_TYPE("[_.._]",a);
507 9113 : if (typ(b)!=t_INT) pari_err_TYPE("[_.._]",b);
508 9107 : if (cmpii(a,b)>0) return cgetg(1,t_VEC);
509 9102 : l = itos(subii(b,a))+1;
510 9102 : a = setloop(a);
511 9102 : y = cgetg(l+1, t_VEC);
512 18802132 : for (i=1; i<=l; a = incloop(a), i++) gel(y,i) = icopy(a);
513 9102 : 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 77 : 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 77 : switch(typ(A))
535 : {
536 11 : case t_LIST:
537 11 : z = list_data(A);
538 11 : l = z? lg(z): 1;
539 11 : if (list_typ(A)==t_LIST_MAP)
540 : {
541 6 : av = avma;
542 6 : return gc_GEN(av, mapselect_shallow(E, f, A));
543 : }
544 5 : break;
545 61 : case t_VEC: case t_COL: case t_MAT:
546 61 : l = lg(A);
547 61 : z = A;
548 61 : break;
549 5 : default:
550 5 : pari_err_TYPE("select",A);
551 : return NULL;/*LCOV_EXCL_LINE*/
552 : }
553 66 : v = cgetg(l, t_VECSMALL);
554 66 : av = avma;
555 66 : clone_lock(A);
556 657 : for (i = lv = 1; i < l; i++) {
557 591 : if (f(E, gel(z,i))) v[lv++] = i;
558 591 : set_avma(av);
559 : }
560 66 : clone_unlock_deep(A); fixlg(v, lv); return v;
561 : }
562 : static GEN
563 61 : extract_copy(GEN A, GEN v)
564 : {
565 61 : long i, l = lg(v);
566 61 : GEN B = cgetg(l, typ(A));
567 268 : for (i = 1; i < l; i++) gel(B,i) = gcopy(gel(A,v[i]));
568 61 : return B;
569 : }
570 : /* as genselect, but treat A [ t_VEC,t_COL, or t_MAT] as a t_VEC */
571 : GEN
572 0 : vecselect(void *E, long (*f)(void* E, GEN x), GEN A)
573 : {
574 : GEN v;
575 0 : clone_lock(A);
576 0 : v = genindexselect(E, f, A);
577 0 : A = extract_copy(A, v); settyp(A, t_VEC);
578 0 : clone_unlock_deep(A); return A;
579 : }
580 : GEN
581 78 : genselect(void *E, long (*f)(void* E, GEN x), GEN A)
582 : {
583 78 : pari_sp av = avma;
584 : GEN y, z, v;/* v left on stack for efficiency */
585 78 : clone_lock(A);
586 78 : switch(typ(A))
587 : {
588 28 : case t_LIST:
589 28 : z = list_data(A);
590 28 : if (!z) y = mklist();
591 : else
592 : {
593 28 : if (list_typ(A)==t_LIST_MAP)
594 : {
595 12 : long i, l = z? lg(z): 1, lv=1;
596 12 : GEN v1 = cgetg(l, t_COL);
597 12 : GEN v2 = cgetg(l, t_COL);
598 48 : for (i = lv = 1; i < l; i++) {
599 36 : if (f(E, gmael3(z,i,1,2)))
600 : {
601 18 : gel(v1,lv) = gmael3(z,i,1,1);
602 18 : gel(v2,lv) = gmael3(z,i,1,2);
603 18 : lv++;
604 : }
605 : }
606 12 : fixlg(v1, lv); fixlg(v2, lv); y = gtomap(mkmat2(v1,v2));
607 : }
608 : else
609 : {
610 : GEN B;
611 16 : y = cgetg(3, t_LIST);
612 16 : v = genindexselect(E, f, z);
613 16 : B = extract_copy(z, v);
614 16 : y[1] = lg(B)-1;
615 16 : list_data(y) = B;
616 : }
617 : }
618 28 : break;
619 45 : case t_VEC: case t_COL: case t_MAT:
620 45 : v = genindexselect(E, f, A);
621 45 : y = extract_copy(A, v);
622 45 : break;
623 5 : default:
624 5 : pari_err_TYPE("select",A);
625 : return NULL;/*LCOV_EXCL_LINE*/
626 : }
627 73 : clone_unlock_deep(A); return gc_upto(av, y);
628 : }
629 :
630 : static void
631 45332 : check_callgen1(GEN f, const char *s)
632 : {
633 45332 : if (typ(f) != t_CLOSURE || closure_is_variadic(f) || closure_arity(f) < 1)
634 0 : pari_err_TYPE(s, f);
635 45332 : }
636 :
637 : GEN
638 94 : select0(GEN f, GEN x, long flag)
639 : {
640 94 : check_callgen1(f, "select");
641 94 : switch(flag)
642 : {
643 78 : case 0: return genselect((void *) f, gp_callbool, x);
644 16 : case 1: return genindexselect((void *) f, gp_callbool, x);
645 0 : default: pari_err_FLAG("select");
646 : return NULL;/*LCOV_EXCL_LINE*/
647 : }
648 : }
649 :
650 : GEN
651 0 : parselect(GEN C, GEN D, long flag)
652 : {
653 : pari_sp av, av2;
654 0 : long lv, l = lg(D), i, pending = 0, workid;
655 : GEN V, done;
656 : struct pari_mt pt;
657 0 : check_callgen1(C, "parselect");
658 0 : if (!is_vec_t(typ(D))) pari_err_TYPE("parselect",D);
659 0 : V = cgetg(l, t_VECSMALL); av = avma;
660 0 : mt_queue_start_lim(&pt, C, l-1);
661 0 : av2 = avma;
662 0 : for (i=1; i<l || pending; i++)
663 : {
664 0 : mt_queue_submit(&pt, i, i<l? mkvec(gel(D,i)): NULL);
665 0 : done = mt_queue_get(&pt, &workid, &pending);
666 0 : if (done) V[workid] = !gequal0(done);
667 0 : set_avma(av2);
668 : }
669 0 : mt_queue_end(&pt);
670 0 : set_avma(av);
671 0 : for (lv=1, i=1; i<l; i++)
672 0 : if (V[i]) V[lv++]=i;
673 0 : fixlg(V, lv);
674 0 : return flag? V: extract_copy(D, V);
675 : }
676 :
677 : GEN
678 0 : veccatapply(void *E, GEN (*f)(void*, GEN), GEN x)
679 : {
680 0 : pari_sp av = avma;
681 0 : GEN v = vecapply(E, f, x);
682 0 : return lg(v) == 1? v: gc_GEN(av, shallowconcat1(v));
683 : }
684 :
685 : static GEN
686 6 : ser_apply(void *E, GEN (*f)(void*, GEN), GEN x)
687 24 : { pari_APPLY_ser(f(E, gel(x,i))); }
688 : static GEN
689 6 : RgX_apply(void *E, GEN (*f)(void*, GEN), GEN x)
690 24 : { pari_APPLY_pol(f(E, gel(x,i))); }
691 : static GEN
692 157486 : RgV_apply(void *E, GEN (*f)(void*, GEN), GEN x)
693 899796 : { pari_APPLY_same( f(E, gel(x,i)) ); }
694 : static GEN
695 36 : RgM_apply(void *E, GEN (*f)(void*, GEN), GEN x)
696 108 : { pari_APPLY_same( RgV_apply(E,f,gel(x,i)) ); }
697 :
698 : static GEN
699 54 : map_apply_i(void *E, GEN (*f)(void*, GEN), GEN x)
700 54 : { retmkvec2(mkvec2(gcopy(gmael(x,1,1)), f(E, gmael(x,1,2))),
701 : gcopy(gel(x, 2))); }
702 : static GEN
703 18 : map_apply(void *E, GEN (*f)(void* E, GEN x), GEN x)
704 72 : { pari_APPLY_same(map_apply_i(E, f, gel(x,i))); }
705 :
706 : /* as genapply, but treat A [ t_VEC,t_COL, or t_MAT] as a t_VEC */
707 : GEN
708 113208 : vecapply(void *E, GEN (*f)(void* E, GEN x), GEN x)
709 : {
710 : GEN y;
711 113208 : clone_lock(x); y = RgV_apply(E,f,x);
712 113208 : clone_unlock_deep(x); settyp(y, t_VEC); return y;
713 : }
714 : GEN
715 44278 : genapply(void *E, GEN (*f)(void* E, GEN x), GEN x)
716 : {
717 44278 : long tx = typ(x);
718 : GEN y, z;
719 :
720 44278 : if (is_scalar_t(tx)) return f(E, x);
721 44278 : clone_lock(x);
722 44278 : switch(tx) {
723 6 : case t_POL: y = RgX_apply(E,f,x); break;
724 6 : case t_SER:
725 6 : y = ser_isexactzero(x)? gcopy(x): ser_apply(E,f,x);
726 6 : break;
727 36 : case t_LIST:
728 : {
729 36 : long t = list_typ(x);
730 36 : z = list_data(x);
731 36 : if (!z)
732 6 : y = mklist_typ(t);
733 : else
734 : {
735 30 : y = cgetg(3, t_LIST);
736 30 : y[1] = evaltyp(t)|_evallg(lg(z)-1);
737 : switch(t)
738 : {
739 12 : case t_LIST_RAW: list_data(y) = RgV_apply(E,f,z); break;
740 18 : case t_LIST_MAP: list_data(y) = map_apply(E,f,z); break;
741 : }
742 : }
743 : }
744 36 : break;
745 36 : case t_MAT: y = RgM_apply(E, f, x); break;
746 44194 : case t_VEC: case t_COL: y = RgV_apply(E,f,x); break;
747 0 : default:
748 0 : pari_err_TYPE("apply",x);
749 : return NULL;/*LCOV_EXCL_LINE*/
750 : }
751 44278 : clone_unlock_deep(x); return y;
752 : }
753 :
754 : GEN
755 44278 : apply0(GEN f, GEN x)
756 : {
757 44278 : check_callgen1(f, "apply");
758 44278 : return genapply((void *) f, gp_call, x);
759 : }
760 :
761 : GEN
762 363 : vecselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
763 : GEN (*fun)(void* E, GEN x), GEN A)
764 : {
765 : GEN y;
766 363 : long i, l = lg(A), nb=1;
767 363 : clone_lock(A); y = cgetg(l, t_VEC);
768 18714387 : for (i=1; i<l; i++)
769 18714024 : if (pred(Epred, gel(A,i))) gel(y,nb++) = fun(Efun, gel(A,i));
770 363 : fixlg(y,nb); clone_unlock_deep(A); return y;
771 : }
772 :
773 : GEN
774 0 : veccatselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
775 : GEN (*fun)(void* E, GEN x), GEN A)
776 : {
777 0 : pari_sp av = avma;
778 0 : GEN v = vecselapply(Epred, pred, Efun, fun, A);
779 0 : return lg(v) == 1? v: gc_GEN(av, shallowconcat1(v));
780 : }
781 :
782 : /* suitable for gc_upto */
783 : GEN
784 10 : parapply_slice_worker(GEN x, GEN worker)
785 12510 : { pari_APPLY_same(closure_callgen1(worker, gel(x,i))); }
786 :
787 : /* B <- {A[i] : i = r (mod m)}, 1 <= r <= m */
788 : static void
789 10 : arithprogset(GEN B, GEN A, long r, long m)
790 : {
791 10 : long i, k, l = lg(A);
792 10 : if (typ(B)==t_VECSMALL)
793 0 : for (k = 1, i = r; i < l; i += m, k++) uel(B, k) = uel(A, i);
794 : else
795 12510 : for (k = 1, i = r; i < l; i += m, k++) gel(B, k) = gel(A, i);
796 10 : setlg(B, k);
797 10 : }
798 : GEN
799 10 : gen_parapply_slice(GEN worker, GEN D, long mmin)
800 : {
801 10 : long l, r, n = lg(D)-1, m = minss(mmin, n), pending = 0;
802 10 : GEN L = cgetg(n / m + 2, t_VEC), va = mkvec(L), V = cgetg_copy(D, &l);
803 : struct pari_mt pt;
804 10 : mt_queue_start_lim(&pt, worker, m);
805 20 : for (r = 1; r <= m || pending; r++)
806 : {
807 : long workid;
808 : GEN done;
809 10 : if (r <= m) arithprogset(L, D, r, m);
810 10 : mt_queue_submit(&pt, r, r <= m? va: NULL);
811 10 : done = mt_queue_get(&pt, &workid, &pending);
812 10 : if (done)
813 : {
814 10 : long j, k, J = lg(done)-1;
815 12510 : for (j = 1, k = workid; j <= J; j++, k +=m) gel(V, k) = gel(done, j);
816 : }
817 : }
818 10 : mt_queue_end(&pt); return V;
819 : }
820 :
821 : GEN
822 5754 : gen_parapply_slice_zv(GEN worker, GEN D, long mmin)
823 : {
824 5754 : long l, r, n = lg(D)-1, m = minss(mmin, n), pending = 0;
825 : struct pari_mt pt;
826 : GEN L, va, V;
827 5754 : if (m == 1) return closure_callgen1(worker, D);
828 0 : L = cgetg(n / m + 2, t_VECSMALL);
829 0 : va = mkvec(L); V = cgetg_copy(D, &l);
830 0 : mt_queue_start_lim(&pt, worker, m);
831 0 : for (r = 1; r <= m || pending; r++)
832 : {
833 : long workid;
834 : GEN done;
835 0 : if (r <= m) arithprogset(L, D, r, m);
836 0 : mt_queue_submit(&pt, r, r <= m? va: NULL);
837 0 : done = mt_queue_get(&pt, &workid, &pending);
838 0 : if (done)
839 : {
840 0 : long j, k, J = lg(done)-1;
841 0 : for (j = 1, k = workid; j <= J; j++, k +=m) V[k] = done[j];
842 : }
843 : }
844 0 : mt_queue_end(&pt); return V;
845 : }
846 :
847 : GEN
848 320862 : gen_parapply_percent(GEN worker, GEN D, long percent)
849 : {
850 320862 : long l = lg(D), i, pending = 0, cnt = 0, lper = -1, lcnt = 0;
851 : GEN W, V;
852 : struct pari_mt pt;
853 :
854 320862 : if (l == 1) return cgetg(1, typ(D));
855 294826 : W = cgetg(2, t_VEC);
856 294826 : V = cgetg(l, typ(D));
857 294826 : mt_queue_start_lim(&pt, worker, l-1);
858 2204056 : for (i = 1; i < l || pending; i++)
859 : {
860 : long workid;
861 : GEN done;
862 1909230 : if (i < l) gel(W,1) = gel(D,i);
863 1909230 : mt_queue_submit(&pt, i, i<l? W: NULL);
864 1909230 : done = mt_queue_get(&pt, &workid, &pending);
865 1909230 : if (done)
866 : {
867 1909230 : gel(V,workid) = done;
868 1909230 : if (percent && (++cnt)-lcnt>=percent)
869 : {
870 0 : long per = (long)(cnt*100./(l-1));
871 0 : lcnt = cnt;
872 0 : if (per > lper) { err_printf("%ld%% ",per); lper = per; }
873 : }
874 : }
875 : }
876 294826 : mt_queue_end(&pt); return V;
877 : }
878 :
879 : GEN
880 318372 : gen_parapply(GEN worker, GEN D)
881 318372 : { return gen_parapply_percent(worker, D, 0); }
882 :
883 : GEN
884 960 : parapply(GEN C, GEN D)
885 : {
886 960 : pari_sp av = avma;
887 960 : check_callgen1(C, "parapply");
888 960 : if (!is_vec_t(typ(D))) pari_err_TYPE("parapply",D);
889 960 : return gc_upto(av, gen_parapply(C, D));
890 : }
891 :
892 : GEN
893 24 : genfold(void *E, GEN (*f)(void* E, GEN x, GEN y), GEN x)
894 : {
895 24 : pari_sp av = avma;
896 : GEN z;
897 24 : long i, l = lg(x);
898 24 : if (!is_vec_t(typ(x))|| l==1 ) pari_err_TYPE("fold",x);
899 24 : clone_lock(x);
900 24 : z = gel(x,1);
901 102 : for (i=2; i<l; i++)
902 : {
903 78 : z = f(E,z,gel(x,i));
904 78 : if (gc_needed(av, 2))
905 : {
906 0 : if (DEBUGMEM>1) pari_warn(warnmem,"fold");
907 0 : z = gc_GEN(av, z);
908 : }
909 : }
910 24 : clone_unlock_deep(x);
911 24 : return gc_GEN(av, z);
912 : }
913 :
914 : GEN
915 24 : fold0(GEN f, GEN x)
916 : {
917 24 : if (typ(f) != t_CLOSURE || closure_arity(f) < 2) pari_err_TYPE("fold",f);
918 24 : return genfold((void *) f, gp_call2, x);
919 : }
920 : /*******************************************************************/
921 : /* */
922 : /* SCALAR-MATRIX OPERATIONS */
923 : /* */
924 : /*******************************************************************/
925 : GEN
926 322471 : gtomat(GEN x)
927 : {
928 : long lx, i;
929 : GEN y;
930 :
931 322471 : if (!x) return cgetg(1, t_MAT);
932 322432 : switch(typ(x))
933 : {
934 24 : case t_LIST:
935 24 : if (list_typ(x)==t_LIST_MAP)
936 12 : return maptomat(x);
937 12 : x = list_data(x);
938 12 : if (!x) return cgetg(1, t_MAT);
939 : /* fall through */
940 : case t_VEC: {
941 5767 : lx=lg(x); y=cgetg(lx,t_MAT);
942 5767 : if (lx == 1) break;
943 5767 : if (typ(gel(x,1)) == t_COL) {
944 3789 : long h = lgcols(x);
945 137316 : for (i=2; i<lx; i++) {
946 133527 : if (typ(gel(x,i)) != t_COL || lg(gel(x,i)) != h) break;
947 : }
948 3789 : if (i == lx) { /* matrix with h-1 rows */
949 3789 : y = cgetg(lx, t_MAT);
950 141105 : for (i=1 ; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
951 3789 : return y;
952 : }
953 : }
954 5928 : for (i=1; i<lx; i++) gel(y,i) = mkcolcopy(gel(x,i));
955 1978 : break;
956 : }
957 103878 : case t_COL:
958 103878 : lx = lg(x);
959 103878 : if (lx == 1) return cgetg(1, t_MAT);
960 103866 : if (typ(gel(x,1)) == t_VEC) {
961 6 : long j, h = lg(gel(x,1));
962 12 : for (i=2; i<lx; i++) {
963 6 : if (typ(gel(x,i)) != t_VEC || lg(gel(x,i)) != h) break;
964 : }
965 6 : if (i == lx) { /* matrix with h cols */
966 6 : y = cgetg(h, t_MAT);
967 24 : for (j=1 ; j<h; j++) {
968 18 : gel(y,j) = cgetg(lx, t_COL);
969 54 : for (i=1; i<lx; i++) gcoeff(y,i,j) = gcopy(gmael(x,i,j));
970 : }
971 6 : return y;
972 : }
973 : }
974 103860 : y = mkmatcopy(x); break;
975 211965 : case t_MAT:
976 211965 : y = gcopy(x); break;
977 6 : case t_QFB: {
978 : GEN b;
979 6 : y = cgetg(3,t_MAT); b = gmul2n(gel(x,2),-1);
980 6 : gel(y,1) = mkcol2(icopy(gel(x,1)), b);
981 6 : gel(y,2) = mkcol2(b, icopy(gel(x,3)));
982 6 : break;
983 : }
984 798 : default:
985 798 : y = cgetg(2,t_MAT); gel(y,1) = mkcolcopy(x);
986 798 : break;
987 : }
988 318607 : return y;
989 : }
990 :
991 : /* create the diagonal matrix, whose diagonal is given by x */
992 : GEN
993 1101569 : diagonal(GEN x)
994 : {
995 1101569 : long j, lx, tx = typ(x);
996 : GEN y;
997 :
998 1101569 : if (! is_matvec_t(tx)) return scalarmat(x,1);
999 1101563 : if (tx==t_MAT)
1000 : {
1001 12 : if (RgM_isdiagonal(x)) return gcopy(x);
1002 6 : pari_err_TYPE("diagonal",x);
1003 : }
1004 1101551 : lx=lg(x); y=cgetg(lx,t_MAT);
1005 3956236 : for (j=1; j<lx; j++)
1006 : {
1007 2854685 : gel(y,j) = zerocol(lx-1);
1008 2854685 : gcoeff(y,j,j) = gcopy(gel(x,j));
1009 : }
1010 1101551 : return y;
1011 : }
1012 : /* same, assuming x is a t_VEC/t_COL. Not memory clean. */
1013 : GEN
1014 297878 : diagonal_shallow(GEN x)
1015 : {
1016 297878 : long j, lx = lg(x);
1017 297878 : GEN y = cgetg(lx,t_MAT);
1018 :
1019 876537 : for (j=1; j<lx; j++)
1020 : {
1021 578659 : gel(y,j) = zerocol(lx-1);
1022 578659 : gcoeff(y,j,j) = gel(x,j);
1023 : }
1024 297878 : return y;
1025 : }
1026 :
1027 : GEN
1028 576 : zv_diagonal(GEN x)
1029 : {
1030 576 : long j, l = lg(x), n = l-1;
1031 576 : GEN y = cgetg(l,t_MAT);
1032 :
1033 2676 : for (j = 1; j < l; j++)
1034 : {
1035 2100 : gel(y,j) = zero_Flv(n);
1036 2100 : ucoeff(y,j,j) = uel(x,j);
1037 : }
1038 576 : return y;
1039 : }
1040 :
1041 : /* compute x*diagonal(d) */
1042 : GEN
1043 60 : matmuldiagonal(GEN x, GEN d)
1044 : {
1045 60 : if (typ(x)!=t_MAT) pari_err_TYPE("matmuldiagonal",x);
1046 60 : if (! is_vec_t(typ(d))) pari_err_TYPE("matmuldiagonal",d);
1047 60 : if (lg(d) != lg(x)) pari_err_OP("operation 'matmuldiagonal'", x,d);
1048 300 : pari_APPLY_same(RgC_Rg_mul(gel(x,i), gel(d,i)));
1049 : }
1050 :
1051 : /* compute A*B assuming the result is a diagonal matrix */
1052 : GEN
1053 6 : matmultodiagonal(GEN A, GEN B)
1054 : {
1055 6 : long i, j, hA, hB, lA = lg(A), lB = lg(B);
1056 6 : GEN y = matid(lB-1);
1057 :
1058 6 : if (typ(A) != t_MAT) pari_err_TYPE("matmultodiagonal",A);
1059 6 : if (typ(B) != t_MAT) pari_err_TYPE("matmultodiagonal",B);
1060 6 : hA = (lA == 1)? lB: lgcols(A);
1061 6 : hB = (lB == 1)? lA: lgcols(B);
1062 6 : if (lA != hB || lB != hA) pari_err_OP("operation 'matmultodiagonal'", A,B);
1063 48 : for (i=1; i<lB; i++)
1064 : {
1065 42 : GEN z = gen_0;
1066 336 : for (j=1; j<lA; j++) z = gadd(z, gmul(gcoeff(A,i,j),gcoeff(B,j,i)));
1067 42 : gcoeff(y,i,i) = z;
1068 : }
1069 6 : return y;
1070 : }
1071 :
1072 : /* [m[1,1], ..., m[l,l]], internal */
1073 : GEN
1074 836426 : RgM_diagonal_shallow(GEN m)
1075 : {
1076 836426 : long i, lx = lg(m);
1077 836426 : GEN y = cgetg(lx,t_VEC);
1078 2786852 : for (i=1; i<lx; i++) gel(y, i) = gcoeff(m,i,i);
1079 836426 : return y;
1080 : }
1081 :
1082 : /* same, public function */
1083 : GEN
1084 0 : RgM_diagonal(GEN m)
1085 : {
1086 0 : long i, lx = lg(m);
1087 0 : GEN y = cgetg(lx,t_VEC);
1088 0 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gcoeff(m,i,i));
1089 0 : return y;
1090 : }
|
