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 153203 : vecsum(GEN v)
34 : {
35 153203 : pari_sp av = avma;
36 153203 : long i, l, t = typ(v);
37 : GEN p;
38 153203 : if (t==t_LIST && list_typ(v)==t_LIST_RAW)
39 : {
40 14 : v = list_data(v);
41 14 : if (!v) return gen_0;
42 : }
43 153189 : else if (!is_vec_t(t))
44 7 : pari_err_TYPE("vecsum", v);
45 153189 : l = lg(v);
46 153189 : if (l == 1) return gen_0;
47 153182 : p = gel(v,1);
48 153182 : if (l == 2) return gcopy(p);
49 263469 : for (i=2; i<l; i++)
50 : {
51 179290 : p = gadd(p, gel(v,i));
52 179290 : if (gc_needed(av, 2))
53 : {
54 0 : if (DEBUGMEM>1) pari_warn(warnmem,"sum");
55 0 : p = gc_upto(av, p);
56 : }
57 : }
58 84179 : return gc_upto(av, p);
59 : }
60 :
61 : /*******************************************************************/
62 : /* */
63 : /* TRANSPOSE */
64 : /* */
65 : /*******************************************************************/
66 : /* A[x0,]~ */
67 : static GEN
68 27246475 : row_transpose(GEN A, long x0)
69 : {
70 27246475 : long i, lB = lg(A);
71 27246475 : GEN B = cgetg(lB, t_COL);
72 209813715 : for (i=1; i<lB; i++) gel(B, i) = gcoeff(A, x0, i);
73 27246475 : return B;
74 : }
75 : static GEN
76 18730 : row_transposecopy(GEN A, long x0)
77 : {
78 18730 : long i, lB = lg(A);
79 18730 : GEN B = cgetg(lB, t_COL);
80 162531 : for (i=1; i<lB; i++) gel(B, i) = gcopy(gcoeff(A, x0, i));
81 18730 : return B;
82 : }
83 :
84 : /* No copy*/
85 : GEN
86 7455779 : shallowtrans(GEN x)
87 : {
88 : long i, dx, lx;
89 : GEN y;
90 7455779 : switch(typ(x))
91 : {
92 1260 : case t_VEC: y = leafcopy(x); settyp(y,t_COL); break;
93 182 : case t_COL: y = leafcopy(x); settyp(y,t_VEC); break;
94 7454337 : case t_MAT:
95 7454337 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
96 7453000 : dx = lgcols(x); y = cgetg(dx,t_MAT);
97 34699475 : for (i = 1; i < dx; i++) gel(y,i) = row_transpose(x,i);
98 7453000 : break;
99 0 : default: pari_err_TYPE("shallowtrans",x);
100 : return NULL;/*LCOV_EXCL_LINE*/
101 : }
102 7454442 : return y;
103 : }
104 :
105 : GEN
106 43744 : gtrans(GEN x)
107 : {
108 : long i, dx, lx;
109 : GEN y;
110 43744 : switch(typ(x))
111 : {
112 36883 : case t_VEC: y = gcopy(x); settyp(y,t_COL); break;
113 4886 : case t_COL: y = gcopy(x); settyp(y,t_VEC); break;
114 1968 : case t_MAT:
115 1968 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
116 1961 : dx = lgcols(x); y = cgetg(dx,t_MAT);
117 20691 : for (i = 1; i < dx; i++) gel(y,i) = row_transposecopy(x,i);
118 1961 : break;
119 7 : default: pari_err_TYPE("gtrans",x);
120 : return NULL;/*LCOV_EXCL_LINE*/
121 : }
122 43730 : return y;
123 : }
124 :
125 : /*******************************************************************/
126 : /* */
127 : /* EXTRACTION */
128 : /* */
129 : /*******************************************************************/
130 :
131 : static long
132 182 : str_to_long(char *s, char **pt)
133 : {
134 182 : long a = atol(s);
135 182 : while (isspace((unsigned char)*s)) s++;
136 182 : if (*s == '-' || *s == '+') s++;
137 385 : while (isdigit((unsigned char)*s) || isspace((unsigned char)*s)) s++;
138 182 : *pt = s; return a;
139 : }
140 :
141 : static int
142 112 : get_range(char *s, long *a, long *b, long *cmpl, long lx)
143 : {
144 112 : long max = lx - 1;
145 :
146 112 : *a = 1; *b = max;
147 112 : if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0;
148 112 : if (!*s) return 0;
149 112 : if (*s != '.')
150 : {
151 105 : *a = str_to_long(s, &s);
152 105 : if (*a < 0) *a += lx;
153 105 : if (*a<1 || *a>max) return 0;
154 : }
155 112 : if (*s == '.')
156 : {
157 105 : s++; if (*s != '.') return 0;
158 105 : do s++; while (isspace((unsigned char)*s));
159 105 : if (*s)
160 : {
161 77 : *b = str_to_long(s, &s);
162 77 : if (*b < 0) *b += lx;
163 77 : if (*b<1 || *b>max || *s) return 0;
164 : }
165 98 : return 1;
166 : }
167 7 : if (*s) return 0;
168 7 : *b = *a; return 1;
169 : }
170 :
171 : static int
172 35 : extract_selector_ok(long lx, GEN L)
173 : {
174 : long i, l;
175 35 : switch (typ(L))
176 : {
177 7 : case t_INT: {
178 : long maxj;
179 7 : if (!signe(L)) return 1;
180 7 : l = lgefint(L)-1;
181 7 : maxj = BITS_IN_LONG - bfffo(*int_MSW(L));
182 7 : return ((l-2) * BITS_IN_LONG + maxj < lx);
183 : }
184 7 : case t_STR: {
185 : long first, last, cmpl;
186 7 : return get_range(GSTR(L), &first, &last, &cmpl, lx);
187 : }
188 14 : case t_VEC: case t_COL:
189 14 : l = lg(L);
190 28 : for (i=1; i<l; i++)
191 : {
192 21 : long j = itos(gel(L,i));
193 21 : if (j>=lx || j<=0) return 0;
194 : }
195 7 : return 1;
196 7 : case t_VECSMALL:
197 7 : l = lg(L);
198 21 : for (i=1; i<l; i++)
199 : {
200 14 : long j = L[i];
201 14 : if (j>=lx || j<=0) return 0;
202 : }
203 7 : return 1;
204 : }
205 0 : return 0;
206 : }
207 :
208 : GEN
209 13216 : shallowmatextract(GEN x, GEN l1, GEN l2)
210 : {
211 13216 : long i, j, n1 = lg(l1), n2 = lg(l2);
212 13216 : GEN M = cgetg(n2, t_MAT);
213 76839 : for(i=1; i < n2; i++)
214 : {
215 63623 : long ii = l2[i];
216 63623 : GEN C = cgetg(n1, t_COL);
217 1026844 : for (j=1; j < n1; j++)
218 : {
219 963221 : long jj = l1[j];
220 963221 : gel(C, j) = gmael(x, ii, jj);
221 : }
222 63623 : gel(M, i) = C;
223 : }
224 13216 : return M;
225 : }
226 :
227 : GEN
228 52687 : shallowextract(GEN x, GEN L)
229 : {
230 52687 : long i,j, tl = typ(L), tx = typ(x), lx = lg(x);
231 : GEN y;
232 :
233 52687 : switch(tx)
234 : {
235 52680 : case t_VEC:
236 : case t_COL:
237 : case t_MAT:
238 52680 : case t_VECSMALL: break;
239 7 : default: pari_err_TYPE("extract",x);
240 :
241 : }
242 52680 : if (tl==t_INT)
243 : { /* extract components of x as per the bits of mask L */
244 : long k, l, ix, iy, maxj;
245 : GEN Ld;
246 3318 : if (!signe(L)) return cgetg(1,tx);
247 3311 : y = new_chunk(lx);
248 3311 : l = lgefint(L)-1; ix = iy = 1;
249 3311 : maxj = BITS_IN_LONG - bfffo(*int_MSW(L));
250 3311 : if ((l-2) * BITS_IN_LONG + maxj >= lx)
251 7 : pari_err_TYPE("vecextract [mask too large]", L);
252 3706 : for (k = 2, Ld = int_LSW(L); k < l; k++, Ld = int_nextW(Ld))
253 : {
254 402 : ulong B = *Ld;
255 21906 : for (j = 0; j < BITS_IN_LONG; j++, B >>= 1, ix++)
256 21504 : if (B & 1) y[iy++] = x[ix];
257 : }
258 : { /* k = l */
259 3304 : ulong B = *Ld;
260 29302 : for (j = 0; j < maxj; j++, B >>= 1, ix++)
261 25998 : if (B & 1) y[iy++] = x[ix];
262 : }
263 3304 : y[0] = evaltyp(tx) | evallg(iy);
264 3304 : return y;
265 : }
266 49362 : if (tl==t_STR)
267 : {
268 105 : char *s = GSTR(L);
269 : long first, last, cmpl, d;
270 105 : if (! get_range(s, &first, &last, &cmpl, lx))
271 7 : pari_err_TYPE("vecextract [incorrect range]", L);
272 98 : if (lx == 1) return cgetg(1,tx);
273 98 : d = last - first;
274 98 : if (cmpl)
275 : {
276 21 : if (d >= 0)
277 : {
278 14 : y = cgetg(lx - (1+d),tx);
279 469 : for (j=1; j<first; j++) gel(y,j) = gel(x,j);
280 266 : for (i=last+1; i<lx; i++,j++) gel(y,j) = gel(x,i);
281 : }
282 : else
283 : {
284 7 : y = cgetg(lx - (1-d),tx);
285 14 : for (j=1,i=lx-1; i>first; i--,j++) gel(y,j) = gel(x,i);
286 14 : for (i=last-1; i>0; i--,j++) gel(y,j) = gel(x,i);
287 : }
288 : }
289 : else
290 : {
291 77 : if (d >= 0)
292 : {
293 35 : y = cgetg(d+2,tx);
294 112 : for (i=first,j=1; i<=last; i++,j++) gel(y,j) = gel(x,i);
295 : }
296 : else
297 : {
298 42 : y = cgetg(2-d,tx);
299 203 : for (i=first,j=1; i>=last; i--,j++) gel(y,j) = gel(x,i);
300 : }
301 : }
302 98 : return y;
303 : }
304 :
305 49257 : if (is_vec_t(tl))
306 : {
307 77 : long ll=lg(L); y=cgetg(ll,tx);
308 196 : for (i=1; i<ll; i++)
309 : {
310 133 : j = itos(gel(L,i));
311 133 : if (j<=0) pari_err_COMPONENT("vecextract","<=",gen_0,stoi(j));
312 126 : if (j>=lx) pari_err_COMPONENT("vecextract",">=",stoi(lx),stoi(j));
313 119 : gel(y,i) = gel(x,j);
314 : }
315 63 : return y;
316 : }
317 49180 : if (tl == t_VECSMALL)
318 : {
319 49173 : long ll=lg(L); y=cgetg(ll,tx);
320 212827 : for (i=1; i<ll; i++)
321 : {
322 163654 : j = L[i];
323 163654 : if (j<=0) pari_err_COMPONENT("vecextract","<=",gen_0,stoi(j));
324 163654 : if (j>=lx) pari_err_COMPONENT("vecextract",">=",stoi(lx),stoi(j));
325 163654 : gel(y,i) = gel(x,j);
326 : }
327 49173 : return y;
328 : }
329 7 : pari_err_TYPE("vecextract [mask]", L);
330 : return NULL; /* LCOV_EXCL_LINE */
331 : }
332 :
333 : /* does the component selector l select 0 component ? */
334 : static int
335 85 : select_0(GEN l)
336 : {
337 85 : switch(typ(l))
338 : {
339 14 : case t_INT:
340 14 : return (!signe(l));
341 50 : case t_VEC: case t_COL: case t_VECSMALL:
342 50 : return (lg(l) == 1);
343 : }
344 21 : return 0;
345 : }
346 :
347 : GEN
348 39914 : extract0(GEN x, GEN l1, GEN l2)
349 : {
350 39914 : pari_sp av = avma, av2;
351 : GEN y;
352 39914 : if (! l2)
353 : {
354 39829 : y = shallowextract(x, l1);
355 39787 : if (lg(y) == 1 || typ(y) == t_VECSMALL) return y;
356 39780 : av2 = avma;
357 39780 : y = gcopy(y);
358 : }
359 : else
360 : {
361 85 : if (typ(x) != t_MAT) pari_err_TYPE("extract",x);
362 85 : y = shallowextract(x,l2);
363 85 : if (select_0(l1)) { set_avma(av); return zeromat(0, lg(y)-1); }
364 71 : if (lg(y) == 1 && lg(x) > 1)
365 : {
366 35 : if (!extract_selector_ok(lgcols(x), l1))
367 7 : pari_err_TYPE("vecextract [incorrect mask]", l1);
368 28 : retgc_const(av, cgetg(1, t_MAT));
369 : }
370 36 : y = shallowextract(shallowtrans(y), l1);
371 36 : av2 = avma;
372 36 : y = gtrans(y);
373 : }
374 39816 : stackdummy(av, av2);
375 39816 : return y;
376 : }
377 :
378 : static long
379 2065 : vecslice_parse_arg(long lA, long *y1, long *y2, long *skip)
380 : {
381 2065 : *skip=0;
382 2065 : if (*y1==LONG_MAX)
383 : {
384 252 : if (*y2!=LONG_MAX)
385 : {
386 140 : if (*y2<0) *y2 += lA;
387 140 : if (*y2<0 || *y2==LONG_MAX || *y2>=lA)
388 0 : pari_err_DIM("_[..]");
389 140 : *skip=*y2;
390 : }
391 252 : *y1 = 1; *y2 = lA-1;
392 : }
393 1813 : else if (*y2==LONG_MAX) *y2 = *y1;
394 2065 : if (*y1<=0) *y1 += lA;
395 2065 : if (*y2<0) *y2 += lA;
396 2065 : if (*y1<=0 || *y1>*y2+1 || *y2>=lA) pari_err_DIM("_[..]");
397 2051 : return *y2 - *y1 + 2 - !!*skip;
398 : }
399 :
400 : static GEN
401 2702 : vecslice_i(GEN A, long t, long lB, long y1, long skip)
402 : {
403 2702 : GEN B = cgetg(lB, t);
404 : long i;
405 27986 : for (i=1; i<lB; i++, y1++)
406 : {
407 25284 : if (y1 == skip) { i--; continue; }
408 25144 : gel(B,i) = gcopy(gel(A,y1));
409 : }
410 2702 : return B;
411 : }
412 :
413 : static GEN
414 14 : rowslice_i(GEN A, long lB, long x1, long y1, long skip)
415 : {
416 14 : GEN B = cgetg(lB, t_VEC);
417 : long i;
418 77 : for (i=1; i<lB; i++, y1++)
419 : {
420 63 : if (y1 == skip) { i--; continue; }
421 56 : gel(B,i) = gcopy(gcoeff(A,x1,y1));
422 : }
423 14 : return B;
424 : }
425 :
426 : static GEN
427 0 : rowsmallslice_i(GEN A, long lB, long x1, long y1, long skip)
428 : {
429 0 : GEN B = cgetg(lB, t_VECSMALL);
430 : long i;
431 0 : for (i=1; i<lB; i++, y1++)
432 : {
433 0 : if (y1 == skip) { i--; continue; }
434 0 : B[i] = coeff(A,x1,y1);
435 : }
436 0 : return B;
437 : }
438 :
439 : static GEN
440 28 : vecsmallslice_i(GEN A, long t, long lB, long y1, long skip)
441 : {
442 28 : GEN B = cgetg(lB, t);
443 : long i;
444 126 : for (i=1; i<lB; i++, y1++)
445 : {
446 98 : if (y1 == skip) { i--; continue; }
447 91 : B[i] = A[y1];
448 : }
449 28 : return B;
450 : }
451 : GEN
452 1666 : vecslice0(GEN A, long y1, long y2)
453 : {
454 1666 : long skip, lB, t = typ(A);
455 1666 : switch(t)
456 : {
457 1568 : case t_VEC: case t_COL:
458 1568 : lB = vecslice_parse_arg(lg(A), &y1, &y2, &skip);
459 1554 : return vecslice_i(A, t,lB,y1,skip);
460 28 : case t_VECSMALL:
461 28 : lB = vecslice_parse_arg(lg(A), &y1, &y2, &skip);
462 28 : return vecsmallslice_i(A, t,lB,y1,skip);
463 63 : case t_LIST:
464 63 : if (list_typ(A) == t_LIST_RAW)
465 : {
466 63 : GEN y, z = list_data(A);
467 63 : long l = z? lg(z): 1;
468 63 : lB = vecslice_parse_arg(l, &y1, &y2, &skip);
469 63 : y = mklist(); if (!z) return y;
470 63 : list_data(y) = vecslice_i(z, t_VEC,lB,y1,skip);
471 63 : return y;
472 : }
473 : default:
474 7 : pari_err_TYPE("_[_.._]",A);
475 : return NULL;/*LCOV_EXCL_LINE*/
476 : }
477 : }
478 :
479 : GEN
480 210 : matslice0(GEN A, long x1, long x2, long y1, long y2)
481 : {
482 : GEN B;
483 210 : long i, lB, lA = lg(A), rA, t, skip, rskip, rlB;
484 210 : long is_col = y1!=LONG_MAX && y2==LONG_MAX;
485 210 : long is_row = x1!=LONG_MAX && x2==LONG_MAX;
486 : GEN (*slice)(GEN A, long t, long lB, long y1, long skip);
487 210 : if (typ(A)!=t_MAT) pari_err_TYPE("_[_.._,_.._]",A);
488 210 : lB = vecslice_parse_arg(lA, &y1, &y2, &skip);
489 210 : if (is_col) return vecslice0(gel(A, y1), x1, x2);
490 196 : rA = lg(A)==1 ? 1: lgcols(A);
491 196 : rlB = vecslice_parse_arg(rA, &x1, &x2, &rskip);
492 196 : t = lg(A)==1 ? t_COL: typ(gel(A,1));
493 196 : if (is_row) return t == t_COL ? rowslice_i(A, lB, x1, y1, skip):
494 0 : rowsmallslice_i(A, lB, x1, y1, skip);
495 182 : slice = t == t_COL? &vecslice_i: &vecsmallslice_i;
496 :
497 182 : B = cgetg(lB, t_MAT);
498 1281 : for (i=1; i<lB; i++, y1++)
499 : {
500 1099 : if (y1 == skip) { i--; continue; }
501 1085 : gel(B,i) = slice(gel(A,y1),t,rlB, x1, rskip);
502 : }
503 182 : return B;
504 : }
505 :
506 : GEN
507 10870 : vecrange(GEN a, GEN b)
508 : {
509 : GEN y;
510 : long i, l;
511 10870 : if (typ(a)!=t_INT) pari_err_TYPE("[_.._]",a);
512 10863 : if (typ(b)!=t_INT) pari_err_TYPE("[_.._]",b);
513 10856 : if (cmpii(a,b)>0) return cgetg(1,t_VEC);
514 10849 : l = itos(subii(b,a))+1;
515 10849 : a = setloop(a);
516 10849 : y = cgetg(l+1, t_VEC);
517 26388385 : for (i=1; i<=l; a = incloop(a), i++) gel(y,i) = icopy(a);
518 10849 : return y;
519 : }
520 :
521 : GEN
522 0 : vecrangess(long a, long b)
523 : {
524 : GEN y;
525 : long i, l;
526 0 : if (a>b) return cgetg(1,t_VEC);
527 0 : l = b-a+1;
528 0 : y = cgetg(l+1, t_VEC);
529 0 : for (i=1; i<=l; a++, i++) gel(y,i) = stoi(a);
530 0 : return y;
531 : }
532 :
533 : GEN
534 110 : genindexselect(void *E, long (*f)(void* E, GEN x), GEN A)
535 : {
536 : long l, i, lv;
537 : GEN v, z;
538 : pari_sp av;
539 110 : switch(typ(A))
540 : {
541 14 : case t_LIST:
542 14 : z = list_data(A);
543 14 : l = z? lg(z): 1;
544 14 : if (list_typ(A)==t_LIST_MAP)
545 : {
546 7 : av = avma;
547 7 : return gc_GEN(av, mapselect_shallow(E, f, A));
548 : }
549 7 : break;
550 89 : case t_VEC: case t_COL: case t_MAT:
551 89 : l = lg(A);
552 89 : z = A;
553 89 : break;
554 7 : default:
555 7 : pari_err_TYPE("select",A);
556 : return NULL;/*LCOV_EXCL_LINE*/
557 : }
558 96 : v = cgetg(l, t_VECSMALL);
559 96 : av = avma;
560 96 : clone_lock(A);
561 12878 : for (i = lv = 1; i < l; i++) {
562 12782 : if (f(E, gel(z,i))) v[lv++] = i;
563 12782 : set_avma(av);
564 : }
565 96 : clone_unlock_deep(A); fixlg(v, lv); return v;
566 : }
567 : static GEN
568 101 : extract_copy(GEN A, GEN v)
569 : {
570 101 : long i, l = lg(v);
571 101 : GEN B = cgetg(l, typ(A));
572 4609 : for (i = 1; i < l; i++) gel(B,i) = gcopy(gel(A,v[i]));
573 101 : return B;
574 : }
575 : /* as genselect, but treat A [ t_VEC,t_COL, or t_MAT] as a t_VEC */
576 : GEN
577 0 : vecselect(void *E, long (*f)(void* E, GEN x), GEN A)
578 : {
579 : GEN v;
580 0 : clone_lock(A);
581 0 : v = genindexselect(E, f, A);
582 0 : A = extract_copy(A, v); settyp(A, t_VEC);
583 0 : clone_unlock_deep(A); return A;
584 : }
585 : GEN
586 104 : genselect(void *E, long (*f)(void* E, GEN x), GEN A)
587 : {
588 104 : pari_sp av = avma;
589 : GEN y, z, v;/* v left on stack for efficiency */
590 104 : clone_lock(A);
591 104 : switch(typ(A))
592 : {
593 35 : case t_LIST:
594 35 : z = list_data(A);
595 35 : if (!z) y = mklist();
596 : else
597 : {
598 35 : if (list_typ(A)==t_LIST_MAP)
599 : {
600 14 : long i, l = z? lg(z): 1, lv=1;
601 14 : GEN v1 = cgetg(l, t_COL);
602 14 : GEN v2 = cgetg(l, t_COL);
603 56 : for (i = lv = 1; i < l; i++) {
604 42 : if (f(E, gmael3(z,i,1,2)))
605 : {
606 21 : gel(v1,lv) = gmael3(z,i,1,1);
607 21 : gel(v2,lv) = gmael3(z,i,1,2);
608 21 : lv++;
609 : }
610 : }
611 14 : fixlg(v1, lv); fixlg(v2, lv); y = gtomap(mkmat2(v1,v2));
612 : }
613 : else
614 : {
615 : GEN B;
616 21 : y = cgetg(3, t_LIST);
617 21 : v = genindexselect(E, f, z);
618 21 : B = extract_copy(z, v);
619 21 : y[1] = lg(B)-1;
620 21 : list_data(y) = B;
621 : }
622 : }
623 35 : break;
624 62 : case t_VEC: case t_COL: case t_MAT:
625 62 : v = genindexselect(E, f, A);
626 62 : y = extract_copy(A, v);
627 62 : break;
628 7 : default:
629 7 : pari_err_TYPE("select",A);
630 : return NULL;/*LCOV_EXCL_LINE*/
631 : }
632 97 : clone_unlock_deep(A); return gc_upto(av, y);
633 : }
634 :
635 : static void
636 55424 : check_callgen1(GEN f, const char *s)
637 : {
638 55424 : if (typ(f) != t_CLOSURE || closure_is_variadic(f) || closure_arity(f) < 1)
639 0 : pari_err_TYPE(s, f);
640 55424 : }
641 :
642 : GEN
643 131 : select0(GEN f, GEN x, long flag)
644 : {
645 131 : check_callgen1(f, "select");
646 131 : switch(flag)
647 : {
648 104 : case 0: return genselect((void *) f, gp_callbool, x);
649 27 : case 1: return genindexselect((void *) f, gp_callbool, x);
650 0 : default: pari_err_FLAG("select");
651 : return NULL;/*LCOV_EXCL_LINE*/
652 : }
653 : }
654 :
655 : GEN
656 30 : parselect(GEN C, GEN D, long flag)
657 : {
658 : pari_sp av, av2;
659 30 : long lv, l = lg(D), i, pending = 0, workid;
660 : GEN V, done;
661 : struct pari_mt pt;
662 30 : check_callgen1(C, "parselect");
663 30 : if (!is_vec_t(typ(D))) pari_err_TYPE("parselect",D);
664 30 : V = cgetg(l, t_VECSMALL); av = avma;
665 30 : mt_queue_start_lim(&pt, C, l-1);
666 30 : av2 = avma;
667 30204 : for (i=1; i<l || pending; i++)
668 : {
669 30174 : mt_queue_submit(&pt, i, i<l? mkvec(gel(D,i)): NULL);
670 30174 : done = mt_queue_get(&pt, &workid, &pending);
671 30174 : if (done) V[workid] = !gequal0(done);
672 30174 : set_avma(av2);
673 : }
674 30 : mt_queue_end(&pt);
675 30 : set_avma(av);
676 30054 : for (lv=1, i=1; i<l; i++)
677 30024 : if (V[i]) V[lv++]=i;
678 30 : fixlg(V, lv);
679 30 : return flag? V: extract_copy(D, V);
680 : }
681 :
682 : GEN
683 0 : veccatapply(void *E, GEN (*f)(void*, GEN), GEN x)
684 : {
685 0 : pari_sp av = avma;
686 0 : GEN v = vecapply(E, f, x);
687 0 : return lg(v) == 1? v: gc_GEN(av, shallowconcat1(v));
688 : }
689 :
690 : static GEN
691 7 : ser_apply(void *E, GEN (*f)(void*, GEN), GEN x)
692 28 : { pari_APPLY_ser(f(E, gel(x,i))); }
693 : static GEN
694 7 : RgX_apply(void *E, GEN (*f)(void*, GEN), GEN x)
695 28 : { pari_APPLY_pol(f(E, gel(x,i))); }
696 : static GEN
697 185858 : RgV_apply(void *E, GEN (*f)(void*, GEN), GEN x)
698 1063923 : { pari_APPLY_same( f(E, gel(x,i)) ); }
699 : static GEN
700 42 : RgM_apply(void *E, GEN (*f)(void*, GEN), GEN x)
701 126 : { pari_APPLY_same( RgV_apply(E,f,gel(x,i)) ); }
702 :
703 : static GEN
704 63 : map_apply_i(void *E, GEN (*f)(void*, GEN), GEN x)
705 63 : { retmkvec2(mkvec2(gcopy(gmael(x,1,1)), f(E, gmael(x,1,2))),
706 : gcopy(gel(x, 2))); }
707 : static GEN
708 21 : map_apply(void *E, GEN (*f)(void* E, GEN x), GEN x)
709 84 : { pari_APPLY_same(map_apply_i(E, f, gel(x,i))); }
710 :
711 : /* as genapply, but treat A [ t_VEC,t_COL, or t_MAT] as a t_VEC */
712 : GEN
713 132339 : vecapply(void *E, GEN (*f)(void* E, GEN x), GEN x)
714 : {
715 : GEN y;
716 132339 : clone_lock(x); y = RgV_apply(E,f,x);
717 132339 : clone_unlock_deep(x); settyp(y, t_VEC); return y;
718 : }
719 : GEN
720 53519 : genapply(void *E, GEN (*f)(void* E, GEN x), GEN x)
721 : {
722 53519 : long tx = typ(x);
723 : GEN y, z;
724 :
725 53519 : if (is_scalar_t(tx)) return f(E, x);
726 53519 : clone_lock(x);
727 53519 : switch(tx) {
728 7 : case t_POL: y = RgX_apply(E,f,x); break;
729 7 : case t_SER:
730 7 : y = ser_isexactzero(x)? gcopy(x): ser_apply(E,f,x);
731 7 : break;
732 42 : case t_LIST:
733 : {
734 42 : long t = list_typ(x);
735 42 : z = list_data(x);
736 42 : if (!z)
737 7 : y = mklist_typ(t);
738 : else
739 : {
740 35 : y = cgetg(3, t_LIST);
741 35 : y[1] = evaltyp(t)|_evallg(lg(z)-1);
742 35 : switch(t)
743 : {
744 14 : case t_LIST_RAW: list_data(y) = RgV_apply(E,f,z); break;
745 21 : case t_LIST_MAP: list_data(y) = map_apply(E,f,z); break;
746 : }
747 : }
748 : }
749 42 : break;
750 42 : case t_MAT: y = RgM_apply(E, f, x); break;
751 53421 : case t_VEC: case t_COL: y = RgV_apply(E,f,x); break;
752 0 : default:
753 0 : pari_err_TYPE("apply",x);
754 : return NULL;/*LCOV_EXCL_LINE*/
755 : }
756 53519 : clone_unlock_deep(x); return y;
757 : }
758 :
759 : GEN
760 53519 : apply0(GEN f, GEN x)
761 : {
762 53519 : check_callgen1(f, "apply");
763 53519 : return genapply((void *) f, gp_call, x);
764 : }
765 :
766 : GEN
767 469 : vecselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
768 : GEN (*fun)(void* E, GEN x), GEN A)
769 : {
770 : GEN y;
771 469 : long i, l = lg(A), nb=1;
772 469 : clone_lock(A); y = cgetg(l, t_VEC);
773 26166700 : for (i=1; i<l; i++)
774 26166231 : if (pred(Epred, gel(A,i))) gel(y,nb++) = fun(Efun, gel(A,i));
775 469 : fixlg(y,nb); clone_unlock_deep(A); return y;
776 : }
777 :
778 : GEN
779 70 : eqselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
780 : GEN (*fun)(void* E, GEN x), GEN A)
781 : {
782 : GEN y;
783 70 : if (!pred || pred(Epred, A))
784 : {
785 28 : y = cgetg(2, t_VEC);
786 28 : if (fun)
787 : {
788 28 : clone_lock(A);
789 28 : gel(y,1) = fun(Efun, A);
790 28 : clone_unlock_deep(A); return y;
791 : } else
792 0 : gel(y,1) = A;
793 : } else
794 42 : y = cgetg(1, t_VEC);
795 42 : return y;
796 : }
797 :
798 : GEN
799 0 : veccatselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
800 : GEN (*fun)(void* E, GEN x), GEN A)
801 : {
802 0 : pari_sp av = avma;
803 0 : GEN v = vecselapply(Epred, pred, Efun, fun, A);
804 0 : return lg(v) == 1? v: gc_GEN(av, shallowconcat1(v));
805 : }
806 :
807 : GEN
808 0 : eqcatselapply(void *Epred, long (*pred)(void* E, GEN x), void *Efun,
809 : GEN (*fun)(void* E, GEN x), GEN A)
810 : {
811 0 : pari_sp av = avma;
812 0 : GEN y = NULL;
813 0 : clone_lock(A);
814 0 : if (!pred || pred(Epred, A))
815 : {
816 0 : if (fun)
817 : {
818 0 : clone_lock(A);
819 0 : y = fun(Efun, A);
820 0 : clone_unlock_deep(A); return y;
821 : } else
822 0 : y = A;
823 : }
824 0 : clone_unlock_deep(A);
825 0 : if (!y) retgc_const(av, cgetg(1, t_VEC));
826 0 : return gc_GEN(av, y);
827 : }
828 :
829 : /* suitable for gc_upto */
830 : GEN
831 44 : parapply_slice_worker(GEN x, GEN worker)
832 17544 : { pari_APPLY_same(closure_callgen1(worker, gel(x,i))); }
833 :
834 : /* B <- {A[i] : i = r (mod m)}, 1 <= r <= m */
835 : static void
836 92704 : arithprogset(GEN B, GEN A, long r, long m)
837 : {
838 92704 : long i, k, l = lg(A);
839 92704 : if (typ(B)==t_VECSMALL)
840 197515 : for (k = 1, i = r; i < l; i += m, k++) uel(B, k) = uel(A, i);
841 : else
842 182127 : for (k = 1, i = r; i < l; i += m, k++) gel(B, k) = gel(A, i);
843 92704 : setlg(B, k);
844 92704 : }
845 : GEN
846 7750 : gen_parapply_slice(GEN worker, GEN D, long mmin)
847 : {
848 7750 : long l, r, n = lg(D)-1, m = minss(mmin, n), pending = 0;
849 7750 : GEN L = cgetg(n / m + 2, t_VEC), va = mkvec(L), V = cgetg_copy(D, &l);
850 : struct pari_mt pt;
851 7750 : mt_queue_start_lim(&pt, worker, m);
852 155880 : for (r = 1; r <= m || pending; r++)
853 : {
854 : long workid;
855 : GEN done;
856 148130 : if (r <= m) arithprogset(L, D, r, m);
857 148130 : mt_queue_submit(&pt, r, r <= m? va: NULL);
858 148130 : done = mt_queue_get(&pt, &workid, &pending);
859 148130 : if (done)
860 : {
861 77940 : long j, k, J = lg(done)-1;
862 182127 : for (j = 1, k = workid; j <= J; j++, k +=m) gel(V, k) = gel(done, j);
863 : }
864 : }
865 7750 : mt_queue_end(&pt); return V;
866 : }
867 :
868 : GEN
869 6825 : gen_parapply_slice_zv(GEN worker, GEN D, long mmin)
870 : {
871 6825 : long l, r, n = lg(D)-1, m = minss(mmin, n), pending = 0;
872 : struct pari_mt pt;
873 : GEN L, va, V;
874 6825 : if (m == 1) return closure_callgen1(worker, D);
875 974 : L = cgetg(n / m + 2, t_VECSMALL);
876 974 : va = mkvec(L); V = cgetg_copy(D, &l);
877 974 : mt_queue_start_lim(&pt, worker, m);
878 29528 : for (r = 1; r <= m || pending; r++)
879 : {
880 : long workid;
881 : GEN done;
882 28554 : if (r <= m) arithprogset(L, D, r, m);
883 28554 : mt_queue_submit(&pt, r, r <= m? va: NULL);
884 28554 : done = mt_queue_get(&pt, &workid, &pending);
885 28554 : if (done)
886 : {
887 14764 : long j, k, J = lg(done)-1;
888 197515 : for (j = 1, k = workid; j <= J; j++, k +=m) V[k] = done[j];
889 : }
890 : }
891 974 : mt_queue_end(&pt); return V;
892 : }
893 :
894 : GEN
895 382062 : gen_parapply_percent(GEN worker, GEN D, long percent)
896 : {
897 382062 : long l = lg(D), i, pending = 0, cnt = 0, lper = -1, lcnt = 0;
898 382062 : long W[] = {evaltyp(t_VEC) | _evallg(2), 0};
899 : GEN V;
900 : struct pari_mt pt;
901 :
902 382062 : if (l == 1) return cgetg(1, typ(D));
903 351612 : V = cgetg(l, typ(D));
904 351612 : mt_queue_start_lim(&pt, worker, l-1);
905 2678764 : for (i = 1; i < l || pending; i++)
906 : {
907 : long workid;
908 : GEN done;
909 2327152 : if (i < l) gel(W,1) = gel(D,i);
910 2327152 : mt_queue_submit(&pt, i, i<l? W: NULL);
911 2327152 : done = mt_queue_get(&pt, &workid, &pending);
912 2327152 : if (done)
913 : {
914 2270888 : gel(V,workid) = done;
915 2270888 : if (percent && (++cnt)-lcnt>=percent)
916 : {
917 0 : long per = (long)(cnt*100./(l-1));
918 0 : lcnt = cnt;
919 0 : if (per > lper) { err_printf("%ld%% ",per); lper = per; }
920 : }
921 : }
922 : }
923 351612 : mt_queue_end(&pt); return V;
924 : }
925 :
926 : GEN
927 372063 : gen_parapply(GEN worker, GEN D)
928 372063 : { return gen_parapply_percent(worker, D, 0); }
929 :
930 : GEN
931 1744 : parapply(GEN C, GEN D)
932 : {
933 1744 : pari_sp av = avma;
934 1744 : check_callgen1(C, "parapply");
935 1744 : if (!is_vec_t(typ(D))) pari_err_TYPE("parapply",D);
936 1744 : return gc_upto(av, gen_parapply(C, D));
937 : }
938 :
939 : GEN
940 0 : gen_parpairwiseop_percent(GEN worker, GEN D, long percent)
941 : {
942 0 : long l = lg(D)-1, i, pending = 0, cnt = 0, lper = -1, lcnt = 0, n = l>>1;
943 0 : long W[] = {evaltyp(t_VEC) | _evallg(3), 0, 0};
944 : GEN V;
945 : struct pari_mt pt;
946 :
947 0 : if (l == 0) return cgetg(1, typ(D));
948 0 : V = cgetg(((l+1)>>1)+1, t_VEC);
949 0 : if (odd(l))
950 0 : gel(V,n+1) = gcopy(gel(D, l));
951 0 : mt_queue_start_lim(&pt, worker, n);
952 0 : for (i = 1; i <= n || pending; i++)
953 : {
954 : long workid;
955 : GEN done;
956 0 : if (i <= n) { gel(W,1) = gel(D,2*i-1); gel(W,2) = gel(D,2*i); }
957 0 : mt_queue_submit(&pt, i, i<=n ? W: NULL);
958 0 : done = mt_queue_get(&pt, &workid, &pending);
959 0 : if (done)
960 : {
961 0 : gel(V,workid) = done;
962 0 : if (percent && (++cnt)-lcnt>=percent)
963 : {
964 0 : long per = (long)(cnt*100./n);
965 0 : lcnt = cnt;
966 0 : if (per > lper) { err_printf("%ld%% ",per); lper = per; }
967 : }
968 : }
969 : }
970 0 : mt_queue_end(&pt); return V;
971 : }
972 :
973 : GEN
974 0 : gen_parpairwiseop(GEN worker, GEN D)
975 0 : { return gen_parpairwiseop_percent(worker, D, 0); }
976 :
977 : GEN
978 28 : genfold(void *E, GEN (*f)(void* E, GEN x, GEN y), GEN x)
979 : {
980 28 : pari_sp av = avma;
981 : GEN z;
982 28 : long i, l = lg(x);
983 28 : if (!is_vec_t(typ(x))|| l==1 ) pari_err_TYPE("fold",x);
984 28 : clone_lock(x);
985 28 : z = gel(x,1);
986 119 : for (i=2; i<l; i++)
987 : {
988 91 : z = f(E,z,gel(x,i));
989 91 : if (gc_needed(av, 2))
990 : {
991 0 : if (DEBUGMEM>1) pari_warn(warnmem,"fold");
992 0 : z = gc_GEN(av, z);
993 : }
994 : }
995 28 : clone_unlock_deep(x);
996 28 : return gc_GEN(av, z);
997 : }
998 :
999 : GEN
1000 28 : fold0(GEN f, GEN x)
1001 : {
1002 28 : if (typ(f) != t_CLOSURE || closure_arity(f) < 2) pari_err_TYPE("fold",f);
1003 28 : return genfold((void *) f, gp_call2, x);
1004 : }
1005 : /*******************************************************************/
1006 : /* */
1007 : /* SCALAR-MATRIX OPERATIONS */
1008 : /* */
1009 : /*******************************************************************/
1010 : GEN
1011 377829 : gtomat(GEN x)
1012 : {
1013 : long lx, i;
1014 : GEN y;
1015 :
1016 377829 : if (!x) return cgetg(1, t_MAT);
1017 377780 : switch(typ(x))
1018 : {
1019 28 : case t_LIST:
1020 28 : if (list_typ(x)==t_LIST_MAP)
1021 14 : return maptomat(x);
1022 14 : x = list_data(x);
1023 14 : if (!x) return cgetg(1, t_MAT);
1024 : /* fall through */
1025 : case t_VEC: {
1026 6748 : lx=lg(x); y=cgetg(lx,t_MAT);
1027 6748 : if (lx == 1) break;
1028 6748 : if (typ(gel(x,1)) == t_COL) {
1029 4431 : long h = lgcols(x);
1030 160237 : for (i=2; i<lx; i++) {
1031 155806 : if (typ(gel(x,i)) != t_COL || lg(gel(x,i)) != h) break;
1032 : }
1033 4431 : if (i == lx) { /* matrix with h-1 rows */
1034 4431 : y = cgetg(lx, t_MAT);
1035 164668 : for (i=1 ; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
1036 4431 : return y;
1037 : }
1038 : }
1039 6944 : for (i=1; i<lx; i++) gel(y,i) = mkcolcopy(gel(x,i));
1040 2317 : break;
1041 : }
1042 121364 : case t_COL:
1043 121364 : lx = lg(x);
1044 121364 : if (lx == 1) return cgetg(1, t_MAT);
1045 121350 : if (typ(gel(x,1)) == t_VEC) {
1046 7 : long j, h = lg(gel(x,1));
1047 14 : for (i=2; i<lx; i++) {
1048 7 : if (typ(gel(x,i)) != t_VEC || lg(gel(x,i)) != h) break;
1049 : }
1050 7 : if (i == lx) { /* matrix with h cols */
1051 7 : y = cgetg(h, t_MAT);
1052 28 : for (j=1 ; j<h; j++) {
1053 21 : gel(y,j) = cgetg(lx, t_COL);
1054 63 : for (i=1; i<lx; i++) gcoeff(y,i,j) = gcopy(gmael(x,i,j));
1055 : }
1056 7 : return y;
1057 : }
1058 : }
1059 121343 : y = mkmatcopy(x); break;
1060 248667 : case t_MAT:
1061 248667 : y = gcopy(x); break;
1062 7 : case t_QFB: {
1063 : GEN b;
1064 7 : y = cgetg(3,t_MAT); b = gmul2n(gel(x,2),-1);
1065 7 : gel(y,1) = mkcol2(icopy(gel(x,1)), b);
1066 7 : gel(y,2) = mkcol2(b, icopy(gel(x,3)));
1067 7 : break;
1068 : }
1069 973 : default:
1070 973 : y = cgetg(2,t_MAT); gel(y,1) = mkcolcopy(x);
1071 973 : break;
1072 : }
1073 373307 : return y;
1074 : }
1075 :
1076 : /* create the diagonal matrix, whose diagonal is given by x */
1077 : GEN
1078 1348726 : diagonal(GEN x)
1079 : {
1080 1348726 : long j, lx, tx = typ(x);
1081 : GEN y;
1082 :
1083 1348726 : if (! is_matvec_t(tx)) return scalarmat(x,1);
1084 1348719 : if (tx==t_MAT)
1085 : {
1086 14 : if (RgM_isdiagonal(x)) return gcopy(x);
1087 7 : pari_err_TYPE("diagonal",x);
1088 : }
1089 1348705 : lx=lg(x); y=cgetg(lx,t_MAT);
1090 4816868 : for (j=1; j<lx; j++)
1091 : {
1092 3468163 : gel(y,j) = zerocol(lx-1);
1093 3468163 : gcoeff(y,j,j) = gcopy(gel(x,j));
1094 : }
1095 1348705 : return y;
1096 : }
1097 : /* same, assuming x is a t_VEC/t_COL. Not memory clean. */
1098 : GEN
1099 349629 : diagonal_shallow(GEN x)
1100 : {
1101 349629 : long j, lx = lg(x);
1102 349629 : GEN y = cgetg(lx,t_MAT);
1103 :
1104 1028608 : for (j=1; j<lx; j++)
1105 : {
1106 678979 : gel(y,j) = zerocol(lx-1);
1107 678979 : gcoeff(y,j,j) = gel(x,j);
1108 : }
1109 349629 : return y;
1110 : }
1111 :
1112 : GEN
1113 672 : zv_diagonal(GEN x)
1114 : {
1115 672 : long j, l = lg(x), n = l-1;
1116 672 : GEN y = cgetg(l,t_MAT);
1117 :
1118 3122 : for (j = 1; j < l; j++)
1119 : {
1120 2450 : gel(y,j) = zero_Flv(n);
1121 2450 : ucoeff(y,j,j) = uel(x,j);
1122 : }
1123 672 : return y;
1124 : }
1125 :
1126 : /* compute x*diagonal(d) */
1127 : GEN
1128 70 : matmuldiagonal(GEN x, GEN d)
1129 : {
1130 70 : if (typ(x)!=t_MAT) pari_err_TYPE("matmuldiagonal",x);
1131 70 : if (! is_vec_t(typ(d))) pari_err_TYPE("matmuldiagonal",d);
1132 70 : if (lg(d) != lg(x)) pari_err_OP("operation 'matmuldiagonal'", x,d);
1133 350 : pari_APPLY_same(RgC_Rg_mul(gel(x,i), gel(d,i)));
1134 : }
1135 :
1136 : /* compute A*B assuming the result is a diagonal matrix */
1137 : GEN
1138 7 : matmultodiagonal(GEN A, GEN B)
1139 : {
1140 7 : long i, j, hA, hB, lA = lg(A), lB = lg(B);
1141 7 : GEN y = matid(lB-1);
1142 :
1143 7 : if (typ(A) != t_MAT) pari_err_TYPE("matmultodiagonal",A);
1144 7 : if (typ(B) != t_MAT) pari_err_TYPE("matmultodiagonal",B);
1145 7 : hA = (lA == 1)? lB: lgcols(A);
1146 7 : hB = (lB == 1)? lA: lgcols(B);
1147 7 : if (lA != hB || lB != hA) pari_err_OP("operation 'matmultodiagonal'", A,B);
1148 56 : for (i=1; i<lB; i++)
1149 : {
1150 49 : GEN z = gen_0;
1151 392 : for (j=1; j<lA; j++) z = gadd(z, gmul(gcoeff(A,i,j),gcoeff(B,j,i)));
1152 49 : gcoeff(y,i,i) = z;
1153 : }
1154 7 : return y;
1155 : }
1156 :
1157 : /* [m[1,1], ..., m[l,l]], internal */
1158 : GEN
1159 1059022 : RgM_diagonal_shallow(GEN m)
1160 : {
1161 1059022 : long i, lx = lg(m);
1162 1059022 : GEN y = cgetg(lx,t_VEC);
1163 3495650 : for (i=1; i<lx; i++) gel(y, i) = gcoeff(m,i,i);
1164 1059022 : return y;
1165 : }
1166 :
1167 : /* same, public function */
1168 : GEN
1169 0 : RgM_diagonal(GEN m)
1170 : {
1171 0 : long i, lx = lg(m);
1172 0 : GEN y = cgetg(lx,t_VEC);
1173 0 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gcoeff(m,i,i));
1174 0 : return y;
1175 : }
|
