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