Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - language - sumiter.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24038-ebe36f6c4) Lines: 1106 1143 96.8 %
Date: 2019-07-23 05:53:17 Functions: 93 94 98.9 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13