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

Generated by: LCOV version 1.13