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