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