Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : GEN
19 1274696 : iferrpari(GEN a, GEN b, GEN c)
20 : {
21 : GEN res;
22 : struct pari_evalstate state;
23 1274696 : evalstate_save(&state);
24 1274699 : pari_CATCH(CATCH_ALL)
25 : {
26 : GEN E;
27 69758 : if (!b&&!c) return gnil;
28 34886 : E = evalstate_restore_err(&state);
29 34886 : if (c)
30 : {
31 292 : push_lex(E,c);
32 292 : res = closure_evalnobrk(c);
33 285 : pop_lex(1);
34 285 : if (gequal0(res))
35 7 : pari_err(0, E);
36 : }
37 34872 : if (!b) return gnil;
38 34872 : push_lex(E,b);
39 34872 : res = closure_evalgen(b);
40 34872 : pop_lex(1);
41 34872 : return res;
42 : } pari_TRY {
43 1274699 : res = closure_evalgen(a);
44 1239813 : } pari_ENDCATCH;
45 1239813 : return res;
46 : }
47 :
48 : /********************************************************************/
49 : /** **/
50 : /** ITERATIONS **/
51 : /** **/
52 : /********************************************************************/
53 :
54 : static void
55 5073900 : forparii(GEN a, GEN b, GEN code)
56 : {
57 5073900 : pari_sp av, av0 = avma;
58 : GEN aa;
59 5073900 : if (gcmp(b,a) < 0) return;
60 5013718 : if (typ(b) != t_INFINITY) b = gfloor(b);
61 5013721 : aa = a = setloop(a);
62 5013722 : av=avma;
63 5013722 : push_lex(a,code);
64 104196278 : while (gcmp(a,b) <= 0)
65 : {
66 99408051 : closure_evalvoid(code); if (loop_break()) break;
67 99192572 : a = get_lex(-1);
68 99185991 : if (a == aa)
69 : {
70 99186391 : a = incloop(a);
71 99180540 : if (a != aa) { set_lex(-1,a); aa = a; }
72 : }
73 : else
74 : { /* 'code' modified a ! Be careful (and slow) from now on */
75 24 : a = gaddgs(a,1);
76 28 : if (gc_needed(av,1))
77 : {
78 0 : if (DEBUGMEM>1) pari_warn(warnmem,"forparii");
79 0 : a = gerepileupto(av,a);
80 : }
81 28 : set_lex(-1,a);
82 : }
83 : }
84 5012938 : pop_lex(1); set_avma(av0);
85 : }
86 :
87 : void
88 5073907 : forpari(GEN a, GEN b, GEN code)
89 : {
90 5073907 : pari_sp ltop=avma, av;
91 5073907 : if (typ(a) == t_INT) { forparii(a,b,code); return; }
92 7 : b = gcopy(b); /* Kludge to work-around the a+(a=2) bug */
93 7 : av=avma;
94 7 : push_lex(a,code);
95 28 : while (gcmp(a,b) <= 0)
96 : {
97 21 : closure_evalvoid(code); if (loop_break()) break;
98 21 : a = get_lex(-1); a = gaddgs(a,1);
99 21 : if (gc_needed(av,1))
100 : {
101 0 : if (DEBUGMEM>1) pari_warn(warnmem,"forpari");
102 0 : a = gerepileupto(av,a);
103 : }
104 21 : set_lex(-1, a);
105 : }
106 7 : pop_lex(1); set_avma(ltop);
107 : }
108 :
109 : void
110 357 : foreachpari(GEN x, GEN code)
111 : {
112 : long i, l;
113 357 : switch(typ(x))
114 : {
115 14 : case t_LIST:
116 14 : x = list_data(x); /* FALL THROUGH */
117 14 : if (!x) return;
118 : case t_MAT: case t_VEC: case t_COL:
119 343 : break;
120 7 : default:
121 7 : pari_err_TYPE("foreach",x);
122 : return; /*LCOV_EXCL_LINE*/
123 : }
124 343 : clone_lock(x); l = lg(x);
125 343 : push_lex(gen_0,code);
126 2261 : for (i = 1; i < l; i++)
127 : {
128 1918 : set_lex(-1, gel(x,i));
129 1918 : closure_evalvoid(code); if (loop_break()) break;
130 : }
131 343 : pop_lex(1); clone_unlock_deep(x);
132 : }
133 :
134 : /* 0 < a <= b. Using small consecutive chunks to 1) limit memory use, 2) allow
135 : * cheap early abort */
136 : static int
137 63 : forfactoredpos(ulong a, ulong b, GEN code)
138 : {
139 63 : const ulong step = 1024;
140 63 : pari_sp av = avma;
141 : ulong x1;
142 6888 : for(x1 = a;; x1 += step, set_avma(av))
143 6825 : { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
144 6888 : ulong j, lv, x2 = (b >= 2*step && b - 2*step >= x1)? x1-1 + step: b;
145 6888 : GEN v = vecfactoru_i(x1, x2);
146 6888 : lv = lg(v);
147 7008421 : for (j = 1; j < lv; j++)
148 : {
149 7001547 : ulong n = x1-1 + j;
150 7001547 : set_lex(-1, mkvec2(utoipos(n), Flm_to_ZM(gel(v,j))));
151 7001547 : closure_evalvoid(code);
152 7001547 : if (loop_break()) return 1;
153 : }
154 6874 : if (x2 == b) break;
155 6825 : set_lex(-1, gen_0);
156 : }
157 49 : return 0;
158 : }
159 :
160 : /* vector of primes to squarefree factorization */
161 : static GEN
162 4255559 : zv_to_ZM(GEN v)
163 4255559 : { return mkmat2(zc_to_ZC(v), const_col(lg(v)-1,gen_1)); }
164 : /* vector of primes to negative squarefree factorization */
165 : static GEN
166 4255559 : zv_to_mZM(GEN v)
167 : {
168 4255559 : long i, l = lg(v);
169 4255559 : GEN w = cgetg(l+1, t_COL);
170 15388443 : gel(w,1) = gen_m1; for (i = 1; i < l; i++) gel(w,i+1) = utoipos(v[i]);
171 4255559 : return mkmat2(w, const_col(l,gen_1));
172 : }
173 : /* 0 <= a <= b. Using small consecutive chunks to 1) limit memory use, 2) allow
174 : * cheap early abort */
175 : static void
176 21 : forsquarefreepos(ulong a, ulong b, GEN code)
177 : {
178 21 : const ulong step = 1024;
179 21 : pari_sp av = avma;
180 : ulong x1;
181 6846 : for(x1 = a;; x1 += step, set_avma(av))
182 6825 : { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
183 6846 : ulong j, lv, x2 = (b >= 2*step && b - 2*step >= x1)? x1-1 + step: b;
184 6846 : GEN v = vecfactorsquarefreeu(x1, x2);
185 6846 : lv = lg(v);
186 7006958 : for (j = 1; j < lv; j++) if (gel(v,j))
187 : {
188 4255559 : ulong n = x1-1 + j;
189 4255559 : set_lex(-1, mkvec2(utoipos(n), zv_to_ZM(gel(v,j))));
190 4255559 : closure_evalvoid(code); if (loop_break()) return;
191 : }
192 6846 : if (x2 == b) break;
193 6825 : set_lex(-1, gen_0);
194 : }
195 : }
196 : /* 0 <= a <= b. Loop from -b, ... -a through squarefree integers */
197 : static void
198 21 : forsquarefreeneg(ulong a, ulong b, GEN code)
199 : {
200 21 : const ulong step = 1024;
201 21 : pari_sp av = avma;
202 : ulong x2;
203 6846 : for(x2 = b;; x2 -= step, set_avma(av))
204 6825 : { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
205 6846 : ulong j, x1 = (x2 >= 2*step && x2-2*step >= a)? x2+1 - step: a;
206 6846 : GEN v = vecfactorsquarefreeu(x1, x2);
207 7006958 : for (j = lg(v)-1; j > 0; j--) if (gel(v,j))
208 : {
209 4255559 : ulong n = x1-1 + j;
210 4255559 : set_lex(-1, mkvec2(utoineg(n), zv_to_mZM(gel(v,j))));
211 4255559 : closure_evalvoid(code); if (loop_break()) return;
212 : }
213 6846 : if (x1 == a) break;
214 6825 : set_lex(-1, gen_0);
215 : }
216 : }
217 : void
218 35 : forsquarefree(GEN a, GEN b, GEN code)
219 : {
220 35 : pari_sp av = avma;
221 : long s;
222 35 : if (typ(a) != t_INT) pari_err_TYPE("forsquarefree", a);
223 35 : if (typ(b) != t_INT) pari_err_TYPE("forsquarefree", b);
224 35 : if (cmpii(a,b) > 0) return;
225 35 : s = signe(a); push_lex(NULL,code);
226 35 : if (s < 0)
227 : {
228 21 : if (signe(b) <= 0)
229 14 : forsquarefreeneg(itou(b), itou(a), code);
230 : else
231 : {
232 7 : forsquarefreeneg(1, itou(a), code);
233 7 : forsquarefreepos(1, itou(b), code);
234 : }
235 : }
236 : else
237 14 : forsquarefreepos(itou(a), itou(b), code);
238 35 : pop_lex(1); set_avma(av);
239 : }
240 :
241 : /* convert factoru(n) to factor(-n); M pre-allocated factorization matrix
242 : * with (-1)^1 already set */
243 : static void
244 7001582 : Flm2negfact(GEN v, GEN M)
245 : {
246 7001582 : GEN p = gel(v,1), e = gel(v,2), P = gel(M,1), E = gel(M,2);
247 7001582 : long i, l = lg(p);
248 26980058 : for (i = 1; i < l; i++)
249 : {
250 19978476 : gel(P,i+1) = utoipos(p[i]);
251 19978476 : gel(E,i+1) = utoipos(e[i]);
252 : }
253 7001582 : setlg(P,l+1);
254 7001582 : setlg(E,l+1);
255 7001582 : }
256 : /* 0 < a <= b, from -b to -a */
257 : static int
258 84 : forfactoredneg(ulong a, ulong b, GEN code)
259 : {
260 84 : const ulong step = 1024;
261 : GEN P, E, M;
262 : pari_sp av;
263 : ulong x2;
264 :
265 84 : P = cgetg(18, t_COL); gel(P,1) = gen_m1;
266 84 : E = cgetg(18, t_COL); gel(E,1) = gen_1;
267 84 : M = mkmat2(P,E);
268 84 : av = avma;
269 6909 : for(x2 = b;; x2 -= step, set_avma(av))
270 6825 : { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
271 6909 : ulong j, x1 = (x2 >= 2*step && x2-2*step >= a)? x2+1 - step: a;
272 6909 : GEN v = vecfactoru_i(x1, x2);
273 7008470 : for (j = lg(v)-1; j; j--)
274 : { /* run backward: from factor(x1..x2) to factor(-x2..-x1) */
275 7001582 : ulong n = x1-1 + j;
276 7001582 : Flm2negfact(gel(v,j), M);
277 7001582 : set_lex(-1, mkvec2(utoineg(n), M));
278 7001582 : closure_evalvoid(code); if (loop_break()) return 1;
279 : }
280 6888 : if (x1 == a) break;
281 6825 : set_lex(-1, gen_0);
282 : }
283 63 : return 0;
284 : }
285 : static int
286 70 : eval0(GEN code)
287 : {
288 70 : pari_sp av = avma;
289 70 : set_lex(-1, mkvec2(gen_0, mkmat2(mkcol(gen_0),mkcol(gen_1))));
290 70 : closure_evalvoid(code); set_avma(av);
291 70 : return loop_break();
292 : }
293 : void
294 140 : forfactored(GEN a, GEN b, GEN code)
295 : {
296 140 : pari_sp av = avma;
297 140 : long sa, sb, stop = 0;
298 140 : if (typ(a) != t_INT) pari_err_TYPE("forfactored", a);
299 140 : if (typ(b) != t_INT) pari_err_TYPE("forfactored", b);
300 140 : if (cmpii(a,b) > 0) return;
301 133 : push_lex(NULL,code);
302 133 : sa = signe(a);
303 133 : sb = signe(b);
304 133 : if (sa < 0)
305 : {
306 84 : stop = forfactoredneg((sb < 0)? uel(b,2): 1UL, itou(a), code);
307 84 : if (!stop && sb >= 0) stop = eval0(code);
308 84 : if (!stop && sb > 0) forfactoredpos(1UL, b[2], code);
309 : }
310 : else
311 : {
312 49 : if (!sa) stop = eval0(code);
313 49 : if (!stop && sb) forfactoredpos(sa? uel(a,2): 1UL, itou(b), code);
314 : }
315 133 : pop_lex(1); set_avma(av);
316 : }
317 : void
318 1793741 : whilepari(GEN a, GEN b)
319 : {
320 1793741 : pari_sp av = avma;
321 : for(;;)
322 16907269 : {
323 18701010 : GEN res = closure_evalnobrk(a);
324 18701010 : if (gequal0(res)) break;
325 16907318 : set_avma(av);
326 16907318 : closure_evalvoid(b); if (loop_break()) break;
327 : }
328 1793741 : set_avma(av);
329 1793741 : }
330 :
331 : void
332 222074 : untilpari(GEN a, GEN b)
333 : {
334 222074 : pari_sp av = avma;
335 : for(;;)
336 1456773 : {
337 : GEN res;
338 1678847 : closure_evalvoid(b); if (loop_break()) break;
339 1678847 : res = closure_evalnobrk(a);
340 1678847 : if (!gequal0(res)) break;
341 1456773 : set_avma(av);
342 : }
343 222074 : set_avma(av);
344 222074 : }
345 :
346 28 : static int negcmp(GEN x, GEN y) { return gcmp(y,x); }
347 :
348 : void
349 1638 : forstep(GEN a, GEN b, GEN s, GEN code)
350 : {
351 : long ss, i;
352 1638 : pari_sp av, av0 = avma;
353 1638 : GEN v = NULL;
354 : int (*cmp)(GEN,GEN);
355 :
356 1638 : b = gcopy(b);
357 1638 : s = gcopy(s); av = avma;
358 1638 : switch(typ(s))
359 : {
360 14 : case t_VEC: case t_COL: ss = gsigne(vecsum(s)); v = s; break;
361 21 : case t_INTMOD:
362 21 : if (typ(a) != t_INT) a = gceil(a);
363 21 : a = addii(a, modii(subii(gel(s,2),a), gel(s,1)));
364 21 : s = gel(s,1);
365 1624 : default: ss = gsigne(s);
366 : }
367 1638 : if (!ss) pari_err_DOMAIN("forstep","step","=",gen_0,s);
368 1631 : cmp = (ss > 0)? &gcmp: &negcmp;
369 1631 : i = 0;
370 1631 : push_lex(a,code);
371 49756 : while (cmp(a,b) <= 0)
372 : {
373 48125 : closure_evalvoid(code); if (loop_break()) break;
374 48125 : if (v)
375 : {
376 98 : if (++i >= lg(v)) i = 1;
377 98 : s = gel(v,i);
378 : }
379 48125 : a = get_lex(-1); a = gadd(a,s);
380 :
381 48125 : if (gc_needed(av,1))
382 : {
383 0 : if (DEBUGMEM>1) pari_warn(warnmem,"forstep");
384 0 : a = gerepileupto(av,a);
385 : }
386 48125 : set_lex(-1,a);
387 : }
388 1631 : pop_lex(1); set_avma(av0);
389 1631 : }
390 :
391 : static void
392 28 : _fordiv(GEN a, GEN code, GEN (*D)(GEN))
393 : {
394 28 : pari_sp av = avma;
395 : long i, l;
396 28 : GEN t = D(a);
397 28 : push_lex(gen_0,code); l = lg(t);
398 231 : for (i=1; i<l; i++)
399 : {
400 203 : set_lex(-1,gel(t,i));
401 203 : closure_evalvoid(code); if (loop_break()) break;
402 : }
403 28 : pop_lex(1); set_avma(av);
404 28 : }
405 : void
406 14 : fordiv(GEN a, GEN code) { return _fordiv(a, code, &divisors); }
407 : void
408 14 : fordivfactored(GEN a, GEN code) { return _fordiv(a, code, &divisors_factored); }
409 :
410 : /* Embedded for loops:
411 : * fl = 0: execute ch (a), where a = (ai) runs through all n-uplets in
412 : * [m1,M1] x ... x [mn,Mn]
413 : * fl = 1: impose a1 <= ... <= an
414 : * fl = 2: a1 < ... < an
415 : */
416 : /* increment and return d->a [over integers]*/
417 : static GEN
418 183621 : _next_i(forvec_t *d)
419 : {
420 183621 : long i = d->n;
421 183621 : if (d->first) { d->first = 0; return (GEN)d->a; }
422 : for (;;) {
423 235720 : if (cmpii(d->a[i], d->M[i]) < 0) {
424 183259 : d->a[i] = incloop(d->a[i]);
425 183259 : return (GEN)d->a;
426 : }
427 52461 : d->a[i] = resetloop(d->a[i], d->m[i]);
428 52461 : if (--i <= 0) return NULL;
429 : }
430 : }
431 : /* increment and return d->a [generic]*/
432 : static GEN
433 63 : _next(forvec_t *d)
434 : {
435 63 : long i = d->n;
436 63 : if (d->first) { d->first = 0; return (GEN)d->a; }
437 : for (;;) {
438 98 : d->a[i] = gaddgs(d->a[i], 1);
439 98 : if (gcmp(d->a[i], d->M[i]) <= 0) return (GEN)d->a;
440 49 : d->a[i] = d->m[i];
441 49 : if (--i <= 0) return NULL;
442 : }
443 : }
444 :
445 : /* nondecreasing order [over integers] */
446 : static GEN
447 206 : _next_le_i(forvec_t *d)
448 : {
449 206 : long i = d->n;
450 206 : if (d->first) { d->first = 0; return (GEN)d->a; }
451 : for (;;) {
452 294 : if (cmpii(d->a[i], d->M[i]) < 0)
453 : {
454 152 : d->a[i] = incloop(d->a[i]);
455 : /* m[i] < a[i] <= M[i] <= M[i+1] */
456 233 : while (i < d->n)
457 : {
458 : GEN t;
459 81 : i++;
460 81 : if (cmpii(d->a[i-1], d->a[i]) <= 0) continue;
461 : /* a[i] < a[i-1] <= M[i-1] <= M[i] */
462 81 : t = d->a[i-1]; if (cmpii(t, d->m[i]) < 0) t = d->m[i];
463 81 : d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1],m[i])*/
464 : }
465 152 : return (GEN)d->a;
466 : }
467 142 : d->a[i] = resetloop(d->a[i], d->m[i]);
468 142 : if (--i <= 0) return NULL;
469 : }
470 : }
471 : /* nondecreasing order [generic] */
472 : static GEN
473 154 : _next_le(forvec_t *d)
474 : {
475 154 : long i = d->n;
476 154 : if (d->first) { d->first = 0; return (GEN)d->a; }
477 : for (;;) {
478 266 : d->a[i] = gaddgs(d->a[i], 1);
479 266 : if (gcmp(d->a[i], d->M[i]) <= 0)
480 : {
481 224 : while (i < d->n)
482 : {
483 : GEN c;
484 98 : i++;
485 98 : if (gcmp(d->a[i-1], d->a[i]) <= 0) continue;
486 : /* M[i] >= M[i-1] >= a[i-1] > a[i] */
487 98 : c = gceil(gsub(d->a[i-1], d->a[i]));
488 98 : d->a[i] = gadd(d->a[i], c);
489 : /* a[i-1] <= a[i] < M[i-1] + 1 => a[i] < M[i]+1 => a[i] <= M[i] */
490 : }
491 126 : return (GEN)d->a;
492 : }
493 140 : d->a[i] = d->m[i];
494 140 : if (--i <= 0) return NULL;
495 : }
496 : }
497 : /* strictly increasing order [over integers] */
498 : static GEN
499 1173574 : _next_lt_i(forvec_t *d)
500 : {
501 1173574 : long i = d->n;
502 1173574 : if (d->first) { d->first = 0; return (GEN)d->a; }
503 : for (;;) {
504 1290100 : if (cmpii(d->a[i], d->M[i]) < 0)
505 : {
506 1159954 : d->a[i] = incloop(d->a[i]);
507 : /* m[i] < a[i] <= M[i] < M[i+1] */
508 1276466 : while (i < d->n)
509 : {
510 : pari_sp av;
511 : GEN t;
512 116512 : i++;
513 116512 : if (cmpii(d->a[i-1], d->a[i]) < 0) continue;
514 116512 : av = avma;
515 : /* M[i] > M[i-1] >= a[i-1] */
516 116512 : t = addiu(d->a[i-1],1); if (cmpii(t, d->m[i]) < 0) t = d->m[i];
517 116512 : d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1]+1,m[i]) <= M[i]*/
518 116512 : set_avma(av);
519 : }
520 1159954 : return (GEN)d->a;
521 : }
522 130146 : d->a[i] = resetloop(d->a[i], d->m[i]);
523 130146 : if (--i <= 0) return NULL;
524 : }
525 : }
526 : /* strictly increasing order [generic] */
527 : static GEN
528 84 : _next_lt(forvec_t *d)
529 : {
530 84 : long i = d->n;
531 84 : if (d->first) { d->first = 0; return (GEN)d->a; }
532 : for (;;) {
533 133 : d->a[i] = gaddgs(d->a[i], 1);
534 133 : if (gcmp(d->a[i], d->M[i]) <= 0)
535 : {
536 91 : while (i < d->n)
537 : {
538 : GEN c;
539 35 : i++;
540 35 : if (gcmp(d->a[i-1], d->a[i]) < 0) continue;
541 : /* M[i] > M[i-1] >= a[i-1] >= a[i] */
542 35 : c = addiu(gfloor(gsub(d->a[i-1], d->a[i])), 1); /* > a[i-1] - a[i] */
543 35 : d->a[i] = gadd(d->a[i], c);
544 : /* a[i-1] < a[i] <= M[i-1] + 1 => a[i] < M[i]+1 => a[i] <= M[i] */
545 : }
546 56 : return (GEN)d->a;
547 : }
548 77 : d->a[i] = d->m[i];
549 77 : if (--i <= 0) return NULL;
550 : }
551 : }
552 : /* for forvec(v=[],) */
553 : static GEN
554 14 : _next_void(forvec_t *d)
555 : {
556 14 : if (d->first) { d->first = 0; return (GEN)d->a; }
557 7 : return NULL;
558 : }
559 :
560 : /* Initialize minima (m) and maxima (M); guarantee M[i] - m[i] integer and
561 : * if flag = 1: m[i-1] <= m[i] <= M[i] <= M[i+1]
562 : * if flag = 2: m[i-1] < m[i] <= M[i] < M[i+1],
563 : * for all i */
564 : int
565 7095 : forvec_init(forvec_t *d, GEN x, long flag)
566 : {
567 7095 : long i, tx = typ(x), l = lg(x), t = t_INT;
568 7095 : if (!is_vec_t(tx)) pari_err_TYPE("forvec [not a vector]", x);
569 7095 : d->first = 1;
570 7095 : d->n = l - 1;
571 7095 : d->a = (GEN*)cgetg(l,tx);
572 7095 : d->m = (GEN*)cgetg(l,tx);
573 7095 : d->M = (GEN*)cgetg(l,tx);
574 7095 : if (l == 1) { d->next = &_next_void; return 1; }
575 21502 : for (i = 1; i < l; i++)
576 : {
577 14442 : GEN a, e = gel(x,i), m = gel(e,1), M = gel(e,2);
578 14442 : tx = typ(e);
579 14442 : if (! is_vec_t(tx) || lg(e)!=3)
580 14 : pari_err_TYPE("forvec [expected vector not of type [min,MAX]]",e);
581 14428 : if (typ(m) != t_INT) t = t_REAL;
582 14428 : if (i > 1) switch(flag)
583 : {
584 62 : case 1: /* a >= m[i-1] - m */
585 62 : a = gceil(gsub(d->m[i-1], m));
586 62 : if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
587 62 : if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
588 62 : break;
589 6859 : case 2: /* a > m[i-1] - m */
590 6859 : a = gfloor(gsub(d->m[i-1], m));
591 6859 : if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
592 6859 : a = addiu(a, 1);
593 6859 : if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
594 6859 : break;
595 433 : default: m = gcopy(m);
596 433 : break;
597 : }
598 14428 : M = gadd(m, gfloor(gsub(M,m))); /* ensure M-m is an integer */
599 14421 : if (gcmp(m,M) > 0) { d->a = NULL; d->next = &_next; return 0; }
600 14414 : d->m[i] = m;
601 14414 : d->M[i] = M;
602 : }
603 7122 : if (flag == 1) for (i = l-2; i >= 1; i--)
604 : {
605 62 : GEN M = d->M[i], a = gfloor(gsub(d->M[i+1], M));
606 62 : if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
607 : /* M[i]+a <= M[i+1] */
608 62 : if (signe(a) < 0) d->M[i] = gadd(M, a);
609 : }
610 13871 : else if (flag == 2) for (i = l-2; i >= 1; i--)
611 : {
612 6852 : GEN M = d->M[i], a = gceil(gsub(d->M[i+1], M));
613 6852 : if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
614 6852 : a = subiu(a, 1);
615 : /* M[i]+a < M[i+1] */
616 6852 : if (signe(a) < 0) d->M[i] = gadd(M, a);
617 : }
618 7060 : if (t == t_INT) {
619 21327 : for (i = 1; i < l; i++) {
620 14302 : d->a[i] = setloop(d->m[i]);
621 14302 : if (typ(d->M[i]) != t_INT) d->M[i] = gfloor(d->M[i]);
622 : }
623 : } else {
624 140 : for (i = 1; i < l; i++) d->a[i] = d->m[i];
625 : }
626 7060 : switch(flag)
627 : {
628 188 : case 0: d->next = t==t_INT? &_next_i: &_next; break;
629 41 : case 1: d->next = t==t_INT? &_next_le_i: &_next_le; break;
630 6824 : case 2: d->next = t==t_INT? &_next_lt_i: &_next_lt; break;
631 7 : default: pari_err_FLAG("forvec");
632 : }
633 7053 : return 1;
634 : }
635 : GEN
636 1357716 : forvec_next(forvec_t *d) { return d->next(d); }
637 :
638 : void
639 7014 : forvec(GEN x, GEN code, long flag)
640 : {
641 7014 : pari_sp av = avma;
642 : forvec_t T;
643 : GEN v;
644 7014 : if (!forvec_init(&T, x, flag)) { set_avma(av); return; }
645 6979 : push_lex((GEN)T.a, code);
646 1357104 : while ((v = forvec_next(&T)))
647 : {
648 1350125 : closure_evalvoid(code);
649 1350125 : if (loop_break()) break;
650 : }
651 6979 : pop_lex(1); set_avma(av);
652 : }
653 :
654 : /********************************************************************/
655 : /** **/
656 : /** SUMS **/
657 : /** **/
658 : /********************************************************************/
659 :
660 : GEN
661 70217 : somme(GEN a, GEN b, GEN code, GEN x)
662 : {
663 70217 : pari_sp av, av0 = avma;
664 : GEN p1;
665 :
666 70217 : if (typ(a) != t_INT) pari_err_TYPE("sum",a);
667 70217 : if (!x) x = gen_0;
668 70217 : if (gcmp(b,a) < 0) return gcopy(x);
669 :
670 70217 : b = gfloor(b);
671 70217 : a = setloop(a);
672 70217 : av=avma;
673 70217 : push_lex(a,code);
674 : for(;;)
675 : {
676 1870211 : p1 = closure_evalnobrk(code);
677 1870211 : x=gadd(x,p1); if (cmpii(a,b) >= 0) break;
678 1799994 : a = incloop(a);
679 1799994 : if (gc_needed(av,1))
680 : {
681 0 : if (DEBUGMEM>1) pari_warn(warnmem,"sum");
682 0 : x = gerepileupto(av,x);
683 : }
684 1799994 : set_lex(-1,a);
685 : }
686 70217 : pop_lex(1); return gerepileupto(av0,x);
687 : }
688 :
689 : static GEN
690 28 : sum_init(GEN x0, GEN t)
691 : {
692 28 : long tp = typ(t);
693 : GEN x;
694 28 : if (is_vec_t(tp))
695 : {
696 7 : x = const_vec(lg(t)-1, x0);
697 7 : settyp(x, tp);
698 : }
699 : else
700 21 : x = x0;
701 28 : return x;
702 : }
703 :
704 : GEN
705 28 : suminf_bitprec(void *E, GEN (*eval)(void *, GEN), GEN a, long bit)
706 : {
707 28 : long fl = 0, G = bit + 1;
708 28 : pari_sp av0 = avma, av;
709 28 : GEN x = NULL, _1;
710 :
711 28 : if (typ(a) != t_INT) pari_err_TYPE("suminf",a);
712 28 : a = setloop(a); av = avma;
713 : for(;;)
714 15617 : {
715 15645 : GEN t = eval(E, a);
716 15645 : if (!x) _1 = x = sum_init(real_1_bit(bit), t);
717 :
718 15645 : x = gadd(x,t);
719 15645 : if (!gequal0(t) && gexpo(t) > gexpo(x)-G)
720 15449 : fl = 0;
721 196 : else if (++fl == 3)
722 28 : break;
723 15617 : a = incloop(a);
724 15617 : if (gc_needed(av,1))
725 : {
726 0 : if (DEBUGMEM>1) pari_warn(warnmem,"suminf");
727 0 : gerepileall(av,2, &x, &_1);
728 : }
729 : }
730 28 : return gerepileupto(av0, gsub(x, _1));
731 : }
732 : GEN
733 0 : suminf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
734 0 : { return suminf_bitprec(E, eval, a, prec2nbits(prec)); }
735 : GEN
736 28 : suminf0_bitprec(GEN a, GEN code, long bit)
737 28 : { EXPR_WRAP(code, suminf_bitprec(EXPR_ARG, a, bit)); }
738 :
739 : GEN
740 56 : sumdivexpr(GEN num, GEN code)
741 : {
742 56 : pari_sp av = avma;
743 56 : GEN y = gen_0, t = divisors(num);
744 56 : long i, l = lg(t);
745 :
746 56 : push_lex(gen_0, code);
747 9352 : for (i=1; i<l; i++)
748 : {
749 9296 : set_lex(-1,gel(t,i));
750 9296 : y = gadd(y, closure_evalnobrk(code));
751 : }
752 56 : pop_lex(1); return gerepileupto(av,y);
753 : }
754 :
755 : GEN
756 49 : sumdivmultexpr(void *D, GEN (*fun)(void*, GEN), GEN num)
757 : {
758 49 : pari_sp av = avma;
759 49 : GEN y = gen_1, P,E;
760 49 : int isint = divisors_init(num, &P,&E);
761 49 : long i, l = lg(P);
762 : GEN (*mul)(GEN,GEN);
763 :
764 49 : if (l == 1) return gc_const(av, gen_1);
765 49 : mul = isint? mulii: gmul;
766 224 : for (i=1; i<l; i++)
767 : {
768 175 : GEN p = gel(P,i), q = p, z = gen_1;
769 175 : long j, e = E[i];
770 581 : for (j = 1; j <= e; j++, q = mul(q, p))
771 : {
772 581 : z = gadd(z, fun(D, q));
773 581 : if (j == e) break;
774 : }
775 175 : y = gmul(y, z);
776 : }
777 49 : return gerepileupto(av,y);
778 : }
779 :
780 : GEN
781 49 : sumdivmultexpr0(GEN num, GEN code)
782 49 : { EXPR_WRAP(code, sumdivmultexpr(EXPR_ARG, num)) }
783 :
784 : /********************************************************************/
785 : /** **/
786 : /** PRODUCTS **/
787 : /** **/
788 : /********************************************************************/
789 :
790 : GEN
791 120694 : produit(GEN a, GEN b, GEN code, GEN x)
792 : {
793 120694 : pari_sp av, av0 = avma;
794 : GEN p1;
795 :
796 120694 : if (typ(a) != t_INT) pari_err_TYPE("prod",a);
797 120694 : if (!x) x = gen_1;
798 120694 : if (gcmp(b,a) < 0) return gcopy(x);
799 :
800 115416 : b = gfloor(b);
801 115416 : a = setloop(a);
802 115416 : av=avma;
803 115416 : push_lex(a,code);
804 : for(;;)
805 : {
806 349804 : p1 = closure_evalnobrk(code);
807 349804 : x = gmul(x,p1); if (cmpii(a,b) >= 0) break;
808 234388 : a = incloop(a);
809 234388 : if (gc_needed(av,1))
810 : {
811 0 : if (DEBUGMEM>1) pari_warn(warnmem,"prod");
812 0 : x = gerepileupto(av,x);
813 : }
814 234388 : set_lex(-1,a);
815 : }
816 115416 : pop_lex(1); return gerepileupto(av0,x);
817 : }
818 :
819 : GEN
820 14 : prodinf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
821 : {
822 14 : pari_sp av0 = avma, av;
823 : long fl,G;
824 14 : GEN p1,x = real_1(prec);
825 :
826 14 : if (typ(a) != t_INT) pari_err_TYPE("prodinf",a);
827 14 : a = setloop(a);
828 14 : av = avma;
829 14 : fl=0; G = -prec2nbits(prec)-5;
830 : for(;;)
831 : {
832 1897 : p1 = eval(E, a); if (gequal0(p1)) { x = p1; break; }
833 1897 : x = gmul(x,p1); a = incloop(a);
834 1897 : p1 = gsubgs(p1, 1);
835 1897 : if (gequal0(p1) || gexpo(p1) <= G) { if (++fl==3) break; } else fl=0;
836 1883 : if (gc_needed(av,1))
837 : {
838 0 : if (DEBUGMEM>1) pari_warn(warnmem,"prodinf");
839 0 : x = gerepileupto(av,x);
840 : }
841 : }
842 14 : return gerepilecopy(av0,x);
843 : }
844 : GEN
845 7 : prodinf1(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
846 : {
847 7 : pari_sp av0 = avma, av;
848 : long fl,G;
849 7 : GEN p1,p2,x = real_1(prec);
850 :
851 7 : if (typ(a) != t_INT) pari_err_TYPE("prodinf1",a);
852 7 : a = setloop(a);
853 7 : av = avma;
854 7 : fl=0; G = -prec2nbits(prec)-5;
855 : for(;;)
856 : {
857 952 : p2 = eval(E, a); p1 = gaddgs(p2,1);
858 952 : if (gequal0(p1)) { x = p1; break; }
859 952 : x = gmul(x,p1); a = incloop(a);
860 952 : if (gequal0(p2) || gexpo(p2) <= G) { if (++fl==3) break; } else fl=0;
861 945 : if (gc_needed(av,1))
862 : {
863 0 : if (DEBUGMEM>1) pari_warn(warnmem,"prodinf1");
864 0 : x = gerepileupto(av,x);
865 : }
866 : }
867 7 : return gerepilecopy(av0,x);
868 : }
869 : GEN
870 28 : prodinf0(GEN a, GEN code, long flag, long prec)
871 : {
872 28 : switch(flag)
873 : {
874 14 : case 0: EXPR_WRAP(code, prodinf (EXPR_ARG, a, prec));
875 7 : case 1: EXPR_WRAP(code, prodinf1(EXPR_ARG, a, prec));
876 : }
877 7 : pari_err_FLAG("prodinf");
878 : return NULL; /* LCOV_EXCL_LINE */
879 : }
880 :
881 : GEN
882 14 : prodeuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
883 : {
884 14 : pari_sp av, av0 = avma;
885 14 : GEN x = real_1(prec), prime;
886 : forprime_t T;
887 :
888 14 : av = avma;
889 14 : if (!forprime_init(&T, a,b)) return gc_const(av, x);
890 :
891 14 : av = avma;
892 8645 : while ( (prime = forprime_next(&T)) )
893 : {
894 8631 : x = gmul(x, eval(E, prime));
895 8631 : if (gc_needed(av,1))
896 : {
897 0 : if (DEBUGMEM>1) pari_warn(warnmem,"prodeuler");
898 0 : x = gerepilecopy(av, x);
899 : }
900 : }
901 14 : return gerepilecopy(av0,x);
902 : }
903 : GEN
904 14 : prodeuler0(GEN a, GEN b, GEN code, long prec)
905 14 : { EXPR_WRAP(code, prodeuler(EXPR_ARG, a, b, prec)); }
906 : GEN
907 133 : direuler0(GEN a, GEN b, GEN code, GEN c)
908 133 : { EXPR_WRAP(code, direuler(EXPR_ARG, a, b, c)); }
909 :
910 : /********************************************************************/
911 : /** **/
912 : /** VECTORS & MATRICES **/
913 : /** **/
914 : /********************************************************************/
915 :
916 : INLINE GEN
917 2766899 : copyupto(GEN z, GEN t)
918 : {
919 2766899 : if (is_universal_constant(z) || (z>(GEN)pari_mainstack->bot && z<=t))
920 2766892 : return z;
921 : else
922 7 : return gcopy(z);
923 : }
924 :
925 : GEN
926 115525 : vecexpr0(GEN vec, GEN code, GEN pred)
927 : {
928 115525 : switch(typ(vec))
929 : {
930 21 : case t_LIST:
931 : {
932 21 : if (list_typ(vec)==t_LIST_MAP)
933 7 : vec = mapdomain_shallow(vec);
934 : else
935 14 : vec = list_data(vec);
936 21 : if (!vec) return cgetg(1, t_VEC);
937 14 : break;
938 : }
939 7 : case t_VECSMALL:
940 7 : vec = vecsmall_to_vec(vec);
941 7 : break;
942 115497 : case t_VEC: case t_COL: case t_MAT: break;
943 0 : default: pari_err_TYPE("[_|_<-_,_]",vec);
944 : }
945 115518 : if (pred && code)
946 448 : EXPR_WRAP(code,vecselapply((void*)pred,&gp_evalbool,EXPR_ARGUPTO,vec))
947 115070 : else if (code)
948 115070 : EXPR_WRAP(code,vecapply(EXPR_ARGUPTO,vec))
949 : else
950 0 : EXPR_WRAP(pred,vecselect(EXPR_ARGBOOL,vec))
951 : }
952 :
953 : GEN
954 175 : vecexpr1(GEN vec, GEN code, GEN pred)
955 : {
956 175 : GEN v = vecexpr0(vec, code, pred);
957 175 : return lg(v) == 1? v: shallowconcat1(v);
958 : }
959 :
960 : GEN
961 2363591 : vecteur(GEN nmax, GEN code)
962 : {
963 : GEN y, c;
964 2363591 : long i, m = gtos(nmax);
965 :
966 2363591 : if (m < 0) pari_err_DOMAIN("vector", "dimension", "<", gen_0, stoi(m));
967 2363577 : if (!code) return zerovec(m);
968 17443 : c = cgetipos(3); /* left on stack */
969 17443 : y = cgetg(m+1,t_VEC); push_lex(c, code);
970 853749 : for (i=1; i<=m; i++)
971 : {
972 836320 : c[2] = i;
973 836320 : gel(y,i) = copyupto(closure_evalnobrk(code), y);
974 836306 : set_lex(-1,c);
975 : }
976 17429 : pop_lex(1); return y;
977 : }
978 :
979 : GEN
980 770 : vecteursmall(GEN nmax, GEN code)
981 : {
982 : pari_sp av;
983 : GEN y, c;
984 770 : long i, m = gtos(nmax);
985 :
986 770 : if (m < 0) pari_err_DOMAIN("vectorsmall", "dimension", "<", gen_0, stoi(m));
987 763 : if (!code) return zero_zv(m);
988 742 : c = cgetipos(3); /* left on stack */
989 742 : y = cgetg(m+1,t_VECSMALL); push_lex(c,code);
990 742 : av = avma;
991 16002 : for (i = 1; i <= m; i++)
992 : {
993 15267 : c[2] = i;
994 15267 : y[i] = gtos(closure_evalnobrk(code));
995 15260 : set_avma(av);
996 15260 : set_lex(-1,c);
997 : }
998 735 : pop_lex(1); return y;
999 : }
1000 :
1001 : GEN
1002 504 : vvecteur(GEN nmax, GEN n)
1003 : {
1004 504 : GEN y = vecteur(nmax,n);
1005 497 : settyp(y,t_COL); return y;
1006 : }
1007 :
1008 : GEN
1009 160727 : matrice(GEN nlig, GEN ncol, GEN code)
1010 : {
1011 : GEN c1, c2, y;
1012 : long i, m, n;
1013 :
1014 160727 : n = gtos(nlig);
1015 160727 : m = ncol? gtos(ncol): n;
1016 160727 : if (m < 0) pari_err_DOMAIN("matrix", "nbcols", "<", gen_0, stoi(m));
1017 160720 : if (n < 0) pari_err_DOMAIN("matrix", "nbrows", "<", gen_0, stoi(n));
1018 160713 : if (!m) return cgetg(1,t_MAT);
1019 160643 : if (!code || !n) return zeromatcopy(n, m);
1020 158438 : c1 = cgetipos(3); push_lex(c1,code);
1021 158438 : c2 = cgetipos(3); push_lex(c2,NULL); /* c1,c2 left on stack */
1022 158438 : y = cgetg(m+1,t_MAT);
1023 595721 : for (i = 1; i <= m; i++)
1024 : {
1025 437283 : GEN z = cgetg(n+1,t_COL);
1026 : long j;
1027 437283 : c2[2] = i; gel(y,i) = z;
1028 2367869 : for (j = 1; j <= n; j++)
1029 : {
1030 1930586 : c1[2] = j;
1031 1930586 : gel(z,j) = copyupto(closure_evalnobrk(code), y);
1032 1930586 : set_lex(-2,c1);
1033 1930586 : set_lex(-1,c2);
1034 : }
1035 : }
1036 158438 : pop_lex(2); return y;
1037 : }
1038 :
1039 : /********************************************************************/
1040 : /** **/
1041 : /** SUMMING SERIES **/
1042 : /** **/
1043 : /********************************************************************/
1044 : /* h = (2+2x)g'- g; g has t_INT coeffs */
1045 : static GEN
1046 1295 : delt(GEN g, long n)
1047 : {
1048 1295 : GEN h = cgetg(n+3,t_POL);
1049 : long k;
1050 1295 : h[1] = g[1];
1051 1295 : gel(h,2) = gel(g,2);
1052 359954 : for (k=1; k<n; k++)
1053 358659 : gel(h,k+2) = addii(mului(k+k+1,gel(g,k+2)), mului(k<<1,gel(g,k+1)));
1054 1295 : gel(h,n+2) = mului(n<<1, gel(g,n+1)); return h;
1055 : }
1056 :
1057 : #ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */
1058 : #pragma optimize("g",off)
1059 : #endif
1060 : /* P = polzagier(n,m)(-X), unnormalized (P(0) != 1) */
1061 : static GEN
1062 84 : polzag1(long n, long m)
1063 : {
1064 84 : long d = n - m, i, k, d2, r, D;
1065 84 : pari_sp av = avma;
1066 : GEN g, T;
1067 :
1068 84 : if (d <= 0 || m < 0) return pol_0(0);
1069 77 : d2 = d << 1; r = (m+1) >> 1, D = (d+1) >> 1;
1070 77 : g = cgetg(d+2, t_POL);
1071 77 : g[1] = evalsigne(1)|evalvarn(0);
1072 77 : T = cgetg(d+1,t_VEC);
1073 : /* T[k+1] = binomial(2d,2k+1), 0 <= k < d */
1074 77 : gel(T,1) = utoipos(d2);
1075 1344 : for (k = 1; k < D; k++)
1076 : {
1077 1267 : long k2 = k<<1;
1078 1267 : gel(T,k+1) = diviiexact(mulii(gel(T,k), muluu(d2-k2+1, d2-k2)),
1079 1267 : muluu(k2,k2+1));
1080 : }
1081 1365 : for (; k < d; k++) gel(T,k+1) = gel(T,d-k);
1082 77 : gel(g,2) = gel(T,d); /* binomial(2d, 2(d-1)+1) */
1083 2632 : for (i = 1; i < d; i++)
1084 : {
1085 2555 : pari_sp av2 = avma;
1086 2555 : GEN s, t = gel(T,d-i); /* binomial(2d, 2(d-1-i)+1) */
1087 2555 : s = t;
1088 180635 : for (k = d-i; k < d; k++)
1089 : {
1090 178080 : long k2 = k<<1;
1091 178080 : t = diviiexact(mulii(t, muluu(d2-k2+1, d-k)), muluu(k2+1,k-(d-i)+1));
1092 178080 : s = addii(s, t);
1093 : }
1094 : /* g_i = sum_{d-1-i <= k < d}, binomial(2*d, 2*k+1)*binomial(k,d-1-i) */
1095 2555 : gel(g,i+2) = gerepileuptoint(av2, s);
1096 : }
1097 : /* sum_{0 <= i < d} g_i x^i * (x+x^2)^r */
1098 77 : g = RgX_mulXn(gmul(g, gpowgs(deg1pol(gen_1,gen_1,0),r)), r);
1099 77 : if (!odd(m)) g = delt(g, n);
1100 1337 : for (i = 1; i <= r; i++)
1101 : {
1102 1260 : g = delt(ZX_deriv(g), n);
1103 1260 : if (gc_needed(av,4))
1104 : {
1105 0 : if (DEBUGMEM>1) pari_warn(warnmem,"polzag, i = %ld/%ld", i,r);
1106 0 : g = gerepilecopy(av, g);
1107 : }
1108 : }
1109 77 : return g;
1110 : }
1111 : GEN
1112 35 : polzag(long n, long m)
1113 : {
1114 35 : pari_sp av = avma;
1115 35 : GEN g = polzag1(n,m);
1116 35 : if (lg(g) == 2) return g;
1117 28 : g = ZX_z_unscale(polzag1(n,m), -1);
1118 28 : return gerepileupto(av, RgX_Rg_div(g,gel(g,2)));
1119 : }
1120 :
1121 : /*0.39322 > 1/log_2(3+sqrt(8))*/
1122 : static ulong
1123 154 : sumalt_N(long prec)
1124 154 : { return (ulong)(0.39322*(prec2nbits(prec) + 7)); }
1125 :
1126 : GEN
1127 84 : sumalt(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
1128 : {
1129 : ulong k, N;
1130 84 : pari_sp av = avma, av2;
1131 : GEN s, az, c, d;
1132 :
1133 84 : if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
1134 84 : N = sumalt_N(prec);
1135 84 : d = powru(addsr(3, sqrtr(utor(8,prec))), N);
1136 84 : d = shiftr(addrr(d, invr(d)),-1);
1137 84 : a = setloop(a);
1138 84 : az = gen_m1; c = d;
1139 84 : s = gen_0;
1140 84 : av2 = avma;
1141 84 : for (k=0; ; k++) /* k < N */
1142 : {
1143 10752 : c = addir(az,c); s = gadd(s, gmul(c, eval(E, a)));
1144 10752 : if (k==N-1) break;
1145 10668 : az = diviuuexact(muluui((N-k)<<1,N+k,az), k+1, (k<<1)+1);
1146 10668 : a = incloop(a); /* in place! */
1147 10668 : if (gc_needed(av,4))
1148 : {
1149 0 : if (DEBUGMEM>1) pari_warn(warnmem,"sumalt, k = %ld/%ld", k,N-1);
1150 0 : gerepileall(av2, 3, &az,&c,&s);
1151 : }
1152 : }
1153 84 : return gerepileupto(av, gdiv(s,d));
1154 : }
1155 :
1156 : GEN
1157 7 : sumalt2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
1158 : {
1159 : long k, N;
1160 7 : pari_sp av = avma, av2;
1161 : GEN s, dn, pol;
1162 :
1163 7 : if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
1164 7 : N = (long)(0.307073*(prec2nbits(prec) + 5)); /*0.307073 > 1/log_2(\beta_B)*/
1165 7 : pol = ZX_div_by_X_1(polzag1(N,N>>1), &dn);
1166 7 : a = setloop(a);
1167 7 : N = degpol(pol);
1168 7 : s = gen_0;
1169 7 : av2 = avma;
1170 280 : for (k=0; k<=N; k++)
1171 : {
1172 280 : GEN t = itor(gel(pol,k+2), prec+EXTRAPRECWORD);
1173 280 : s = gadd(s, gmul(t, eval(E, a)));
1174 280 : if (k == N) break;
1175 273 : a = incloop(a); /* in place! */
1176 273 : if (gc_needed(av,4))
1177 : {
1178 0 : if (DEBUGMEM>1) pari_warn(warnmem,"sumalt2, k = %ld/%ld", k,N-1);
1179 0 : s = gerepileupto(av2, s);
1180 : }
1181 : }
1182 7 : return gerepileupto(av, gdiv(s,dn));
1183 : }
1184 :
1185 : GEN
1186 28 : sumalt0(GEN a, GEN code, long flag, long prec)
1187 : {
1188 28 : switch(flag)
1189 : {
1190 14 : case 0: EXPR_WRAP(code, sumalt (EXPR_ARG,a,prec));
1191 7 : case 1: EXPR_WRAP(code, sumalt2(EXPR_ARG,a,prec));
1192 7 : default: pari_err_FLAG("sumalt");
1193 : }
1194 : return NULL; /* LCOV_EXCL_LINE */
1195 : }
1196 :
1197 : /* For k > 0, set S[k*2^i] <- g(k*2^i), k*2^i <= N = #S.
1198 : * Only needed with k odd (but also works for g even). */
1199 : static void
1200 8953 : binsum(GEN S, ulong k, void *E, GEN (*f)(void *, GEN), GEN a,
1201 : long G, long prec)
1202 : {
1203 8953 : long e, i, N = lg(S)-1, l = expu(N / k); /* k 2^l <= N < k 2^(l+1) */
1204 8953 : pari_sp av = avma;
1205 8953 : GEN t = real_0(prec); /* unused unless f(a + k <<l) = 0 */
1206 :
1207 8953 : G -= l;
1208 8953 : if (!signe(a)) a = NULL;
1209 8953 : for (e = 0;; e++)
1210 5389657 : { /* compute g(k 2^l) with absolute error ~ 2^(G-l) */
1211 5398610 : GEN u, r = shifti(utoipos(k), l+e);
1212 5398610 : if (a) r = addii(r, a);
1213 5398610 : u = gtofp(f(E, r), prec);
1214 5398610 : if (typ(u) != t_REAL) pari_err_TYPE("sumpos",u);
1215 5398610 : if (!signe(u)) break;
1216 5398421 : if (!e)
1217 8764 : t = u;
1218 : else {
1219 5389657 : shiftr_inplace(u, e);
1220 5389657 : t = addrr(t,u); if (expo(u) < G) break;
1221 5380893 : if ((e & 0x1ff) == 0) t = gerepileuptoleaf(av, t);
1222 : }
1223 : }
1224 8953 : gel(S, k << l) = t = gerepileuptoleaf(av, t);
1225 : /* g(j) = 2g(2j) + f(a+j) for all j > 0 */
1226 17906 : for(i = l-1; i >= 0; i--)
1227 : { /* t ~ g(2 * k*2^i) with error ~ 2^(G-i-1) */
1228 : GEN u;
1229 8953 : av = avma; u = gtofp(f(E, a? addiu(a, k << i): utoipos(k << i)), prec);
1230 8953 : if (typ(u) != t_REAL) pari_err_TYPE("sumpos",u);
1231 8953 : t = addrr(gtofp(u,prec), mpshift(t,1)); /* ~ g(k*2^i) */
1232 8953 : gel(S, k << i) = t = gerepileuptoleaf(av, t);
1233 : }
1234 8953 : }
1235 : /* For k > 0, let g(k) := \sum_{e >= 0} 2^e f(a + k*2^e).
1236 : * Return [g(k), 1 <= k <= N] */
1237 : static GEN
1238 84 : sumpos_init(void *E, GEN (*f)(void *, GEN), GEN a, long N, long prec)
1239 : {
1240 84 : GEN S = cgetg(N+1,t_VEC);
1241 84 : long k, G = -prec2nbits(prec) - 5;
1242 9037 : for (k=1; k<=N; k+=2) binsum(S,k, E,f, a,G,prec);
1243 84 : return S;
1244 : }
1245 :
1246 : GEN
1247 70 : sumpos(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
1248 : {
1249 : ulong k, N;
1250 70 : pari_sp av = avma;
1251 : GEN s, az, c, d, S;
1252 :
1253 70 : if (typ(a) != t_INT) pari_err_TYPE("sumpos",a);
1254 70 : a = subiu(a,1);
1255 70 : N = sumalt_N(prec);
1256 70 : if (odd(N)) N++; /* extra precision for free */
1257 70 : d = powru(addsr(3, sqrtr(utor(8,prec))), N);
1258 70 : d = shiftr(addrr(d, invr(d)),-1);
1259 70 : az = gen_m1; c = d;
1260 :
1261 70 : S = sumpos_init(E, eval, a, N, prec);
1262 70 : s = gen_0;
1263 13454 : for (k=0; k<N; k++)
1264 : {
1265 : GEN t;
1266 13454 : c = addir(az,c);
1267 13454 : t = mulrr(gel(S,k+1), c);
1268 13454 : s = odd(k)? mpsub(s, t): mpadd(s, t);
1269 13454 : if (k == N-1) break;
1270 13384 : az = diviuuexact(muluui((N-k)<<1,N+k,az), k+1, (k<<1)+1);
1271 : }
1272 70 : return gerepileupto(av, gdiv(s,d));
1273 : }
1274 :
1275 : GEN
1276 14 : sumpos2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
1277 : {
1278 : ulong k, N;
1279 14 : pari_sp av = avma;
1280 : GEN s, pol, dn, S;
1281 :
1282 14 : if (typ(a) != t_INT) pari_err_TYPE("sumpos2",a);
1283 14 : a = subiu(a,1);
1284 14 : N = (ulong)(0.31*(prec2nbits(prec) + 5));
1285 :
1286 14 : if (odd(N)) N++; /* extra precision for free */
1287 14 : S = sumpos_init(E, eval, a, N, prec);
1288 14 : pol = ZX_div_by_X_1(polzag1(N,N>>1), &dn);
1289 14 : s = gen_0;
1290 4466 : for (k=0; k<N; k++)
1291 : {
1292 4452 : GEN t = mulri(gel(S,k+1), gel(pol,k+2));
1293 4452 : s = odd(k)? mpsub(s,t): mpadd(s,t);
1294 : }
1295 14 : return gerepileupto(av, gdiv(s,dn));
1296 : }
1297 :
1298 : GEN
1299 91 : sumpos0(GEN a, GEN code, long flag, long prec)
1300 : {
1301 91 : switch(flag)
1302 : {
1303 70 : case 0: EXPR_WRAP(code, sumpos (EXPR_ARG,a,prec));
1304 14 : case 1: EXPR_WRAP(code, sumpos2(EXPR_ARG,a,prec));
1305 7 : default: pari_err_FLAG("sumpos");
1306 : }
1307 : return NULL; /* LCOV_EXCL_LINE */
1308 : }
1309 :
1310 : /********************************************************************/
1311 : /** **/
1312 : /** SEARCH FOR REAL ZEROS of an expression **/
1313 : /** **/
1314 : /********************************************************************/
1315 : /* Brent's method, [a,b] bracketing interval */
1316 : GEN
1317 23373 : zbrent(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
1318 : {
1319 : long sig, iter, itmax, bit, bit0;
1320 23373 : pari_sp av = avma;
1321 : GEN c, d, e, fa, fb, fc;
1322 :
1323 23373 : if (typ(a) == t_INFINITY && typ(b) != t_INFINITY) swap(a,b);
1324 23373 : if (typ(a) == t_INFINITY && typ(b) == t_INFINITY)
1325 : {
1326 7 : long s = gsigne(eval(E, real_0(prec))), r = 0;
1327 7 : if (gidentical(gel(a,1), gel(b,1)))
1328 0 : pari_err_DOMAIN("solve", "a and b", "=", a, mkvec2(a, b));
1329 7 : a = real_m1(prec); /* domain = R */
1330 7 : b = real_1(prec);
1331 : for(;;)
1332 : {
1333 7 : fa = eval(E, a);
1334 7 : fb = eval(E, b);
1335 7 : if (gsigne(fa) != s)
1336 : {
1337 0 : if (r) b[1] = evalsigne(-1) | _evalexpo(r-1); else b = real_0(prec);
1338 0 : break;
1339 : }
1340 7 : if (gsigne(fb) != s)
1341 : {
1342 7 : if (r) a[1] = evalsigne(1) | _evalexpo(r-1); else a = real_0(prec);
1343 7 : break;
1344 : }
1345 0 : r++; setexpo(a, r); setexpo(b, r);
1346 : }
1347 7 : c = b;
1348 7 : goto SOLVE;
1349 : }
1350 23366 : if (typ(b) == t_INFINITY)
1351 : { /* a real, b == [+-]oo */
1352 28 : long s, r, minf = inf_get_sign(b) < 0;
1353 : GEN inc;
1354 28 : if (typ(a) != t_REAL || realprec(a) < prec) a = gtofp(a, prec);
1355 28 : fa = eval(E, a);
1356 28 : s = gsigne(fa);
1357 28 : inc = minf ? real_m1(prec) : real_1(prec);
1358 28 : r = gsigne(a) ? expo(a) : 0;
1359 : for(;;)
1360 : {
1361 577 : setexpo(inc, r);
1362 577 : b = addrr(a, inc); fb = eval(E, b);
1363 563 : if (gsigne(fb) != s) break;
1364 549 : a = b; fa = fb; r++;
1365 : }
1366 14 : if (minf) { c = a; swap(a, b); swap(fa, fb);} else c = b;
1367 14 : goto SOLVE;
1368 : }
1369 23338 : if (typ(a) != t_REAL || realprec(a) < prec) a = gtofp(a, prec);
1370 23338 : if (typ(b) != t_REAL || realprec(b) < prec) b = gtofp(b, prec);
1371 23338 : sig = cmprr(b, a);
1372 23338 : if (!sig) return gerepileupto(av, a);
1373 23338 : if (sig < 0) { c = a; swap(a, b); } else c = b;
1374 23338 : fa = eval(E, a);
1375 23338 : fb = eval(E, b);
1376 23338 : if (gsigne(fa)*gsigne(fb) > 0)
1377 7 : pari_err_DOMAIN("solve", "f(a)f(b)", ">", gen_0, mkvec2(fa, fb));
1378 23331 : SOLVE:
1379 23352 : bit0 = -prec2nbits(prec); bit = 3+bit0; itmax = 1 - 2*bit0;
1380 23352 : fc = fb; e = d = NULL;
1381 225474 : for (iter = 1; iter <= itmax; ++iter)
1382 : { /* b = current best guess, a and c auxiliary points */
1383 : long bit2;
1384 : GEN m;
1385 225474 : if (gsigne(fb)*gsigne(fc) > 0) { c = a; fc = fa; e = d = subrr(b, a); }
1386 225474 : if (gcmp(gabs(fc, 0), gabs(fb, 0)) < 0)
1387 50975 : { a = b; b = c; c = a; fa = fb; fb = fc; fc = fa; }
1388 225474 : bit2 = bit-1 + maxss(bit, expo(b));
1389 225474 : m = subrr(c, b); shiftr_inplace(m, -1);
1390 225474 : if (expo(m) <= expo(b)+bit0 || gequal0(fb)) break; /*SUCCESS*/
1391 :
1392 202122 : if (expo(e) > bit2 && gexpo(fa) > gexpo(fb))
1393 159442 : { /* interpolation, m != 0, |f(c)| >= |f(b)|, f(b)f(c) < 0 */
1394 159442 : GEN min1, min2, p, q, s = gdiv(fb, fa);
1395 159442 : if (a == c || equalrr(a,c))
1396 : {
1397 134012 : p = gmul2n(gmul(m, s), 1);
1398 134012 : q = gsubsg(1, s);
1399 : }
1400 : else
1401 : {
1402 25430 : GEN r = gdiv(fb, fc);
1403 25430 : q = gdiv(fa, fc);
1404 25430 : p = gmul2n(gmul(gsub(q, r), gmul(m, q)), 1);
1405 25430 : p = gmul(s, gsub(p, gmul(subrr(b, a), gsubgs(r, 1))));
1406 25430 : q = gmul(gmul(gsubgs(q, 1), gsubgs(r, 1)), gsubgs(s, 1));
1407 : }
1408 159442 : if (gsigne(p) > 0) q = gneg_i(q); else p = gneg_i(p);
1409 159442 : min1 = gsub(gmulsg(3, gmul(m,q)), gmul2n(gabs(q,0), bit2));
1410 159442 : min2 = gabs(gmul(e, q), 0);
1411 159442 : if (gcmp(gmul2n(p, 1), gmin_shallow(min1, min2)) < 0)
1412 158421 : { e = d; d = gdiv(p, q); } /* interpolation OK */
1413 : else
1414 1021 : e = d = m; /* failed, use bisection */
1415 : }
1416 42680 : else e = d = m; /* bound decreasing too slowly, use bisection */
1417 202122 : a = b; fa = fb;
1418 202122 : if (d == m) { b = addrr(c, b); shiftr_inplace(b,-1); }
1419 158421 : else if (gexpo(d) > bit2) b = gadd(b, d);
1420 23985 : else if (gsigne(m) > 0) b = addrr(b, real2n(bit2, LOWDEFAULTPREC));
1421 11862 : else b = subrr(b, real2n(bit2, LOWDEFAULTPREC));
1422 202122 : if (equalrr(a, b)) break;
1423 202122 : if (realprec(b) < prec) b = rtor(b, prec);
1424 202122 : fb = eval(E, b);
1425 : }
1426 23352 : if (iter > itmax) pari_err_IMPL("solve recovery [too many iterations]");
1427 23352 : return gerepileuptoleaf(av, rcopy(b));
1428 : }
1429 :
1430 : GEN
1431 63 : zbrent0(GEN a, GEN b, GEN code, long prec)
1432 63 : { EXPR_WRAP(code, zbrent(EXPR_ARG, a, b, prec)); }
1433 :
1434 : /* Find zeros of a function in the real interval [a,b] by interval splitting */
1435 : GEN
1436 119 : solvestep(void *E, GEN (*f)(void *,GEN), GEN a, GEN b, GEN step, long flag, long prec)
1437 : {
1438 119 : const long ITMAX = 10;
1439 119 : pari_sp av = avma;
1440 : GEN fa, a0, b0;
1441 119 : long sa0, it, bit = prec2nbits(prec) / 2, ct = 0, s = gcmp(a,b);
1442 :
1443 119 : if (!s) return gequal0(f(E, a)) ? gcopy(mkvec(a)): cgetg(1,t_VEC);
1444 119 : if (s > 0) swap(a, b);
1445 119 : if (flag&4)
1446 : {
1447 84 : if (gcmpgs(step,1)<=0) pari_err_DOMAIN("solvestep","step","<=",gen_1,step);
1448 84 : if (gsigne(a) <= 0) pari_err_DOMAIN("solvestep","a","<=",gen_0,a);
1449 : }
1450 35 : else if (gsigne(step) <= 0)
1451 7 : pari_err_DOMAIN("solvestep","step","<=",gen_0,step);
1452 112 : a0 = a = gtofp(a, prec); fa = f(E, a);
1453 112 : b0 = b = gtofp(b, prec); step = gtofp(step, prec);
1454 112 : sa0 = gsigne(fa);
1455 112 : if (gexpo(fa) < -bit) sa0 = 0;
1456 119 : for (it = 0; it < ITMAX; it++)
1457 : {
1458 119 : pari_sp av2 = avma;
1459 119 : GEN v = cgetg(1, t_VEC);
1460 119 : long sa = sa0;
1461 119 : a = a0; b = b0;
1462 37520 : while (gcmp(a,b) < 0)
1463 : {
1464 37401 : GEN fc, c = (flag&4)? gmul(a, step): gadd(a, step);
1465 : long sc;
1466 37401 : if (gcmp(c,b) > 0) c = b;
1467 37401 : fc = f(E, c); sc = gsigne(fc);
1468 37401 : if (gexpo(fc) < -bit) sc = 0;
1469 37401 : if (!sc || sa*sc < 0)
1470 : {
1471 22813 : GEN z = sc? zbrent(E, f, a, c, prec): c;
1472 : long e;
1473 22813 : (void)grndtoi(z, &e);
1474 22813 : if (e <= -bit) ct = 1;
1475 22813 : if ((flag&1) && ((!(flag&8)) || ct)) return gerepileupto(av, z);
1476 22813 : v = shallowconcat(v, z);
1477 : }
1478 37401 : a = c; fa = fc; sa = sc;
1479 37401 : if (gc_needed(av2,1))
1480 : {
1481 52 : if (DEBUGMEM>1) pari_warn(warnmem,"solvestep");
1482 52 : gerepileall(av2, 4, &a ,&fa, &v, &step);
1483 : }
1484 : }
1485 119 : if ((!(flag&2) || lg(v) > 1) && (!(flag&8) || ct))
1486 112 : return gerepilecopy(av, v);
1487 7 : step = (flag&4)? sqrtnr(step,4): gmul2n(step, -2);
1488 7 : gerepileall(av2, 2, &fa, &step);
1489 : }
1490 0 : pari_err_IMPL("solvestep recovery [too many iterations]");
1491 : return NULL;/*LCOV_EXCL_LINE*/
1492 : }
1493 :
1494 : GEN
1495 35 : solvestep0(GEN a, GEN b, GEN step, GEN code, long flag, long prec)
1496 35 : { EXPR_WRAP(code, solvestep(EXPR_ARG, a,b, step, flag, prec)); }
1497 :
1498 : /********************************************************************/
1499 : /** Numerical derivation **/
1500 : /********************************************************************/
1501 :
1502 : struct deriv_data
1503 : {
1504 : GEN code;
1505 : GEN args;
1506 : GEN def;
1507 : };
1508 :
1509 322 : static GEN deriv_eval(void *E, GEN x, long prec)
1510 : {
1511 322 : struct deriv_data *data=(struct deriv_data *)E;
1512 322 : gel(data->args,1)=x;
1513 322 : uel(data->def,1)=1;
1514 322 : return closure_callgenvecdefprec(data->code, data->args, data->def, prec);
1515 : }
1516 :
1517 : /* Rationale: (f(2^-e) - f(-2^-e) + O(2^-b)) / (2 * 2^-e) = f'(0) + O(2^-2e)
1518 : * since 2nd derivatives cancel.
1519 : * prec(LHS) = b - e
1520 : * prec(RHS) = 2e, equal when b = 3e = 3/2 b0 (b0 = required final bitprec)
1521 : *
1522 : * For f'(x), x far from 0: prec(LHS) = b - e - expo(x)
1523 : * --> pr = 3/2 b0 + expo(x) */
1524 : GEN
1525 966 : derivnum(void *E, GEN (*eval)(void *, GEN, long), GEN x, long prec)
1526 : {
1527 966 : long newprec, e, ex = gexpo(x), p = precision(x);
1528 966 : long b0 = prec2nbits(p? p: prec), b = (long)ceil(b0 * 1.5 + maxss(0,ex));
1529 : GEN eps, u, v, y;
1530 966 : pari_sp av = avma;
1531 966 : newprec = nbits2prec(b + BITS_IN_LONG);
1532 966 : switch(typ(x))
1533 : {
1534 385 : case t_REAL:
1535 : case t_COMPLEX:
1536 385 : x = gprec_w(x, newprec);
1537 : }
1538 966 : e = b0/2; /* 1/2 required prec (in sig. bits) */
1539 966 : b -= e; /* >= b0 */
1540 966 : eps = real2n(-e, ex < -e? newprec: nbits2prec(b));
1541 966 : u = eval(E, gsub(x, eps), newprec);
1542 966 : v = eval(E, gadd(x, eps), newprec);
1543 966 : y = gmul2n(gsub(v,u), e-1);
1544 966 : return gerepilecopy(av, gprec_wtrunc(y, nbits2prec(b0)));
1545 : }
1546 :
1547 : /* Fornberg interpolation algorithm for finite differences coefficients
1548 : * using 2N+1 equidistant grid points around 0 [ assume 2N even >= M ].
1549 : * Compute \delta[m]_{N,i} for all derivation orders m = 0..M such that
1550 : * h^m * f^{(m)}(0) = \sum_{i = 0}^n delta[m]_{N,i} f(a_i) + O(h^{N-m+1}),
1551 : * for step size h.
1552 : * Return a = [0,-1,1...,-N,N] and vector of vectors d: d[m+1][i+1]
1553 : * = w'(a_i) delta[m]_{2N,i}, i = 0..2N */
1554 : static void
1555 112 : FD(long M, long N2, GEN *pd, GEN *pa)
1556 : {
1557 : GEN d, a, b, W, F;
1558 112 : long N = N2>>1, m, i;
1559 :
1560 112 : F = cgetg(N2+2, t_VEC);
1561 112 : a = cgetg(N2+2, t_VEC);
1562 112 : b = cgetg(N+1, t_VEC);
1563 112 : gel(a,1) = gen_0;
1564 623 : for (i = 1; i <= N; i++)
1565 : {
1566 511 : gel(a,2*i) = utoineg(i);
1567 511 : gel(a,2*i+1) = utoipos(i);
1568 511 : gel(b,i) = sqru(i);
1569 : }
1570 : /* w = \prod (X - a[i]) = x W(x^2) */
1571 112 : W = roots_to_pol(b, 0);
1572 112 : gel(F,1) = RgX_inflate(W,2);
1573 623 : for (i = 1; i <= N; i++)
1574 : {
1575 511 : pari_sp av = avma;
1576 : GEN r, U, S;
1577 511 : U = RgX_inflate(RgX_div_by_X_x(W, gel(b,i), &r), 2);
1578 511 : U = RgXn_red_shallow(U, M); /* higher terms not needed */
1579 511 : U = RgX_shift_shallow(U,1); /* w(X) / (X^2-a[i]^2) mod X^(M+1) */
1580 511 : S = ZX_sub(RgX_shift_shallow(U,1),
1581 511 : ZX_Z_mul(U, gel(a,2*i+1)));
1582 511 : S = gerepileupto(av, S);
1583 511 : gel(F,2*i) = S;
1584 511 : gel(F,2*i+1) = ZX_z_unscale(S, -1);
1585 : }
1586 : /* F[i] = w(X) / (X-a[i]) + O(X^(M+1)) in Z[X] */
1587 112 : d = cgetg(M+2, t_VEC);
1588 581 : for (m = 0; m <= M; m++)
1589 : {
1590 469 : GEN v = cgetg(N2+2, t_VEC); /* coeff(F[i],X^m) */
1591 11550 : for (i = 0; i <= N2; i++) gel(v, i+1) = gmael(F, i+1, m+2);
1592 469 : gel(d,m+1) = v;
1593 : }
1594 112 : *pd = d;
1595 112 : *pa = a;
1596 112 : }
1597 :
1598 : static void
1599 329 : chk_ord(long m)
1600 : {
1601 329 : if (m < 0)
1602 14 : pari_err_DOMAIN("derivnumk", "derivation order", "<", gen_0, stoi(m));
1603 315 : }
1604 : /* m! / N! for m in ind; vecmax(ind) <= N */
1605 : static GEN
1606 112 : vfact(GEN ind, long N, long prec)
1607 : {
1608 : GEN v, iN;
1609 : long i, l;
1610 112 : ind = vecsmall_uniq(ind); chk_ord(ind[1]); l = lg(ind);
1611 105 : iN = invr(itor(mulu_interval(ind[1] + 1, N), prec));
1612 105 : v = const_vec(ind[l-1], NULL); gel(v, ind[1]) = iN;
1613 196 : for (i = 2; i < l; i++)
1614 91 : gel(v, ind[i]) = iN = mulri(iN, mulu_interval(ind[i-1] + 1, ind[i]));
1615 105 : return v;
1616 : }
1617 :
1618 : static GEN
1619 175 : chk_ind(GEN ind, long *M)
1620 : {
1621 175 : *M = 0;
1622 175 : switch(typ(ind))
1623 : {
1624 63 : case t_INT: ind = mkvecsmall(itos(ind)); break;
1625 0 : case t_VECSMALL:
1626 0 : if (lg(ind) == 1) return NULL;
1627 0 : break;
1628 105 : case t_VEC: case t_COL:
1629 105 : if (lg(ind) == 1) return NULL;
1630 98 : if (RgV_is_ZV(ind)) { ind = ZV_to_zv(ind); break; }
1631 : /* fall through */
1632 : default:
1633 7 : pari_err_TYPE("derivnum", ind);
1634 : return NULL; /*LCOV_EXCL_LINE*/
1635 : }
1636 161 : *M = vecsmall_max(ind); chk_ord(*M); return ind;
1637 : }
1638 : GEN
1639 140 : derivnumk(void *E, GEN (*eval)(void *, GEN, long), GEN x, GEN ind0, long prec)
1640 : {
1641 : GEN A, C, D, DM, T, X, F, v, ind, t;
1642 : long M, N, N2, fpr, p, i, pr, l, lA, e, ex, emin, emax, newprec;
1643 140 : pari_sp av = avma;
1644 140 : int allodd = 1;
1645 :
1646 140 : ind = chk_ind(ind0, &M); if (!ind) return cgetg(1, t_VEC);
1647 126 : l = lg(ind); F = cgetg(l, t_VEC);
1648 126 : if (!M) /* silly degenerate case */
1649 : {
1650 14 : X = eval(E, x, prec);
1651 28 : for (i = 1; i < l; i++) { chk_ord(ind[i]); gel(F,i) = X; }
1652 7 : if (typ(ind0) == t_INT) F = gel(F,1);
1653 7 : return gerepilecopy(av, F);
1654 : }
1655 112 : N2 = 3*M - 1; if (odd(N2)) N2++;
1656 112 : N = N2 >> 1;
1657 112 : FD(M, N2, &D,&A); /* optimal if 'eval' uses quadratic time */
1658 112 : C = vecbinomial(N2); DM = gel(D,M);
1659 112 : T = cgetg(N2+2, t_VEC);
1660 : /* (2N)! / w'(i) = (2N)! / w'(-i) = (-1)^(N-i) binom(2*N, N-i) */
1661 112 : t = gel(C, N+1);
1662 112 : gel(T,1) = odd(N)? negi(t): t;
1663 623 : for (i = 1; i <= N; i++)
1664 : {
1665 511 : t = gel(C, N-i+1);
1666 511 : gel(T,2*i) = gel(T,2*i+1) = odd(N-i)? negi(t): t;
1667 : }
1668 112 : N = N2 >> 1; emin = LONG_MAX; emax = 0;
1669 623 : for (i = 1; i <= N; i++)
1670 : {
1671 511 : e = expi(gel(DM,i)) + expi(gel(T,i));
1672 511 : if (e < 0) continue; /* 0 */
1673 448 : if (e < emin) emin = e;
1674 252 : else if (e > emax) emax = e;
1675 : }
1676 :
1677 112 : p = precision(x);
1678 112 : fpr = p ? prec2nbits(p): prec2nbits(prec);
1679 112 : e = (fpr + 3*M*log2((double)M)) / (2*M);
1680 112 : ex = gexpo(x);
1681 112 : if (ex < 0) ex = 0; /* near 0 */
1682 112 : pr = (long)ceil(fpr + e * M); /* ~ 3fpr/2 */
1683 112 : newprec = nbits2prec(pr + (emax - emin) + ex + BITS_IN_LONG);
1684 112 : switch(typ(x))
1685 : {
1686 28 : case t_REAL:
1687 : case t_COMPLEX:
1688 28 : x = gprec_w(x, newprec);
1689 : }
1690 112 : lA = lg(A); X = cgetg(lA, t_VEC);
1691 161 : for (i = 1; i < l; i++)
1692 119 : if (!odd(ind[i])) { allodd = 0; break; }
1693 : /* if only odd derivation orders, the value at 0 (A[1]) is not needed */
1694 112 : gel(X, 1) = gen_0;
1695 1204 : for (i = allodd? 2: 1; i < lA; i++)
1696 : {
1697 1092 : GEN t = eval(E, gadd(x, gmul2n(gel(A,i), -e)), newprec);
1698 1092 : t = gmul(t, gel(T,i));
1699 1092 : if (!gprecision(t))
1700 224 : t = is_scalar_t(typ(t))? gtofp(t, newprec): gmul(t, real_1(newprec));
1701 1092 : gel(X,i) = t;
1702 : }
1703 :
1704 112 : v = vfact(ind, N2, nbits2prec(fpr + 32));
1705 301 : for (i = 1; i < l; i++)
1706 : {
1707 196 : long m = ind[i];
1708 196 : GEN t = RgV_dotproduct(gel(D,m+1), X);
1709 196 : gel(F,i) = gmul(t, gmul2n(gel(v, m), e*m));
1710 : }
1711 105 : if (typ(ind0) == t_INT) F = gel(F,1);
1712 105 : return gerepilecopy(av, F);
1713 : }
1714 : /* v(t') */
1715 : static long
1716 14 : rfrac_val_deriv(GEN t)
1717 : {
1718 14 : long v = varn(gel(t,2));
1719 14 : return gvaluation(deriv(t, v), pol_x(v));
1720 : }
1721 :
1722 : GEN
1723 1190 : derivfunk(void *E, GEN (*eval)(void *, GEN, long), GEN x, GEN ind0, long prec)
1724 : {
1725 : pari_sp av;
1726 : GEN ind, xp, ixp, F, G;
1727 : long i, l, vx, M;
1728 1190 : if (!ind0) return derivfun(E, eval, x, prec);
1729 203 : switch(typ(x))
1730 : {
1731 140 : case t_REAL: case t_INT: case t_FRAC: case t_COMPLEX:
1732 140 : return derivnumk(E,eval, x, ind0, prec);
1733 21 : case t_POL:
1734 21 : ind = chk_ind(ind0,&M); if (!ind) return cgetg(1,t_VEC);
1735 21 : xp = RgX_deriv(x);
1736 21 : x = RgX_to_ser(x, precdl+2 + M * (1+RgX_val(xp)));
1737 21 : break;
1738 7 : case t_RFRAC:
1739 7 : ind = chk_ind(ind0,&M); if (!ind) return cgetg(1,t_VEC);
1740 7 : x = rfrac_to_ser_i(x, precdl+2 + M * (1+rfrac_val_deriv(x)));
1741 7 : xp = derivser(x);
1742 7 : break;
1743 7 : case t_SER:
1744 7 : ind = chk_ind(ind0,&M); if (!ind) return cgetg(1,t_VEC);
1745 7 : xp = derivser(x);
1746 7 : break;
1747 28 : default: pari_err_TYPE("numerical derivation",x);
1748 : return NULL; /*LCOV_EXCL_LINE*/
1749 : }
1750 35 : av = avma; vx = varn(x);
1751 35 : ixp = M? ginv(xp): NULL;
1752 35 : F = cgetg(M+2, t_VEC);
1753 35 : gel(F,1) = eval(E, x, prec);
1754 126 : for (i = 1; i <= M; i++) gel(F,i+1) = gmul(deriv(gel(F,i),vx), ixp);
1755 35 : l = lg(ind); G = cgetg(l, t_VEC);
1756 70 : for (i = 1; i < l; i++)
1757 : {
1758 35 : long m = ind[i]; chk_ord(m);
1759 35 : gel(G,i) = gel(F,m+1);
1760 : }
1761 35 : if (typ(ind0) == t_INT) G = gel(G,1);
1762 35 : return gerepilecopy(av, G);
1763 : }
1764 :
1765 : GEN
1766 987 : derivfun(void *E, GEN (*eval)(void *, GEN, long), GEN x, long prec)
1767 : {
1768 987 : pari_sp av = avma;
1769 : GEN xp;
1770 : long vx;
1771 987 : switch(typ(x))
1772 : {
1773 966 : case t_REAL: case t_INT: case t_FRAC: case t_COMPLEX:
1774 966 : return derivnum(E,eval, x, prec);
1775 7 : case t_POL:
1776 7 : xp = RgX_deriv(x);
1777 7 : x = RgX_to_ser(x, precdl+2+ (1 + RgX_val(xp)));
1778 7 : break;
1779 7 : case t_RFRAC:
1780 7 : x = rfrac_to_ser_i(x, precdl+2+ (1 + rfrac_val_deriv(x)));
1781 : /* fall through */
1782 14 : case t_SER:
1783 14 : xp = derivser(x);
1784 14 : break;
1785 0 : default: pari_err_TYPE("formal derivation",x);
1786 : return NULL; /*LCOV_EXCL_LINE*/
1787 : }
1788 21 : vx = varn(x);
1789 21 : return gerepileupto(av, gdiv(deriv(eval(E, x, prec),vx), xp));
1790 : }
1791 :
1792 : GEN
1793 21 : laurentseries(void *E, GEN (*f)(void*,GEN x, long), long M, long v, long prec)
1794 : {
1795 21 : pari_sp av = avma;
1796 : long d;
1797 :
1798 21 : if (v < 0) v = 0;
1799 21 : d = maxss(M+1,1);
1800 : for (;;)
1801 14 : {
1802 : long i, dr, vr;
1803 : GEN s;
1804 35 : s = cgetg(d+2, t_SER); s[1] = evalsigne(1) | evalvalp(1) | evalvarn(v);
1805 245 : gel(s, 2) = gen_1; for (i = 3; i <= d+1; i++) gel(s, i) = gen_0;
1806 35 : s = f(E, s, prec);
1807 35 : if (typ(s) != t_SER || varn(s) != v) pari_err_TYPE("laurentseries", s);
1808 35 : vr = valp(s);
1809 35 : if (M < vr) { set_avma(av); return zeroser(v, M); }
1810 35 : dr = lg(s) + vr - 3 - M;
1811 35 : if (dr >= 0) return gerepileupto(av, s);
1812 14 : set_avma(av); d -= dr;
1813 : }
1814 : }
1815 : static GEN
1816 35 : _evalclosprec(void *E, GEN x, long prec)
1817 : {
1818 : GEN s;
1819 35 : push_localprec(prec); s = closure_callgen1((GEN)E, x);
1820 35 : pop_localprec(); return s;
1821 : }
1822 : #define CLOS_ARGPREC __E, &_evalclosprec
1823 : GEN
1824 35 : laurentseries0(GEN f, long M, long v, long prec)
1825 : {
1826 35 : if (typ(f) != t_CLOSURE || closure_arity(f) != 1 || closure_is_variadic(f))
1827 14 : pari_err_TYPE("laurentseries",f);
1828 21 : EXPR_WRAP(f, laurentseries(CLOS_ARGPREC,M,v,prec));
1829 : }
1830 :
1831 : GEN
1832 1085 : derivnum0(GEN a, GEN code, GEN ind, long prec)
1833 1085 : { EXPR_WRAP(code, derivfunk(EXPR_ARGPREC,a,ind,prec)); }
1834 :
1835 : GEN
1836 105 : derivfun0(GEN args, GEN def, GEN code, long k, long prec)
1837 : {
1838 105 : pari_sp av = avma;
1839 : struct deriv_data E;
1840 : GEN z;
1841 105 : E.code=code; E.args=args; E.def=def;
1842 105 : z = gel(derivfunk((void*)&E, deriv_eval, gel(args,1), mkvecs(k), prec),1);
1843 77 : return gerepilecopy(av, z);
1844 : }
1845 :
1846 : /********************************************************************/
1847 : /** Numerical extrapolation **/
1848 : /********************************************************************/
1849 : /* [u(n), u <= N] */
1850 : static GEN
1851 140 : get_u(void *E, GEN (*f)(void *, GEN, long), long N, long prec)
1852 : {
1853 : long n;
1854 : GEN u;
1855 140 : if (f)
1856 : {
1857 126 : GEN v = f(E, utoipos(N), prec);
1858 126 : u = cgetg(N+1, t_VEC);
1859 126 : if (typ(v) != t_VEC || lg(v) != N+1) { gel(u,N) = v; v = NULL; }
1860 : else
1861 : {
1862 14 : GEN w = f(E, gen_1, LOWDEFAULTPREC);
1863 14 : if (typ(w) != t_VEC || lg(w) != 2) { gel(u,N) = v; v = NULL; }
1864 : }
1865 126 : if (v) u = v;
1866 : else
1867 9702 : for (n = 1; n < N; n++) gel(u,n) = f(E, utoipos(n), prec);
1868 : }
1869 : else
1870 : {
1871 14 : GEN v = (GEN)E;
1872 14 : long t = lg(v)-1;
1873 14 : if (t < N) pari_err_COMPONENT("limitnum","<",stoi(N), stoi(t));
1874 14 : u = vecslice(v, 1, N);
1875 : }
1876 12236 : for (n = 1; n <= N; n++)
1877 : {
1878 12096 : GEN un = gel(u,n);
1879 12096 : if (is_rational_t(typ(un))) gel(u,n) = gtofp(un, prec);
1880 : }
1881 140 : return u;
1882 : }
1883 :
1884 : struct limit
1885 : {
1886 : long prec; /* working accuracy */
1887 : long N; /* number of terms */
1888 : GEN na; /* [n^alpha, n <= N] */
1889 : GEN coef; /* or NULL (alpha != 1) */
1890 : };
1891 :
1892 : static GEN
1893 20822 : _gi(void *E, GEN x)
1894 : {
1895 20822 : GEN A = (GEN)E, y = gsubgs(x, 1);
1896 20822 : if (gequal0(y)) return A;
1897 20808 : return gdiv(gsubgs(gpow(x, A, LOWDEFAULTPREC), 1), y);
1898 : }
1899 : static GEN
1900 166 : _g(void *E, GEN x)
1901 : {
1902 166 : GEN D = (GEN)E, A = gel(D,1), T = gel(D,2);
1903 166 : const long prec = LOWDEFAULTPREC;
1904 166 : return gadd(glog(x,prec), intnum((void*)A, _gi, gen_0, gaddgs(x,1), T, prec));
1905 : }
1906 :
1907 : /* solve log(b) + int_0^{b+1} (x^(1/a)-1) / (x-1) dx = 0, b in [0,1]
1908 : * return -log_2(b), rounded up */
1909 : static double
1910 140 : get_accu(GEN a)
1911 : {
1912 140 : pari_sp av = avma;
1913 140 : const long prec = LOWDEFAULTPREC;
1914 140 : const double We2 = 1.844434455794; /* (W(1/e) + 1) / log(2) */
1915 : GEN b, T;
1916 140 : if (!a) return We2;
1917 49 : if (typ(a) == t_INT) switch(itos_or_0(a))
1918 : {
1919 0 : case 1: return We2;
1920 21 : case 2: return 1.186955309668;
1921 0 : case 3: return 0.883182331990;
1922 : }
1923 28 : else if (typ(a) == t_FRAC && equali1(gel(a,1))) switch(itos_or_0(gel(a,2)))
1924 : {
1925 14 : case 2: return 2.644090500290;
1926 0 : case 3: return 3.157759214459;
1927 0 : case 4: return 3.536383237500;
1928 : }
1929 14 : T = intnuminit(gen_0, gen_1, 0, prec);
1930 14 : b = zbrent((void*)mkvec2(ginv(a), T), &_g, dbltor(1E-5), gen_1, prec);
1931 14 : return gc_double(av, -dbllog2r(b));
1932 : }
1933 :
1934 : static double
1935 147 : get_c(GEN a)
1936 : {
1937 147 : double A = a? gtodouble(a): 1.0;
1938 147 : if (A <= 0) pari_err_DOMAIN("limitnum","alpha","<=",gen_0, a);
1939 140 : if (A >= 2) return 0.2270;
1940 105 : if (A >= 1) return 0.3318;
1941 14 : if (A >= 0.5) return 0.6212;
1942 0 : if (A >= 0.3333) return 1.2;
1943 0 : return 3; /* only tested for A >= 0.25 */
1944 : }
1945 : static void
1946 133 : limit_Nprec(struct limit *L, GEN alpha, long prec)
1947 : {
1948 133 : long bit = prec2nbits(prec);
1949 133 : L->N = ceil(get_c(alpha) * bit);
1950 126 : L->prec = nbits2prec(bit + (long)ceil(get_accu(alpha) * L->N));
1951 126 : }
1952 : /* solve x - a log(x) = b; a, b >= 3 */
1953 : static double
1954 14 : solvedivlog(double a, double b) { return dbllemma526(a,1,1,b); }
1955 :
1956 : /* #u > 1, prod_{j != i} u[i] - u[j] */
1957 : static GEN
1958 3003 : proddiff(GEN u, long i)
1959 : {
1960 3003 : pari_sp av = avma;
1961 3003 : long l = lg(u), j;
1962 3003 : GEN p = NULL;
1963 3003 : if (i == 1)
1964 : {
1965 28 : p = gsub(gel(u,1), gel(u,2));
1966 2975 : for (j = 3; j < l; j++)
1967 2947 : p = gmul(p, gsub(gel(u,i), gel(u,j)));
1968 : }
1969 : else
1970 : {
1971 2975 : p = gsub(gel(u,i), gel(u,1));
1972 367346 : for (j = 2; j < l; j++)
1973 364371 : if (j != i) p = gmul(p, gsub(gel(u,i), gel(u,j)));
1974 : }
1975 3003 : return gerepileupto(av, p);
1976 : }
1977 : static GEN
1978 14 : vecpows(GEN u, long N)
1979 : {
1980 : long i, l;
1981 14 : GEN v = cgetg_copy(u, &l);
1982 1883 : for (i = 1; i < l; i++) gel(v,i) = gpowgs(gel(u,i), N);
1983 14 : return v;
1984 : }
1985 :
1986 : static void
1987 140 : limit_init(struct limit *L, GEN alpha, int asymp)
1988 : {
1989 140 : long n, N = L->N, a = 0;
1990 140 : GEN c, v, T = NULL;
1991 :
1992 140 : if (!alpha) a = 1;
1993 49 : else if (typ(alpha) == t_INT)
1994 : {
1995 21 : a = itos_or_0(alpha);
1996 21 : if (a > 2) a = 0;
1997 : }
1998 28 : else if (typ(alpha) == t_FRAC)
1999 : {
2000 14 : long na = itos_or_0(gel(alpha,1)), da = itos_or_0(gel(alpha,2));
2001 14 : if (da && na && da <= 4 && na <= 4)
2002 : { /* don't bother with other cases */
2003 14 : long e = (N-1) % da, k = (N-1) / da;
2004 14 : if (e) { N += da - e; k++; } /* N = 1 (mod d) => simpler ^ (n/d)(N-1) */
2005 14 : L->N = N;
2006 14 : T = vecpowuu(N, na * k);
2007 : }
2008 : }
2009 140 : L->coef = v = cgetg(N+1, t_VEC);
2010 140 : if (!a)
2011 : {
2012 28 : long prec2 = gprecision(alpha);
2013 : GEN u;
2014 28 : if (prec2 && prec2 < L->prec) alpha = gprec_w(alpha, L->prec);
2015 28 : L->na = u = vecpowug(N, alpha, L->prec);
2016 28 : if (!T) T = vecpows(u, N-1);
2017 3031 : for (n = 1; n <= N; n++) gel(v,n) = gdiv(gel(T,n), proddiff(u,n));
2018 28 : return;
2019 : }
2020 112 : L->na = asymp? vecpowuu(N, a): NULL;
2021 112 : c = mpfactr(N-1, L->prec);
2022 112 : if (a == 1)
2023 : {
2024 91 : c = invr(c);
2025 91 : gel(v,1) = c; if (!odd(N)) togglesign(c);
2026 7651 : for (n = 2; n <= N; n++) gel(v,n) = divru(mulrs(gel(v,n-1), n-1-N), n);
2027 : }
2028 : else
2029 : { /* a = 2 */
2030 21 : c = invr(mulru(sqrr(c), (N*(N+1)) >> 1));
2031 21 : gel(v,1) = c; if (!odd(N)) togglesign(c);
2032 1442 : for (n = 2; n <= N; n++) gel(v,n) = divru(mulrs(gel(v,n-1), n-1-N), N+n);
2033 : }
2034 112 : T = vecpowuu(N, a*N);
2035 9093 : for (n = 2; n <= N; n++) gel(v,n) = mulri(gel(v,n), gel(T,n));
2036 : }
2037 :
2038 : /* Zagier/Lagrange extrapolation */
2039 : static GEN
2040 983 : limitnum_i(struct limit *L, GEN u, long prec)
2041 983 : { return gprec_w(RgV_dotproduct(u,L->coef), prec); }
2042 : GEN
2043 84 : limitnum(void *E, GEN (*f)(void *, GEN, long), GEN alpha, long prec)
2044 : {
2045 84 : pari_sp av = avma;
2046 : struct limit L;
2047 : GEN u;
2048 84 : limit_Nprec(&L, alpha, prec);
2049 77 : limit_init(&L, alpha, 0);
2050 77 : u = get_u(E, f, L.N, L.prec);
2051 77 : return gerepilecopy(av, limitnum_i(&L, u, prec));
2052 : }
2053 : typedef GEN (*LIMIT_FUN)(void*,GEN,long);
2054 161 : static LIMIT_FUN get_fun(GEN u, const char *s)
2055 : {
2056 161 : switch(typ(u))
2057 : {
2058 14 : case t_COL: case t_VEC: break;
2059 133 : case t_CLOSURE: return gp_callprec;
2060 14 : default: pari_err_TYPE(s, u);
2061 : }
2062 14 : return NULL;
2063 : }
2064 : GEN
2065 91 : limitnum0(GEN u, GEN alpha, long prec)
2066 91 : { return limitnum((void*)u,get_fun(u, "limitnum"), alpha, prec); }
2067 :
2068 : GEN
2069 49 : asympnum(void *E, GEN (*f)(void *, GEN, long), GEN alpha, long prec)
2070 : {
2071 49 : const long MAX = 100;
2072 49 : pari_sp av = avma;
2073 49 : GEN u, A = cgetg(MAX+1, t_VEC);
2074 49 : long i, B = prec2nbits(prec);
2075 49 : double LB = 0.9*expu(B); /* 0.9 and 0.95 below are heuristic */
2076 : struct limit L;
2077 49 : limit_Nprec(&L, alpha, prec);
2078 49 : if (alpha) LB *= gtodouble(alpha);
2079 49 : limit_init(&L, alpha, 1);
2080 49 : u = get_u(E, f, L.N, L.prec);
2081 906 : for(i = 1; i <= MAX; i++)
2082 : {
2083 906 : GEN a, v, q, s = limitnum_i(&L, u, prec);
2084 : long n;
2085 : /* NOT bestappr: lindep properly ignores the lower bits */
2086 906 : v = lindep_bit(mkvec2(gen_1, s), maxss((long)(0.95*floor(B - i*LB)), 32));
2087 906 : if (lg(v) == 1) break;
2088 899 : q = gel(v,2); if (!signe(q)) break;
2089 899 : a = gdiv(negi(gel(v,1)), q);
2090 899 : s = gsub(s, a);
2091 : /* |s|q^2 > eps */
2092 899 : if (!gequal0(s) && gexpo(s) + 2*expi(q) > -17) break;
2093 857 : gel(A,i) = a;
2094 82293 : for (n = 1; n <= L.N; n++) gel(u,n) = gmul(gsub(gel(u,n), a), gel(L.na,n));
2095 : }
2096 49 : setlg(A,i); return gerepilecopy(av, A);
2097 : }
2098 : GEN
2099 14 : asympnumraw(void *E, GEN (*f)(void *,GEN,long), long LIM, GEN alpha, long prec)
2100 : {
2101 14 : pari_sp av = avma;
2102 : double c, d, al;
2103 : long i, B;
2104 : GEN u, A;
2105 : struct limit L;
2106 :
2107 14 : if (LIM < 0) return cgetg(1, t_VEC);
2108 14 : c = get_c(alpha);
2109 14 : d = get_accu(alpha);
2110 14 : al = alpha? gtodouble(alpha): 1.0;
2111 14 : B = prec2nbits(prec);
2112 14 : L.N = ceil(solvedivlog(c * al * LIM / M_LN2, c * B));
2113 14 : L.prec = nbits2prec(ceil(B + L.N / c + d * L.N));
2114 14 : limit_init(&L, alpha, 1);
2115 14 : u = get_u(E, f, L.N, L.prec);
2116 14 : A = cgetg(LIM+2, t_VEC);
2117 217 : for(i = 0; i <= LIM; i++)
2118 : {
2119 203 : GEN a = RgV_dotproduct(u,L.coef);
2120 : long n;
2121 34461 : for (n = 1; n <= L.N; n++)
2122 34258 : gel(u,n) = gprec_wensure(gmul(gsub(gel(u,n), a), gel(L.na,n)), L.prec);
2123 203 : gel(A,i+1) = gprec_wtrunc(a, prec);
2124 : }
2125 14 : return gerepilecopy(av, A);
2126 : }
2127 : GEN
2128 56 : asympnum0(GEN u, GEN alpha, long prec)
2129 56 : { return asympnum((void*)u,get_fun(u, "asympnum"), alpha, prec); }
2130 : GEN
2131 14 : asympnumraw0(GEN u, long LIM, GEN alpha, long prec)
2132 14 : { return asympnumraw((void*)u,get_fun(u, "asympnumraw"), LIM, alpha, prec); }
|
