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 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.10.0 lcov report (development 21196-f12677d) Lines: 970 1011 95.9 %
Date: 2017-10-22 06:23:24 Functions: 79 79 100.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     1270929 : iferrpari(GEN a, GEN b, GEN c)
      19             : {
      20             :   GEN res;
      21             :   struct pari_evalstate state;
      22     1270929 :   evalstate_save(&state);
      23     1270929 :   pari_CATCH(CATCH_ALL)
      24             :   {
      25             :     GEN E;
      26       67294 :     if (!b&&!c) return gnil;
      27       33654 :     E = evalstate_restore_err(&state);
      28       33654 :     if (c)
      29             :     {
      30         236 :       push_lex(E,c);
      31         236 :       res = closure_evalnobrk(c);
      32         229 :       pop_lex(1);
      33         229 :       if (gequal0(res))
      34           7 :         pari_err(0, E);
      35             :     }
      36       33640 :     if (!b) return gnil;
      37       33640 :     push_lex(E,b);
      38       33640 :     res = closure_evalgen(b);
      39       33640 :     pop_lex(1);
      40       33640 :     return res;
      41             :   } pari_TRY {
      42     1270929 :     res = closure_evalgen(a);
      43     1237275 :   } pari_ENDCATCH;
      44     1237275 :   return res;
      45             : }
      46             : 
      47             : /********************************************************************/
      48             : /**                                                                **/
      49             : /**                        ITERATIONS                              **/
      50             : /**                                                                **/
      51             : /********************************************************************/
      52             : 
      53             : static void
      54     4218370 : forparii(GEN a, GEN b, GEN code)
      55             : {
      56     4218370 :   pari_sp av, av0 = avma;
      57             :   GEN aa;
      58     8436684 :   if (gcmp(b,a) < 0) return;
      59     4164022 :   if (typ(b) != t_INFINITY) b = gfloor(b);
      60     4164013 :   aa = a = setloop(a);
      61     4164022 :   av=avma;
      62     4164022 :   push_lex(a,code);
      63    87509575 :   while (gcmp(a,b) <= 0)
      64             :   {
      65    79503313 :     closure_evalvoid(code); if (loop_break()) break;
      66    79202434 :     a = get_lex(-1);
      67    79168136 :     if (a == aa)
      68             :     {
      69    79168108 :       a = incloop(a);
      70    79181509 :       if (a != aa) { set_lex(-1,a); aa = a; }
      71             :     }
      72             :     else
      73             :     { /* 'code' modified a ! Be careful (and slow) from now on */
      74          28 :       a = gaddgs(a,1);
      75          28 :       if (gc_needed(av,1))
      76             :       {
      77           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"forparii");
      78           0 :         a = gerepileupto(av,a);
      79             :       }
      80          28 :       set_lex(-1,a);
      81             :     }
      82             :   }
      83     4101859 :   pop_lex(1);  avma = av0;
      84             : }
      85             : 
      86             : void
      87     4218377 : forpari(GEN a, GEN b, GEN code)
      88             : {
      89     4218377 :   pari_sp ltop=avma, av;
      90     8436694 :   if (typ(a) == t_INT) { forparii(a,b,code); return; }
      91           7 :   b = gcopy(b); /* Kludge to work-around the a+(a=2) bug */
      92           7 :   av=avma;
      93           7 :   push_lex(a,code);
      94          35 :   while (gcmp(a,b) <= 0)
      95             :   {
      96          21 :     closure_evalvoid(code); if (loop_break()) break;
      97          21 :     a = get_lex(-1); a = gaddgs(a,1);
      98          21 :     if (gc_needed(av,1))
      99             :     {
     100           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"forpari");
     101           0 :       a = gerepileupto(av,a);
     102             :     }
     103          21 :     set_lex(-1, a);
     104             :   }
     105           7 :   pop_lex(1); avma = ltop;
     106             : }
     107             : 
     108             : /* 0 < a <= b. Using small consecutive chunks to 1) limit memory use, 2) allow
     109             :  * cheap early abort */
     110             : static int
     111          49 : forfactoredpos(ulong a, ulong b, GEN code)
     112             : {
     113          49 :   const ulong step = 1024;
     114          49 :   pari_sp av = avma;
     115             :   ulong x1;
     116        6874 :   for(x1 = a;; x1 += step, avma = av)
     117             :   { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
     118        6874 :     ulong j, lv, x2 = (b >= 2*step && b - 2*step >= x1)? x1-1 + step: b;
     119        6874 :     GEN v = vecfactoru(x1, x2);
     120        6874 :     lv = lg(v);
     121     7007679 :     for (j = 1; j < lv; j++)
     122             :     {
     123     7000819 :       ulong n = x1-1 + j;
     124     7000819 :       set_lex(-1, mkvec2(utoipos(n), Flm_to_ZM(gel(v,j))));
     125     7000819 :       closure_evalvoid(code);
     126     7000819 :       if (loop_break()) return 1;
     127             :     }
     128        6860 :     if (x2 == b) break;
     129        6825 :     set_lex(-1, gen_0);
     130        6825 :   }
     131          35 :   return 0;
     132             : }
     133             : /* convert factoru(n) to factor(-n); M pre-allocated factorization matrix
     134             :  * with (-1)^1 already set */
     135             : static void
     136     7000882 : Flm2negfact(GEN v, GEN M)
     137             : {
     138     7000882 :   GEN p = gel(v,1), e = gel(v,2), P = gel(M,1), E = gel(M,2);
     139     7000882 :   long i, l = lg(p);
     140    26978161 :   for (i = 1; i < l; i++)
     141             :   {
     142    19977279 :     gel(P,i+1) = utoipos(p[i]);
     143    19977279 :     gel(E,i+1) = utoipos(e[i]);
     144             :   }
     145     7000882 :   setlg(P,l+1);
     146     7000882 :   setlg(E,l+1);
     147     7000882 : }
     148             : /* 0 < a <= b, from -b to -a */
     149             : static int
     150          77 : forfactoredneg(ulong a, ulong b, GEN code)
     151             : {
     152          77 :   const ulong step = 1024;
     153             :   GEN P, E, M;
     154             :   pari_sp av;
     155             :   ulong x2;
     156             : 
     157          77 :   P = cgetg(18, t_COL); gel(P,1) = gen_m1;
     158          77 :   E = cgetg(18, t_COL); gel(E,1) = gen_1;
     159          77 :   M = mkmat2(P,E);
     160          77 :   av = avma;
     161        6902 :   for(x2 = b;; x2 -= step, avma = av)
     162             :   { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
     163        6902 :     ulong j, x1 = (x2 >= 2*step && x2-2*step >= a)? x2+1 - step: a;
     164        6902 :     GEN v = vecfactoru(x1, x2);
     165     7007763 :     for (j = lg(v)-1; j; j--)
     166             :     { /* run backward: from factor(x1..x2) to factor(-x2..-x1) */
     167     7000882 :       ulong n = x1-1 + j;
     168     7000882 :       Flm2negfact(gel(v,j), M);
     169     7000882 :       set_lex(-1, mkvec2(utoineg(n), M));
     170     7000882 :       closure_evalvoid(code);
     171     7000882 :       if (loop_break()) return 1;
     172             :     }
     173        6881 :     if (x1 == a) break;
     174        6825 :     set_lex(-1, gen_0);
     175        6825 :   }
     176          56 :   return 0;
     177             : }
     178             : static int
     179          63 : eval0(GEN code)
     180             : {
     181          63 :   pari_sp av = avma;
     182          63 :   set_lex(-1, mkvec2(gen_0, mkmat2(mkcol(gen_0),mkcol(gen_1))));
     183          63 :   closure_evalvoid(code); avma = av;
     184          63 :   return loop_break();
     185             : }
     186             : void
     187         119 : forfactored(GEN a, GEN b, GEN code)
     188             : {
     189         119 :   pari_sp av = avma;
     190         119 :   long sa, sb, stop = 0;
     191         119 :   if (typ(a) != t_INT) pari_err_TYPE("forfactored", a);
     192         119 :   if (typ(b) != t_INT) pari_err_TYPE("forfactored", b);
     193         238 :   if (cmpii(a,b) > 0) return;
     194         112 :   push_lex(NULL,code);
     195         112 :   sa = signe(a);
     196         112 :   sb = signe(b);
     197         112 :   if (sa < 0)
     198             :   {
     199          77 :     stop = forfactoredneg((sb < 0)? b[2]: 1UL, itou(a), code);
     200          77 :     if (!stop && sb >= 0) stop = eval0(code);
     201          77 :     if (!stop && sb > 0) forfactoredpos(1UL, b[2], code);
     202             :   }
     203             :   else
     204             :   {
     205          35 :     if (!sa) stop = eval0(code);
     206          35 :     if (!stop && sb) forfactoredpos(sa? a[2]: 1UL, itou(b), code);
     207             :   }
     208         112 :   pop_lex(1); avma = av;
     209             : }
     210             : void
     211     1768823 : whilepari(GEN a, GEN b)
     212             : {
     213     1768823 :   pari_sp av = avma;
     214             :   for(;;)
     215             :   {
     216    18638151 :     GEN res = closure_evalnobrk(a);
     217    18638151 :     if (gequal0(res)) break;
     218    16869328 :     avma = av;
     219    16869328 :     closure_evalvoid(b); if (loop_break()) break;
     220    16869328 :   }
     221     1768823 :   avma = av;
     222     1768823 : }
     223             : 
     224             : void
     225      222074 : untilpari(GEN a, GEN b)
     226             : {
     227      222074 :   pari_sp av = avma;
     228             :   for(;;)
     229             :   {
     230             :     GEN res;
     231     2194143 :     closure_evalvoid(b); if (loop_break()) break;
     232     2194143 :     res = closure_evalnobrk(a);
     233     2194143 :     if (!gequal0(res)) break;
     234     1972069 :     avma = av;
     235     1972069 :   }
     236      222074 :   avma = av;
     237      222074 : }
     238             : 
     239          28 : static int negcmp(GEN x, GEN y) { return gcmp(y,x); }
     240             : 
     241             : void
     242        1484 : forstep(GEN a, GEN b, GEN s, GEN code)
     243             : {
     244             :   long ss, i;
     245        1484 :   pari_sp av, av0 = avma;
     246        1484 :   GEN v = NULL;
     247             :   int (*cmp)(GEN,GEN);
     248             : 
     249        1484 :   b = gcopy(b); s = gcopy(s); av=avma;
     250        1484 :   push_lex(a,code);
     251        1484 :   if (is_vec_t(typ(s)))
     252             :   {
     253           7 :     v = s; s = gen_0;
     254           7 :     for (i=lg(v)-1; i; i--) s = gadd(s,gel(v,i));
     255             :   }
     256        1484 :   ss = gsigne(s);
     257        1484 :   if (!ss) pari_err_DOMAIN("forstep","step","=",gen_0,s);
     258        1477 :   cmp = (ss > 0)? &gcmp: &negcmp;
     259        1477 :   i = 0;
     260       50484 :   while (cmp(a,b) <= 0)
     261             :   {
     262       47530 :     closure_evalvoid(code); if (loop_break()) break;
     263       47530 :     if (v)
     264             :     {
     265          42 :       if (++i >= lg(v)) i = 1;
     266          42 :       s = gel(v,i);
     267             :     }
     268       47530 :     a = get_lex(-1); a = gadd(a,s);
     269             : 
     270       47530 :     if (gc_needed(av,1))
     271             :     {
     272           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"forstep");
     273           0 :       a = gerepileupto(av,a);
     274             :     }
     275       47530 :     set_lex(-1,a);
     276             :   }
     277        1477 :   pop_lex(1); avma = av0;
     278        1477 : }
     279             : 
     280             : static void
     281          14 : _fordiv(GEN a, GEN code, GEN (*D)(GEN))
     282             : {
     283             :   long i, l;
     284          14 :   pari_sp av2, av = avma;
     285          14 :   GEN t = D(a);
     286          14 :   push_lex(gen_0,code); l=lg(t); av2 = avma;
     287         105 :   for (i=1; i<l; i++)
     288             :   {
     289          91 :     set_lex(-1,gel(t,i));
     290          91 :     closure_evalvoid(code); if (loop_break()) break;
     291          91 :     avma = av2;
     292             :   }
     293          14 :   pop_lex(1); avma=av;
     294          14 : }
     295             : void
     296           7 : fordiv(GEN a, GEN code) { return _fordiv(a, code, &divisors); }
     297             : void
     298           7 : fordivfactored(GEN a, GEN code) { return _fordiv(a, code, &divisors_factored); }
     299             : 
     300             : /* Embedded for loops:
     301             :  *   fl = 0: execute ch (a), where a = (ai) runs through all n-uplets in
     302             :  *     [m1,M1] x ... x [mn,Mn]
     303             :  *   fl = 1: impose a1 <= ... <= an
     304             :  *   fl = 2:        a1 <  ... <  an
     305             :  */
     306             : /* increment and return d->a [over integers]*/
     307             : static GEN
     308      181472 : _next_i(forvec_t *d)
     309             : {
     310      181472 :   long i = d->n;
     311      181472 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     312             :   for (;;) {
     313      233116 :     if (cmpii(d->a[i], d->M[i]) < 0) {
     314      181236 :       d->a[i] = incloop(d->a[i]);
     315      181236 :       return (GEN)d->a;
     316             :     }
     317       51880 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     318       51880 :     if (--i <= 0) return NULL;
     319       51762 :   }
     320             : }
     321             : /* increment and return d->a [generic]*/
     322             : static GEN
     323          63 : _next(forvec_t *d)
     324             : {
     325          63 :   long i = d->n;
     326          63 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     327             :   for (;;) {
     328          98 :     d->a[i] = gaddgs(d->a[i], 1);
     329          98 :     if (gcmp(d->a[i], d->M[i]) <= 0) return (GEN)d->a;
     330          49 :     d->a[i] = d->m[i];
     331          49 :     if (--i <= 0) return NULL;
     332          42 :   }
     333             : }
     334             : 
     335             : /* non-decreasing order [over integers] */
     336             : static GEN
     337         157 : _next_le_i(forvec_t *d)
     338             : {
     339         157 :   long i = d->n;
     340         157 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     341             :   for (;;) {
     342         231 :     if (cmpii(d->a[i], d->M[i]) < 0)
     343             :     {
     344         117 :       d->a[i] = incloop(d->a[i]);
     345             :       /* m[i] < a[i] <= M[i] <= M[i+1] */
     346         301 :       while (i < d->n)
     347             :       {
     348             :         GEN t;
     349          67 :         i++;
     350          67 :         if (cmpii(d->a[i-1], d->a[i]) <= 0) continue;
     351             :         /* a[i] < a[i-1] <= M[i-1] <= M[i] */
     352          67 :         t = d->a[i-1]; if (cmpii(t, d->m[i]) < 0) t = d->m[i];
     353          67 :         d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1],m[i])*/
     354             :       }
     355         117 :       return (GEN)d->a;
     356             :     }
     357         114 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     358         114 :     if (--i <= 0) return NULL;
     359          94 :   }
     360             : }
     361             : /* non-decreasing order [generic] */
     362             : static GEN
     363         154 : _next_le(forvec_t *d)
     364             : {
     365         154 :   long i = d->n;
     366         154 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     367             :   for (;;) {
     368         266 :     d->a[i] = gaddgs(d->a[i], 1);
     369         266 :     if (gcmp(d->a[i], d->M[i]) <= 0)
     370             :     {
     371         350 :       while (i < d->n)
     372             :       {
     373             :         GEN c;
     374          98 :         i++;
     375          98 :         if (gcmp(d->a[i-1], d->a[i]) <= 0) continue;
     376             :         /* M[i] >= M[i-1] >= a[i-1] > a[i] */
     377          98 :         c = gceil(gsub(d->a[i-1], d->a[i]));
     378          98 :         d->a[i] = gadd(d->a[i], c);
     379             :         /* a[i-1] <= a[i] < M[i-1] + 1 => a[i] < M[i]+1 => a[i] <= M[i] */
     380             :       }
     381         126 :       return (GEN)d->a;
     382             :     }
     383         140 :     d->a[i] = d->m[i];
     384         140 :     if (--i <= 0) return NULL;
     385         126 :   }
     386             : }
     387             : /* strictly increasing order [over integers] */
     388             : static GEN
     389     1173546 : _next_lt_i(forvec_t *d)
     390             : {
     391     1173546 :   long i = d->n;
     392     1173546 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     393             :   for (;;) {
     394     1290065 :     if (cmpii(d->a[i], d->M[i]) < 0)
     395             :     {
     396     1159940 :       d->a[i] = incloop(d->a[i]);
     397             :       /* m[i] < a[i] <= M[i] < M[i+1] */
     398     2436385 :       while (i < d->n)
     399             :       {
     400             :         pari_sp av;
     401             :         GEN t;
     402      116505 :         i++;
     403      116505 :         if (cmpii(d->a[i-1], d->a[i]) < 0) continue;
     404      116505 :         av = avma;
     405             :         /* M[i] > M[i-1] >= a[i-1] */
     406      116505 :         t = addiu(d->a[i-1],1); if (cmpii(t, d->m[i]) < 0) t = d->m[i];
     407      116505 :         d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1]+1,m[i]) <= M[i]*/
     408      116505 :         avma = av;
     409             :       }
     410     1159940 :       return (GEN)d->a;
     411             :     }
     412      130125 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     413      130125 :     if (--i <= 0) return NULL;
     414      123322 :   }
     415             : }
     416             : /* strictly increasing order [generic] */
     417             : static GEN
     418          84 : _next_lt(forvec_t *d)
     419             : {
     420          84 :   long i = d->n;
     421          84 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     422             :   for (;;) {
     423         133 :     d->a[i] = gaddgs(d->a[i], 1);
     424         133 :     if (gcmp(d->a[i], d->M[i]) <= 0)
     425             :     {
     426         147 :       while (i < d->n)
     427             :       {
     428             :         GEN c;
     429          35 :         i++;
     430          35 :         if (gcmp(d->a[i-1], d->a[i]) < 0) continue;
     431             :         /* M[i] > M[i-1] >= a[i-1] >= a[i] */
     432          35 :         c = addiu(gfloor(gsub(d->a[i-1], d->a[i])), 1); /* > a[i-1] - a[i] */
     433          35 :         d->a[i] = gadd(d->a[i], c);
     434             :         /* a[i-1] < a[i] <= M[i-1] + 1 => a[i] < M[i]+1 => a[i] <= M[i] */
     435             :       }
     436          56 :       return (GEN)d->a;
     437             :     }
     438          77 :     d->a[i] = d->m[i];
     439          77 :     if (--i <= 0) return NULL;
     440          63 :   }
     441             : }
     442             : /* for forvec(v=[],) */
     443             : static GEN
     444          14 : _next_void(forvec_t *d)
     445             : {
     446          14 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     447           7 :   return NULL;
     448             : }
     449             : 
     450             : /* Initialize minima (m) and maxima (M); guarantee M[i] - m[i] integer and
     451             :  *   if flag = 1: m[i-1] <= m[i] <= M[i] <= M[i+1]
     452             :  *   if flag = 2: m[i-1] <  m[i] <= M[i] <  M[i+1],
     453             :  * for all i */
     454             : int
     455        7018 : forvec_init(forvec_t *d, GEN x, long flag)
     456             : {
     457        7018 :   long i, tx = typ(x), l = lg(x), t = t_INT;
     458        7018 :   if (!is_vec_t(tx)) pari_err_TYPE("forvec [not a vector]", x);
     459        7018 :   d->first = 1;
     460        7018 :   d->n = l - 1;
     461        7018 :   d->a = (GEN*)cgetg(l,tx);
     462        7018 :   d->m = (GEN*)cgetg(l,tx);
     463        7018 :   d->M = (GEN*)cgetg(l,tx);
     464        7018 :   if (l == 1) { d->next = &_next_void; return 1; }
     465       21285 :   for (i = 1; i < l; i++)
     466             :   {
     467       14302 :     GEN a, e = gel(x,i), m = gel(e,1), M = gel(e,2);
     468       14302 :     tx = typ(e);
     469       14302 :     if (! is_vec_t(tx) || lg(e)!=3)
     470          14 :       pari_err_TYPE("forvec [expected vector not of type [min,MAX]]",e);
     471       14288 :     if (typ(m) != t_INT) t = t_REAL;
     472       14288 :     if (i > 1) switch(flag)
     473             :     {
     474             :       case 1: /* a >= m[i-1] - m */
     475          55 :         a = gceil(gsub(d->m[i-1], m));
     476          55 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     477          55 :         if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
     478          55 :         break;
     479             :       case 2: /* a > m[i-1] - m */
     480        6852 :         a = gfloor(gsub(d->m[i-1], m));
     481        6852 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     482        6852 :         a = addiu(a, 1);
     483        6852 :         if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
     484        6852 :         break;
     485         384 :       default: m = gcopy(m);
     486         384 :         break;
     487             :     }
     488       14288 :     M = gadd(m, gfloor(gsub(M,m))); /* ensure M-m is an integer */
     489       14281 :     if (gcmp(m,M) > 0) { d->a = NULL; d->next = &_next; return 0; }
     490       14274 :     d->m[i] = m;
     491       14274 :     d->M[i] = M;
     492             :   }
     493        7038 :   if (flag == 1) for (i = l-2; i >= 1; i--)
     494             :   {
     495          55 :     GEN M = d->M[i], a = gfloor(gsub(d->M[i+1], M));
     496          55 :     if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     497             :     /* M[i]+a <= M[i+1] */
     498          55 :     if (signe(a) < 0) d->M[i] = gadd(M, a);
     499             :   }
     500       13794 :   else if (flag == 2) for (i = l-2; i >= 1; i--)
     501             :   {
     502        6845 :     GEN M = d->M[i], a = gceil(gsub(d->M[i+1], M));
     503        6845 :     if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     504        6845 :     a = subiu(a, 1);
     505             :     /* M[i]+a < M[i+1] */
     506        6845 :     if (signe(a) < 0) d->M[i] = gadd(M, a);
     507             :   }
     508        6983 :   if (t == t_INT) {
     509       21110 :     for (i = 1; i < l; i++) {
     510       14162 :       d->a[i] = setloop(d->m[i]);
     511       14162 :       if (typ(d->M[i]) != t_INT) d->M[i] = gfloor(d->M[i]);
     512             :     }
     513             :   } else {
     514          35 :     for (i = 1; i < l; i++) d->a[i] = d->m[i];
     515             :   }
     516        6983 :   switch(flag)
     517             :   {
     518         125 :     case 0: d->next = t==t_INT? &_next_i:    &_next; break;
     519          34 :     case 1: d->next = t==t_INT? &_next_le_i: &_next_le; break;
     520        6817 :     case 2: d->next = t==t_INT? &_next_lt_i: &_next_lt; break;
     521           7 :     default: pari_err_FLAG("forvec");
     522             :   }
     523        6976 :   return 1;
     524             : }
     525             : GEN
     526     1355490 : forvec_next(forvec_t *d) { return d->next(d); }
     527             : 
     528             : void
     529        7000 : forvec(GEN x, GEN code, long flag)
     530             : {
     531        7000 :   pari_sp av = avma;
     532             :   forvec_t T;
     533             :   GEN v;
     534       13972 :   if (!forvec_init(&T, x, flag)) { avma = av; return; }
     535        6965 :   push_lex((GEN)T.a, code);
     536        6965 :   while ((v = forvec_next(&T)))
     537             :   {
     538     1348333 :     closure_evalvoid(code);
     539     1348333 :     if (loop_break()) break;
     540             :   }
     541        6965 :   pop_lex(1); avma = av;
     542             : }
     543             : 
     544             : /********************************************************************/
     545             : /**                                                                **/
     546             : /**                              SUMS                              **/
     547             : /**                                                                **/
     548             : /********************************************************************/
     549             : 
     550             : GEN
     551       63392 : somme(GEN a, GEN b, GEN code, GEN x)
     552             : {
     553       63392 :   pari_sp av, av0 = avma;
     554             :   GEN p1;
     555             : 
     556       63392 :   if (typ(a) != t_INT) pari_err_TYPE("sum",a);
     557       63392 :   if (!x) x = gen_0;
     558       63392 :   if (gcmp(b,a) < 0) return gcopy(x);
     559             : 
     560       63392 :   b = gfloor(b);
     561       63392 :   a = setloop(a);
     562       63392 :   av=avma;
     563       63392 :   push_lex(a,code);
     564             :   for(;;)
     565             :   {
     566     1716400 :     p1 = closure_evalnobrk(code);
     567     1716400 :     x=gadd(x,p1); if (cmpii(a,b) >= 0) break;
     568     1653008 :     a = incloop(a);
     569     1653008 :     if (gc_needed(av,1))
     570             :     {
     571           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sum");
     572           0 :       x = gerepileupto(av,x);
     573             :     }
     574     1653008 :     set_lex(-1,a);
     575     1653008 :   }
     576       63392 :   pop_lex(1); return gerepileupto(av0,x);
     577             : }
     578             : 
     579             : static GEN
     580          21 : sum_init(GEN x0, GEN t)
     581             : {
     582          21 :   long tp = typ(t);
     583             :   GEN x;
     584          21 :   if (is_vec_t(tp))
     585             :   {
     586           7 :     x = const_vec(lg(t)-1, x0);
     587           7 :     settyp(x, tp);
     588             :   }
     589             :   else
     590          14 :     x = x0;
     591          21 :   return x;
     592             : }
     593             : 
     594             : GEN
     595          21 : suminf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     596             : {
     597          21 :   long fl = 0, G = prec2nbits(prec) + 5;
     598          21 :   pari_sp av0 = avma, av;
     599          21 :   GEN x = NULL, _1;
     600             : 
     601          21 :   if (typ(a) != t_INT) pari_err_TYPE("suminf",a);
     602          21 :   a = setloop(a);
     603          21 :   av = avma;
     604             :   for(;;)
     605             :   {
     606       15379 :     GEN t = eval(E, a);
     607       15379 :     if (!x) _1 = x = sum_init(real_1(prec), t);
     608             : 
     609       15379 :     x = gadd(x,t);
     610       15379 :     if (!gequal0(t) && gexpo(t) > gexpo(x)-G)
     611       15316 :       fl = 0;
     612          63 :     else if (++fl == 3)
     613          21 :       break;
     614       15358 :     a = incloop(a);
     615       15358 :     if (gc_needed(av,1))
     616             :     {
     617           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"suminf");
     618           0 :       gerepileall(av,2, &x, &_1);
     619             :     }
     620       15358 :   }
     621          21 :   return gerepileupto(av0, gsub(x, _1));
     622             : }
     623             : GEN
     624          21 : suminf0(GEN a, GEN code, long prec)
     625          21 : { EXPR_WRAP(code, suminf(EXPR_ARG, a, prec)); }
     626             : 
     627             : GEN
     628          49 : sumdivexpr(GEN num, GEN code)
     629             : {
     630          49 :   pari_sp av = avma;
     631          49 :   GEN y = gen_0, t = divisors(num);
     632          49 :   long i, l = lg(t);
     633             : 
     634          49 :   push_lex(gen_0, code);
     635        9289 :   for (i=1; i<l; i++)
     636             :   {
     637        9240 :     set_lex(-1,gel(t,i));
     638        9240 :     y = gadd(y, closure_evalnobrk(code));
     639             :   }
     640          49 :   pop_lex(1); return gerepileupto(av,y);
     641             : }
     642             : GEN
     643          42 : sumdivmultexpr(GEN num, GEN code)
     644             : {
     645          42 :   pari_sp av = avma;
     646          42 :   GEN y = gen_1, P,E;
     647          42 :   int isint = divisors_init(num, &P,&E);
     648          42 :   long i, l = lg(P);
     649             :   GEN (*mul)(GEN,GEN);
     650             : 
     651          42 :   if (l == 1) { avma = av; return gen_1; }
     652          42 :   push_lex(gen_0, code);
     653          42 :   mul = isint? mulii: gmul;
     654         196 :   for (i=1; i<l; i++)
     655             :   {
     656         154 :     GEN p = gel(P,i), q = p, z = gen_1;
     657         154 :     long j, e = E[i];
     658         560 :     for (j = 1; j <= e; j++, q = mul(q, p))
     659             :     {
     660         560 :       set_lex(-1, q);
     661         560 :       z = gadd(z, closure_evalnobrk(code));
     662         560 :       if (j == e) break;
     663             :     }
     664         154 :     y = gmul(y, z);
     665             :   }
     666          42 :   pop_lex(1); return gerepileupto(av,y);
     667             : }
     668             : 
     669             : /********************************************************************/
     670             : /**                                                                **/
     671             : /**                           PRODUCTS                             **/
     672             : /**                                                                **/
     673             : /********************************************************************/
     674             : 
     675             : GEN
     676      120148 : produit(GEN a, GEN b, GEN code, GEN x)
     677             : {
     678      120148 :   pari_sp av, av0 = avma;
     679             :   GEN p1;
     680             : 
     681      120148 :   if (typ(a) != t_INT) pari_err_TYPE("prod",a);
     682      120148 :   if (!x) x = gen_1;
     683      120148 :   if (gcmp(b,a) < 0) return gcopy(x);
     684             : 
     685      115206 :   b = gfloor(b);
     686      115206 :   a = setloop(a);
     687      115206 :   av=avma;
     688      115206 :   push_lex(a,code);
     689             :   for(;;)
     690             :   {
     691      348768 :     p1 = closure_evalnobrk(code);
     692      348768 :     x = gmul(x,p1); if (cmpii(a,b) >= 0) break;
     693      233562 :     a = incloop(a);
     694      233562 :     if (gc_needed(av,1))
     695             :     {
     696           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prod");
     697           0 :       x = gerepileupto(av,x);
     698             :     }
     699      233562 :     set_lex(-1,a);
     700      233562 :   }
     701      115206 :   pop_lex(1); return gerepileupto(av0,x);
     702             : }
     703             : 
     704             : GEN
     705           7 : prodinf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     706             : {
     707           7 :   pari_sp av0 = avma, av;
     708             :   long fl,G;
     709           7 :   GEN p1,x = real_1(prec);
     710             : 
     711           7 :   if (typ(a) != t_INT) pari_err_TYPE("prodinf",a);
     712           7 :   a = setloop(a);
     713           7 :   av = avma;
     714           7 :   fl=0; G = -prec2nbits(prec)-5;
     715             :   for(;;)
     716             :   {
     717         952 :     p1 = eval(E, a); if (gequal0(p1)) { x = p1; break; }
     718         952 :     x = gmul(x,p1); a = incloop(a);
     719         952 :     p1 = gsubgs(p1, 1);
     720         952 :     if (gequal0(p1) || gexpo(p1) <= G) { if (++fl==3) break; } else fl=0;
     721         945 :     if (gc_needed(av,1))
     722             :     {
     723           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodinf");
     724           0 :       x = gerepileupto(av,x);
     725             :     }
     726         945 :   }
     727           7 :   return gerepilecopy(av0,x);
     728             : }
     729             : GEN
     730           7 : prodinf1(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     731             : {
     732           7 :   pari_sp av0 = avma, av;
     733             :   long fl,G;
     734           7 :   GEN p1,p2,x = real_1(prec);
     735             : 
     736           7 :   if (typ(a) != t_INT) pari_err_TYPE("prodinf1",a);
     737           7 :   a = setloop(a);
     738           7 :   av = avma;
     739           7 :   fl=0; G = -prec2nbits(prec)-5;
     740             :   for(;;)
     741             :   {
     742         952 :     p2 = eval(E, a); p1 = gaddgs(p2,1);
     743         952 :     if (gequal0(p1)) { x = p1; break; }
     744         952 :     x = gmul(x,p1); a = incloop(a);
     745         952 :     if (gequal0(p2) || gexpo(p2) <= G) { if (++fl==3) break; } else fl=0;
     746         945 :     if (gc_needed(av,1))
     747             :     {
     748           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodinf1");
     749           0 :       x = gerepileupto(av,x);
     750             :     }
     751         945 :   }
     752           7 :   return gerepilecopy(av0,x);
     753             : }
     754             : GEN
     755          14 : prodinf0(GEN a, GEN code, long flag, long prec)
     756             : {
     757          14 :   switch(flag)
     758             :   {
     759           7 :     case 0: EXPR_WRAP(code, prodinf (EXPR_ARG, a, prec));
     760           7 :     case 1: EXPR_WRAP(code, prodinf1(EXPR_ARG, a, prec));
     761             :   }
     762           0 :   pari_err_FLAG("prodinf");
     763             :   return NULL; /* LCOV_EXCL_LINE */
     764             : }
     765             : 
     766             : GEN
     767           7 : prodeuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
     768             : {
     769           7 :   pari_sp av, av0 = avma;
     770           7 :   GEN x = real_1(prec), prime;
     771             :   forprime_t T;
     772             : 
     773           7 :   av = avma;
     774           7 :   if (!forprime_init(&T, a,b)) { avma = av; return x; }
     775             : 
     776           7 :   av = avma;
     777        8617 :   while ( (prime = forprime_next(&T)) )
     778             :   {
     779        8603 :     x = gmul(x, eval(E, prime));
     780        8603 :     if (gc_needed(av,1))
     781             :     {
     782           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodeuler");
     783           0 :       x = gerepilecopy(av, x);
     784             :     }
     785             :   }
     786           7 :   return gerepilecopy(av0,x);
     787             : }
     788             : GEN
     789           7 : prodeuler0(GEN a, GEN b, GEN code, long prec)
     790           7 : { EXPR_WRAP(code, prodeuler(EXPR_ARG, a, b, prec)); }
     791             : GEN
     792         126 : direuler0(GEN a, GEN b, GEN code, GEN c)
     793         126 : { EXPR_WRAP(code, direuler(EXPR_ARG, a, b, c)); }
     794             : 
     795             : /********************************************************************/
     796             : /**                                                                **/
     797             : /**                       VECTORS & MATRICES                       **/
     798             : /**                                                                **/
     799             : /********************************************************************/
     800             : 
     801             : INLINE GEN
     802     3376513 : copyupto(GEN z, GEN t)
     803             : {
     804     3376513 :   if (is_universal_constant(z) || (z>(GEN)pari_mainstack->bot && z<=t))
     805     3376506 :     return z;
     806             :   else
     807           7 :     return gcopy(z);
     808             : }
     809             : 
     810             : GEN
     811      114730 : vecexpr0(GEN vec, GEN code, GEN pred)
     812             : {
     813      114730 :   switch(typ(vec))
     814             :   {
     815             :     case t_LIST:
     816             :     {
     817          14 :       if (list_typ(vec)==t_LIST_MAP)
     818           0 :         vec = mapdomain_shallow(vec);
     819             :       else
     820          14 :         vec = list_data(vec);
     821          14 :       if (!vec) return cgetg(1, t_VEC);
     822           7 :       break;
     823             :     }
     824      114716 :     case t_VEC: case t_COL: case t_MAT: break;
     825           0 :     default: pari_err_TYPE("[_|_<-_,_]",vec);
     826             :   }
     827      114723 :   if (pred && code)
     828         259 :     EXPR_WRAP(code,vecselapply((void*)pred,&gp_evalbool,EXPR_ARGUPTO,vec))
     829      114464 :   else if (code)
     830      114464 :     EXPR_WRAP(code,vecapply(EXPR_ARGUPTO,vec))
     831             :   else
     832           0 :     EXPR_WRAP(pred,vecselect(EXPR_ARGBOOL,vec))
     833             : }
     834             : 
     835             : GEN
     836         175 : vecexpr1(GEN vec, GEN code, GEN pred)
     837             : {
     838         175 :   GEN v = vecexpr0(vec, code, pred);
     839         175 :   return lg(v) == 1? v: shallowconcat1(v);
     840             : }
     841             : 
     842             : GEN
     843     2325019 : vecteur(GEN nmax, GEN code)
     844             : {
     845             :   GEN y, c;
     846     2325019 :   long i, m = gtos(nmax);
     847             : 
     848     2325019 :   if (m < 0)  pari_err_DOMAIN("vector", "dimension", "<", gen_0, stoi(m));
     849     2325005 :   if (!code) return zerovec(m);
     850        5086 :   c = cgetipos(3); /* left on stack */
     851        5086 :   y = cgetg(m+1,t_VEC); push_lex(c, code);
     852     1625684 :   for (i=1; i<=m; i++)
     853             :   {
     854     1620612 :     c[2] = i;
     855     1620612 :     gel(y,i) = copyupto(closure_evalnobrk(code), y);
     856     1620598 :     set_lex(-1,c);
     857             :   }
     858        5072 :   pop_lex(1); return y;
     859             : }
     860             : 
     861             : GEN
     862         763 : vecteursmall(GEN nmax, GEN code)
     863             : {
     864             :   pari_sp av;
     865             :   GEN y, c;
     866         763 :   long i, m = gtos(nmax);
     867             : 
     868         763 :   if (m < 0)  pari_err_DOMAIN("vectorsmall", "dimension", "<", gen_0, stoi(m));
     869         756 :   if (!code) return zero_zv(m);
     870         735 :   c = cgetipos(3); /* left on stack */
     871         735 :   y = cgetg(m+1,t_VECSMALL); push_lex(c,code);
     872         735 :   av = avma;
     873       15974 :   for (i = 1; i <= m; i++)
     874             :   {
     875       15246 :     c[2] = i;
     876       15246 :     y[i] = gtos(closure_evalnobrk(code));
     877       15239 :     avma = av;
     878       15239 :     set_lex(-1,c);
     879             :   }
     880         728 :   pop_lex(1); return y;
     881             : }
     882             : 
     883             : GEN
     884         490 : vvecteur(GEN nmax, GEN n)
     885             : {
     886         490 :   GEN y = vecteur(nmax,n);
     887         483 :   settyp(y,t_COL); return y;
     888             : }
     889             : 
     890             : GEN
     891      136010 : matrice(GEN nlig, GEN ncol, GEN code)
     892             : {
     893             :   GEN c1, c2, y;
     894             :   long i, m, n;
     895             : 
     896      136010 :   n = gtos(nlig);
     897      136010 :   m = ncol? gtos(ncol): n;
     898      136010 :   if (m < 0)  pari_err_DOMAIN("matrix", "nbcols", "<", gen_0, stoi(m));
     899      136003 :   if (n < 0)  pari_err_DOMAIN("matrix", "nbrows", "<", gen_0, stoi(n));
     900      135996 :   if (!m) return cgetg(1,t_MAT);
     901      135954 :   if (!code || !n) return zeromatcopy(n, m);
     902      135856 :   c1 = cgetipos(3); push_lex(c1,code);
     903      135856 :   c2 = cgetipos(3); push_lex(c2,NULL); /* c1,c2 left on stack */
     904      135856 :   y = cgetg(m+1,t_MAT);
     905      524272 :   for (i = 1; i <= m; i++)
     906             :   {
     907      388416 :     GEN z = cgetg(n+1,t_COL);
     908             :     long j;
     909      388416 :     c2[2] = i; gel(y,i) = z;
     910     2144324 :     for (j = 1; j <= n; j++)
     911             :     {
     912     1755908 :       c1[2] = j;
     913     1755908 :       gel(z,j) = copyupto(closure_evalnobrk(code), y);
     914     1755908 :       set_lex(-2,c1);
     915     1755908 :       set_lex(-1,c2);
     916             :     }
     917             :   }
     918      135856 :   pop_lex(2); return y;
     919             : }
     920             : 
     921             : /********************************************************************/
     922             : /**                                                                **/
     923             : /**                         SUMMING SERIES                         **/
     924             : /**                                                                **/
     925             : /********************************************************************/
     926             : /* h = (2+2x)g'- g; g has t_INT coeffs */
     927             : static GEN
     928        1246 : delt(GEN g, long n)
     929             : {
     930        1246 :   GEN h = cgetg(n+3,t_POL);
     931             :   long k;
     932        1246 :   h[1] = g[1];
     933        1246 :   gel(h,2) = gel(g,2);
     934      359597 :   for (k=1; k<n; k++)
     935      358351 :     gel(h,k+2) = addii(mului(k+k+1,gel(g,k+2)), mului(k<<1,gel(g,k+1)));
     936        1246 :   gel(h,n+2) = mului(n<<1, gel(g,n+1)); return h;
     937             : }
     938             : 
     939             : #ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */
     940             : #pragma optimize("g",off)
     941             : #endif
     942             : /* P = polzagier(n,m)(-X), unnormalized (P(0) != 1) */
     943             : static GEN
     944          49 : polzag1(long n, long m)
     945             : {
     946          49 :   const long d = n - m, d2 = d<<1, r = (m+1)>>1, D = (d+1)>>1;
     947             :   long i, k;
     948          49 :   pari_sp av = avma;
     949             :   GEN g, T;
     950             : 
     951          49 :   if (d <= 0 || m < 0) return pol_0(0);
     952          49 :   g = cgetg(d+2, t_POL);
     953          49 :   g[1] = evalsigne(1)|evalvarn(0);
     954          49 :   T = cgetg(d+1,t_VEC);
     955             :   /* T[k+1] = binomial(2d,2k+1), 0 <= k < d */
     956          49 :   gel(T,1) = utoipos(d2);
     957        1267 :   for (k = 1; k < D; k++)
     958             :   {
     959        1218 :     long k2 = k<<1;
     960        2436 :     gel(T,k+1) = diviiexact(mulii(gel(T,k), muluu(d2-k2+1, d2-k2)),
     961        1218 :                             muluu(k2,k2+1));
     962             :   }
     963          49 :   for (; k < d; k++) gel(T,k+1) = gel(T,d-k);
     964          49 :   gel(g,2) = gel(T,d); /* binomial(2d, 2(d-1)+1) */
     965        2499 :   for (i = 1; i < d; i++)
     966             :   {
     967        2450 :     pari_sp av2 = avma;
     968        2450 :     GEN s, t = gel(T,d-i); /* binomial(2d, 2(d-1-i)+1) */
     969        2450 :     s = t;
     970      180082 :     for (k = d-i; k < d; k++)
     971             :     {
     972      177632 :       long k2 = k<<1;
     973      177632 :       t = diviiexact(mulii(t, muluu(d2-k2+1, d-k)), muluu(k2+1,k-(d-i)+1));
     974      177632 :       s = addii(s, t);
     975             :     }
     976             :     /* g_i = sum_{d-1-i <= k < d}, binomial(2*d, 2*k+1)*binomial(k,d-1-i) */
     977        2450 :     gel(g,i+2) = gerepileuptoint(av2, s);
     978             :   }
     979             :   /* sum_{0 <= i < d} g_i x^i * (x+x^2)^r */
     980          49 :   g = RgX_mulXn(gmul(g, gpowgs(deg1pol(gen_1,gen_1,0),r)), r);
     981          49 :   if (!odd(m)) g = delt(g, n);
     982        1274 :   for (i=1; i<=r; i++)
     983             :   {
     984        1225 :     g = delt(ZX_deriv(g), n);
     985        1225 :     if (gc_needed(av,4))
     986             :     {
     987           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"polzag, i = %ld/%ld", i,r);
     988           0 :       g = gerepilecopy(av, g);
     989             :     }
     990             :   }
     991          49 :   return g;
     992             : }
     993             : GEN
     994          28 : polzag(long n, long m)
     995             : {
     996          28 :   pari_sp av = avma;
     997          28 :   GEN g = ZX_z_unscale(polzag1(n,m), -1);
     998          28 :   return gerepileupto(av, RgX_Rg_div(g,gel(g,2)));
     999             : }
    1000             : 
    1001             : GEN
    1002           7 : sumalt(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1003             : {
    1004             :   ulong k, N;
    1005           7 :   pari_sp av = avma, av2;
    1006             :   GEN s, az, c, d;
    1007             : 
    1008           7 :   if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
    1009           7 :   N = (ulong)(0.39322*(prec2nbits(prec) + 7)); /*0.39322 > 1/log_2(3+sqrt(8))*/
    1010           7 :   d = powru(addsr(3, sqrtr(stor(8,prec))), N);
    1011           7 :   d = shiftr(addrr(d, invr(d)),-1);
    1012           7 :   a = setloop(a);
    1013           7 :   az = gen_m1; c = d;
    1014           7 :   s = gen_0;
    1015           7 :   av2 = avma;
    1016         371 :   for (k=0; ; k++) /* k < N */
    1017             :   {
    1018         371 :     c = addir(az,c); s = gadd(s, gmul(c, eval(E, a)));
    1019         371 :     if (k==N-1) break;
    1020         364 :     az = diviuuexact(muluui((N-k)<<1,N+k,az), k+1, (k<<1)+1);
    1021         364 :     a = incloop(a); /* in place! */
    1022         364 :     if (gc_needed(av,4))
    1023             :     {
    1024           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumalt, k = %ld/%ld", k,N-1);
    1025           0 :       gerepileall(av2, 3, &az,&c,&s);
    1026             :     }
    1027         364 :   }
    1028           7 :   return gerepileupto(av, gdiv(s,d));
    1029             : }
    1030             : 
    1031             : GEN
    1032           7 : sumalt2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1033             : {
    1034             :   long k, N;
    1035           7 :   pari_sp av = avma, av2;
    1036             :   GEN s, dn, pol;
    1037             : 
    1038           7 :   if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
    1039           7 :   N = (long)(0.307073*(prec2nbits(prec) + 5)); /*0.307073 > 1/log_2(\beta_B)*/
    1040           7 :   pol = ZX_div_by_X_1(polzag1(N,N>>1), &dn);
    1041           7 :   a = setloop(a);
    1042           7 :   N = degpol(pol);
    1043           7 :   s = gen_0;
    1044           7 :   av2 = avma;
    1045         280 :   for (k=0; k<=N; k++)
    1046             :   {
    1047         280 :     GEN t = itor(gel(pol,k+2), prec+EXTRAPRECWORD);
    1048         280 :     s = gadd(s, gmul(t, eval(E, a)));
    1049         280 :     if (k == N) break;
    1050         273 :     a = incloop(a); /* in place! */
    1051         273 :     if (gc_needed(av,4))
    1052             :     {
    1053           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumalt2, k = %ld/%ld", k,N-1);
    1054           0 :       s = gerepileupto(av2, s);
    1055             :     }
    1056             :   }
    1057           7 :   return gerepileupto(av, gdiv(s,dn));
    1058             : }
    1059             : 
    1060             : GEN
    1061          14 : sumalt0(GEN a, GEN code, long flag, long prec)
    1062             : {
    1063          14 :   switch(flag)
    1064             :   {
    1065           7 :     case 0: EXPR_WRAP(code, sumalt (EXPR_ARG,a,prec));
    1066           7 :     case 1: EXPR_WRAP(code, sumalt2(EXPR_ARG,a,prec));
    1067           0 :     default: pari_err_FLAG("sumalt");
    1068             :   }
    1069             :   return NULL; /* LCOV_EXCL_LINE */
    1070             : }
    1071             : 
    1072             : /* For k > 0, set S[k*2^i] <- g(k*2^i), k*2^i <= N = #S.
    1073             :  * Only needed with k odd (but also works for g even). */
    1074             : static void
    1075        7406 : binsum(GEN S, ulong k, void *E, GEN (*f)(void *, GEN), GEN a,
    1076             :         long G, long prec)
    1077             : {
    1078        7406 :   long e, i, N = lg(S)-1, l = expu(N / k); /* k 2^l <= N < k 2^(l+1) */
    1079             :   pari_sp av;
    1080        7406 :   GEN r, t = gen_0;
    1081             : 
    1082        7406 :   gel(S, k << l) = cgetr(prec); av = avma;
    1083        7406 :   G -= l;
    1084        7406 :   r = utoipos(k<<l);
    1085     3964471 :   for(e=0;;e++) /* compute g(k 2^l) with absolute error ~ 2^(G-l) */
    1086             :   {
    1087     3964471 :     GEN u = gtofp(f(E, addii(a,r)), prec);
    1088     3964471 :     if (typ(u) != t_REAL) pari_err_TYPE("sumpos",u);
    1089     3964471 :     if (!signe(u)) break;
    1090     3964282 :     if (!e)
    1091        7217 :       t = u;
    1092             :     else {
    1093     3957065 :       shiftr_inplace(u, e);
    1094     3957065 :       t = addrr(t,u);
    1095     3957065 :       if (expo(u) < G) break;
    1096             :     }
    1097     3957065 :     r = shifti(r,1);
    1098     3957065 :   }
    1099        7406 :   gel(S, k << l) = t = gerepileuptoleaf(av, t);
    1100             :   /* g(j) = 2g(2j) + f(a+j) for all j > 0 */
    1101       14812 :   for(i = l-1; i >= 0; i--)
    1102             :   { /* t ~ g(2 * k*2^i) with error ~ 2^(G-i-1) */
    1103             :     GEN u;
    1104        7406 :     av = avma; u = gtofp(f(E, addiu(a, k << i)), prec);
    1105        7406 :     if (typ(u) != t_REAL) pari_err_TYPE("sumpos",u);
    1106        7406 :     t = addrr(gtofp(u,prec), shiftr(t,1)); /* ~ g(k*2^i) */
    1107        7406 :     gel(S, k << i) = t = gerepileuptoleaf(av, t);
    1108             :   }
    1109        7406 : }
    1110             : /* For k > 0, let g(k) := \sum_{e >= 0} 2^e f(a + k*2^e).
    1111             :  * Return [g(k), 1 <= k <= N] */
    1112             : static GEN
    1113          70 : sumpos_init(void *E, GEN (*f)(void *, GEN), GEN a, long N, long prec)
    1114             : {
    1115          70 :   GEN S = cgetg(N+1,t_VEC);
    1116          70 :   long k, G = -prec2nbits(prec) - 5;
    1117          70 :   for (k=1; k<=N; k+=2) binsum(S,k, E,f, a,G,prec);
    1118          70 :   return S;
    1119             : }
    1120             : 
    1121             : GEN
    1122          56 : sumpos(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1123             : {
    1124             :   ulong k, N;
    1125          56 :   pari_sp av = avma;
    1126             :   GEN s, az, c, d, S;
    1127             : 
    1128          56 :   if (typ(a) != t_INT) pari_err_TYPE("sumpos",a);
    1129          56 :   a = subiu(a,1);
    1130          56 :   N = (ulong)(0.4*(prec2nbits(prec) + 7));
    1131          56 :   if (odd(N)) N++; /* extra precision for free */
    1132          56 :   d = powru(addsr(3, sqrtr(stor(8,prec))), N);
    1133          56 :   d = shiftr(addrr(d, invr(d)),-1);
    1134          56 :   az = gen_m1; c = d;
    1135             : 
    1136          56 :   S = sumpos_init(E, eval, a, N, prec);
    1137          56 :   s = gen_0;
    1138       10360 :   for (k=0; k<N; k++)
    1139             :   {
    1140             :     GEN t;
    1141       10360 :     c = addir(az,c);
    1142       10360 :     t = mulrr(gel(S,k+1), c);
    1143       10360 :     s = odd(k)? mpsub(s, t): mpadd(s, t);
    1144       10360 :     if (k == N-1) break;
    1145       10304 :     az = diviuuexact(muluui((N-k)<<1,N+k,az), k+1, (k<<1)+1);
    1146             :   }
    1147          56 :   return gerepileupto(av, gdiv(s,d));
    1148             : }
    1149             : 
    1150             : GEN
    1151          14 : sumpos2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1152             : {
    1153             :   ulong k, N;
    1154          14 :   pari_sp av = avma;
    1155             :   GEN s, pol, dn, S;
    1156             : 
    1157          14 :   if (typ(a) != t_INT) pari_err_TYPE("sumpos2",a);
    1158          14 :   a = subiu(a,1);
    1159          14 :   N = (ulong)(0.31*(prec2nbits(prec) + 5));
    1160             : 
    1161          14 :   if (odd(N)) N++; /* extra precision for free */
    1162          14 :   S = sumpos_init(E, eval, a, N, prec);
    1163          14 :   pol = ZX_div_by_X_1(polzag1(N,N>>1), &dn);
    1164          14 :   s = gen_0;
    1165        4466 :   for (k=0; k<N; k++)
    1166             :   {
    1167        4452 :     GEN t = mulri(gel(S,k+1), gel(pol,k+2));
    1168        4452 :     s = odd(k)? mpsub(s,t): mpadd(s,t);
    1169             :   }
    1170          14 :   return gerepileupto(av, gdiv(s,dn));
    1171             : }
    1172             : 
    1173             : GEN
    1174          70 : sumpos0(GEN a, GEN code, long flag, long prec)
    1175             : {
    1176          70 :   switch(flag)
    1177             :   {
    1178          56 :     case 0: EXPR_WRAP(code, sumpos (EXPR_ARG,a,prec));
    1179          14 :     case 1: EXPR_WRAP(code, sumpos2(EXPR_ARG,a,prec));
    1180           0 :     default: pari_err_FLAG("sumpos");
    1181             :   }
    1182             :   return NULL; /* LCOV_EXCL_LINE */
    1183             : }
    1184             : 
    1185             : /********************************************************************/
    1186             : /**                                                                **/
    1187             : /**            SEARCH FOR REAL ZEROS of an expression              **/
    1188             : /**                                                                **/
    1189             : /********************************************************************/
    1190             : /* Brent's method, [a,b] bracketing interval */
    1191             : GEN
    1192         889 : zbrent(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
    1193             : {
    1194             :   long sig, iter, itmax;
    1195         889 :   pari_sp av = avma;
    1196             :   GEN c, d, e, tol, fa, fb, fc;
    1197             : 
    1198         889 :   if (typ(a) != t_REAL || realprec(a) < prec) a = gtofp(a, prec);
    1199         889 :   if (typ(b) != t_REAL || realprec(b) < prec) b = gtofp(b, prec);
    1200         889 :   sig = cmprr(b, a);
    1201         889 :   if (!sig) return gerepileupto(av, a);
    1202         889 :   if (sig < 0) {c = a; a = b; b = c;} else c = b;
    1203         889 :   fa = eval(E, a);
    1204         889 :   fb = eval(E, b);
    1205         889 :   if (gsigne(fa)*gsigne(fb) > 0)
    1206           7 :     pari_err_DOMAIN("solve", "f(a)f(b)", ">", gen_0, mkvec2(fa, fb));
    1207         882 :   itmax = prec2nbits(prec) * 2 + 1;
    1208         882 :   tol = real2n(5-prec2nbits(prec), LOWDEFAULTPREC);
    1209         882 :   fc = fb;
    1210         882 :   e = d = NULL; /* gcc -Wall */
    1211       14036 :   for (iter = 1; iter <= itmax; ++iter)
    1212             :   {
    1213             :     GEN xm, tol1;
    1214       14036 :     if (gsigne(fb)*gsigne(fc) > 0)
    1215             :     {
    1216        7496 :       c = a; fc = fa; e = d = subrr(b, a);
    1217             :     }
    1218       14036 :     if (gcmp(gabs(fc, 0), gabs(fb, 0)) < 0)
    1219             :     {
    1220        3619 :       a = b; b = c; c = a; fa = fb; fb = fc; fc = fa;
    1221             :     }
    1222       14036 :     tol1 = abscmprr(tol, b) > 0? sqrr(tol): mulrr(tol, absr(b));
    1223       14036 :     xm = shiftr(subrr(c, b), -1);
    1224       14036 :     if (abscmprr(xm, tol1) <= 0 || gequal0(fb)) break; /* SUCCESS */
    1225             : 
    1226       13154 :     if (abscmprr(e, tol1) >= 0 && gcmp(gabs(fa, 0), gabs(fb, 0)) > 0)
    1227       11888 :     { /* attempt interpolation */
    1228       11888 :       GEN min1, min2, p, q, s = gdiv(fb, fa);
    1229       11888 :       if (cmprr(a, c) == 0)
    1230             :       {
    1231        7157 :         p = gmul2n(gmul(xm, s), 1);
    1232        7157 :         q = gsubsg(1, s);
    1233             :       }
    1234             :       else
    1235             :       {
    1236        4731 :         GEN r = gdiv(fb, fc);
    1237        4731 :         q = gdiv(fa, fc);
    1238        4731 :         p = gmul2n(gmul(gsub(q, r), gmul(xm, q)), 1);
    1239        4731 :         p = gmul(s, gsub(p, gmul(gsub(b, a), gsubgs(r, 1))));
    1240        4731 :         q = gmul(gmul(gsubgs(q, 1), gsubgs(r, 1)), gsubgs(s, 1));
    1241             :       }
    1242       11888 :       if (gsigne(p) > 0) q = gneg_i(q); else p = gneg_i(p);
    1243       11888 :       min1 = gsub(gmulsg(3, gmul(xm,q)), gabs(gmul(q, tol1), 0));
    1244       11888 :       min2 = gabs(gmul(e, q), 0);
    1245       11888 :       if (gcmp(gmul2n(p, 1), gmin(min1, min2)) < 0)
    1246       10055 :         { e = d; d = gdiv(p, q); } /* interpolation OK */
    1247             :       else
    1248        1833 :         { d = xm; e = d; } /* failed, use bisection */
    1249             :     }
    1250        1266 :     else { d = xm; e = d; } /* bound decreasing too slowly, use bisection */
    1251       13154 :     a = b; fa = fb;
    1252       13154 :     if (gcmp(gabs(d, 0), tol1) > 0) b = gadd(b, d);
    1253         838 :     else if (gsigne(xm) > 0)      b = addrr(b, tol1);
    1254         491 :     else                          b = subrr(b, tol1);
    1255       13154 :     if (realprec(b) < prec) b = rtor(b, prec);
    1256       13154 :     fb = eval(E, b);
    1257             :   }
    1258         882 :   if (iter > itmax) pari_err_IMPL("solve recovery [too many iterations]");
    1259         882 :   return gerepileuptoleaf(av, rcopy(b));
    1260             : }
    1261             : 
    1262             : GEN
    1263          21 : zbrent0(GEN a, GEN b, GEN code, long prec)
    1264          21 : { EXPR_WRAP(code, zbrent(EXPR_ARG, a, b, prec)); }
    1265             : 
    1266             : /* x = solve_start(&D, a, b, prec)
    1267             :  * while (x) {
    1268             :  *   y = ...(x);
    1269             :  *   x = solve_next(&D, y);
    1270             :  * }
    1271             :  * return D.res; */
    1272             : 
    1273             : /* Find zeros of a function in the real interval [a,b] by interval splitting */
    1274             : GEN
    1275          91 : solvestep(void *E, GEN (*f)(void *,GEN), GEN a, GEN b, GEN step, long flag, long prec)
    1276             : {
    1277          91 :   const long ITMAX = 10;
    1278          91 :   pari_sp av = avma;
    1279             :   GEN fa, ainit, binit;
    1280          91 :   long sainit, it, bit = prec2nbits(prec) / 2, ct = 0, s = gcmp(a,b);
    1281             : 
    1282          91 :   if (!s) return gequal0(f(E, a)) ? gcopy(mkvec(a)): cgetg(1,t_VEC);
    1283          91 :   if (s > 0) swap(a, b);
    1284          91 :   if (flag&4)
    1285             :   {
    1286          84 :     if (gcmpgs(step,1)<=0) pari_err_DOMAIN("solvestep","step","<=",gen_1,step);
    1287          84 :     if (gsigne(a) <= 0) pari_err_DOMAIN("solvestep","a","<=",gen_0,a);
    1288             :   }
    1289           7 :   else if (gsigne(step) <= 0)
    1290           0 :     pari_err_DOMAIN("solvestep","step","<=",gen_0,step);
    1291          91 :   ainit = a = gtofp(a, prec); fa = f(E, a);
    1292          91 :   binit = b = gtofp(b, prec); step = gtofp(step, prec);
    1293          91 :   sainit = gsigne(fa);
    1294          91 :   if (gexpo(fa) < -bit) sainit = 0;
    1295          98 :   for (it = 0; it < ITMAX; it++)
    1296             :   {
    1297          98 :     pari_sp av2 = avma;
    1298          98 :     GEN v = cgetg(1, t_VEC);
    1299             :     long sa;
    1300          98 :     a = ainit;
    1301          98 :     b = binit;
    1302          98 :     sa = sainit;
    1303        2457 :     while (gcmp(a,b) < 0)
    1304             :     {
    1305        2261 :       GEN fc, c = (flag&4)? gmul(a, step): gadd(a, step);
    1306             :       long sc;
    1307        2261 :       if (gcmp(c,b) > 0) c = b;
    1308        2261 :       fc = f(E, c);
    1309        2261 :       sc = gsigne(fc);
    1310        2261 :       if (gexpo(fc) < -bit) sc = 0;
    1311        2261 :       if (!sc || sa*sc < 0)
    1312             :       {
    1313             :         long e;
    1314             :         GEN z;
    1315         441 :         z = sc? zbrent(E, f, a, c, prec): c;
    1316         441 :         (void)grndtoi(z, &e);
    1317         441 :         if (e  <= -bit) ct = 1;
    1318         441 :         if ((flag&1) && ((!(flag&8)) || ct)) return gerepileupto(av, z);
    1319         441 :         v = gconcat(v, z);
    1320             :       }
    1321        2261 :       a = c; fa = fc; sa = sc;
    1322             :     }
    1323          98 :     if ((!(flag&2) || lg(v) > 1) && (!(flag&8) || ct))
    1324          91 :       return gerepilecopy(av, v);
    1325           7 :     step = (flag&4)? sqrtr(sqrtr(step)): gmul2n(step, -2);
    1326           7 :     gerepileall(av2, 2, &fa, &step);
    1327             :   }
    1328           0 :   if (it == ITMAX) pari_err_IMPL("solvestep recovery [too many iterations]");
    1329           0 :   return NULL;
    1330             : }
    1331             : 
    1332             : GEN
    1333           7 : solvestep0(GEN a, GEN b, GEN step, GEN code, long flag, long prec)
    1334           7 : { EXPR_WRAP(code, solvestep(EXPR_ARG, a,b, step, flag, prec)); }
    1335             : 
    1336             : /********************************************************************/
    1337             : /**                     Numerical derivation                       **/
    1338             : /********************************************************************/
    1339             : 
    1340             : struct deriv_data
    1341             : {
    1342             :   GEN code;
    1343             :   GEN args;
    1344             : };
    1345             : 
    1346         105 : static GEN deriv_eval(void *E, GEN x, long prec)
    1347             : {
    1348         105 :  struct deriv_data *data=(struct deriv_data *)E;
    1349         105 :  gel(data->args,1)=x;
    1350         105 :  return closure_callgenvecprec(data->code, data->args, prec);
    1351             : }
    1352             : 
    1353             : /* Rationale: (f(2^-e) - f(-2^-e) + O(2^-pr)) / (2 * 2^-e) = f'(0) + O(2^-2e)
    1354             :  * since 2nd derivatives cancel.
    1355             :  *   prec(LHS) = pr - e
    1356             :  *   prec(RHS) = 2e, equal when  pr = 3e = 3/2 fpr (fpr = required final prec)
    1357             :  *
    1358             :  * For f'(x), x far from 0: prec(LHS) = pr - e - expo(x)
    1359             :  * --> pr = 3/2 fpr + expo(x) */
    1360             : GEN
    1361        1008 : derivnum(void *E, GEN (*eval)(void *, GEN, long), GEN x, long prec)
    1362             : {
    1363             :   GEN eps,a,b, y;
    1364             :   long pr, l, e, ex, newprec;
    1365        1008 :   pari_sp av = avma;
    1366        1008 :   long p = precision(x);
    1367        1008 :   long fpr = p ? prec2nbits(p): prec2nbits(prec);
    1368        1008 :   ex = gexpo(x);
    1369        1008 :   if (ex < 0) ex = 0; /* near 0 */
    1370        1008 :   pr = (long)ceil(fpr * 1.5 + ex);
    1371        1008 :   l = nbits2prec(pr);
    1372        1008 :   newprec = nbits2prec(pr + ex + BITS_IN_LONG);
    1373        1008 :   switch(typ(x))
    1374             :   {
    1375             :     case t_REAL:
    1376             :     case t_COMPLEX:
    1377         406 :       x = gprec_w(x, newprec);
    1378             :   }
    1379             : 
    1380        1008 :   e = fpr/2; /* 1/2 required prec (in sig. bits) */
    1381        1008 :   eps = real2n(-e, l);
    1382        1008 :   a = eval(E, gsub(x, eps), newprec);
    1383        1008 :   b = eval(E, gadd(x, eps), newprec);
    1384        1008 :   y = gmul2n(gsub(b,a), e-1);
    1385        1008 :   return gerepileupto(av, gprec_w(y, nbits2prec(fpr)));
    1386             : }
    1387             : 
    1388             : /* Fornberg interpolation algorithm for finite differences coefficients
    1389             : * using N+1 equidistant grid points around 0 [ assume N even >= M ].
    1390             : * Compute \delta[m]_{N,nu} for all derivation orders m = 0..M such that
    1391             : *   h^m * f^{(m)}(0) = \sum_{nu = 0}^n delta[m]_{N,nu}  f(a_nu) + O(h^{N-m+1}),
    1392             : * for step size h.
    1393             : * Return a = [0,-1,1...,-N2,N2] and vector of vectors d: d[m+1][nu+1]
    1394             : * = (M!/m!) * delta[m]_{N,nu}, nu = 0..N */
    1395             : static void
    1396          28 : FD(long M, long N, GEN *pd, GEN *pa)
    1397             : {
    1398             :   GEN d, a, b, W, Wp, t, F, Mfact;
    1399             :   long N2, m, nu, i;
    1400             : 
    1401          28 :   if (odd(N)) N++; /* make it even */
    1402          28 :   N2 = N>>1;
    1403          28 :   F = cgetg(N+2, t_VEC);
    1404          28 :   a = cgetg(N+2, t_VEC);
    1405          28 :   b = cgetg(N2+1, t_VEC);
    1406          28 :   gel(a,1) = gen_0;
    1407         357 :   for (i = 1; i <= N2; i++)
    1408             :   {
    1409         329 :     gel(a,2*i)   = utoineg(i);
    1410         329 :     gel(a,2*i+1) = utoipos(i);
    1411         329 :     gel(b,i) = sqru(i);
    1412             :   }
    1413             :   /* w = \prod (X - a[i]) = x W(x^2) */
    1414          28 :   Mfact = mpfact(M);
    1415          28 :   W = roots_to_pol(b, 0);
    1416          28 :   Wp = ZX_deriv(W);
    1417          28 :   t = gel(W,2); /* w'(0) */
    1418          28 :   t = diviiexact(t, Mfact);
    1419          28 :   gel(F,1) = RgX_Rg_div(RgX_inflate(W,2), t);
    1420         357 :   for (i = 1; i <= N2; i++)
    1421             :   { /* t = w'(a_{2i}) = w'(a_{2i+1}) */
    1422         329 :     GEN r, t = mulii(shifti(gel(b,i),1), poleval(Wp, gel(b,i)));
    1423             :     GEN U, S, T;
    1424         329 :     U = RgX_inflate(RgX_div_by_X_x(W, gel(b,i), &r), 2);
    1425         329 :     U = RgX_shift_shallow(U, 1);
    1426         329 :     U = RgXn_red_shallow(U, M+1); /* higher terms not needed */
    1427         329 :     t = diviiexact(t, Mfact);
    1428         329 :     U = RgX_Rg_div(U, t);
    1429         329 :     S = RgX_shift_shallow(U,1);
    1430         329 :     T = RgX_Rg_mul(U, gel(a,2*i+1));
    1431         329 :     gel(F,2*i)   = RgX_sub(S, T);
    1432         329 :     gel(F,2*i+1) = RgX_add(S, T);
    1433             :   }
    1434             :   /* F[i] = M! w(X) / ((X-a[i])w'(a[i])) + O(X^(M+1)) */
    1435          28 :   d = cgetg(M+2, t_VEC);
    1436         280 :   for (m = 0; m <= M; m++)
    1437             :   {
    1438         252 :     GEN v = cgetg(N+2, t_VEC); /* coeff(F[nu],X^m) */
    1439         252 :     for (nu = 0; nu <= N; nu++) gel(v, nu+1) = gmael(F, nu+1, m+2);
    1440         252 :     gel(d,m+1) = v;
    1441             :   }
    1442          28 :   *pd = d;
    1443          28 :   *pa = a;
    1444          28 : }
    1445             : 
    1446             : static void
    1447         231 : chk_ord(long m)
    1448             : {
    1449         231 :   if (m < 0)
    1450          14 :     pari_err_DOMAIN("derivnumk", "derivation order", "<", gen_0, stoi(m));
    1451         217 : }
    1452             : 
    1453             : GEN
    1454          42 : derivnumk(void *E, GEN (*eval)(void *, GEN, long), GEN x, GEN ind0, long prec)
    1455             : {
    1456             :   GEN A, D, X, F, ind;
    1457             :   long M, fpr, p, i, pr, l, lA, e, ex, eD, newprec;
    1458          42 :   pari_sp av = avma;
    1459          42 :   int allodd = 1;
    1460             : 
    1461          42 :   ind = gtovecsmall(ind0);
    1462          42 :   l = lg(ind);
    1463          42 :   F = cgetg(l, t_VEC);
    1464          42 :   M = vecsmall_max(ind);
    1465          42 :   chk_ord(M);
    1466          42 :   if (!M) /* silly degenerate case */
    1467             :   {
    1468          14 :     X = eval(E, x, prec);
    1469          14 :     for (i = 1; i < l; i++) { chk_ord(ind[i]); gel(F,i) = X; }
    1470           7 :     if (typ(ind0) == t_INT) F = gel(F,1);
    1471           7 :     return gerepilecopy(av, F);
    1472             :   }
    1473          28 :   FD(M, 3*M-1, &D,&A); /* optimal if 'eval' uses quadratic time */
    1474             : 
    1475          28 :   p = precision(x);
    1476          28 :   fpr = p ? prec2nbits(p): prec2nbits(prec);
    1477          28 :   eD = gexpo(gel(D,M));
    1478          28 :   e = (fpr + 3*M*log2((double)M)) / (2*M);
    1479          28 :   ex = gexpo(x);
    1480          28 :   if (ex < 0) ex = 0; /* near 0 */
    1481          28 :   pr = (long)ceil(fpr + e * M); /* ~ 3fpr/2 */
    1482          28 :   newprec = nbits2prec(pr + eD + ex + BITS_IN_LONG);
    1483          28 :   switch(typ(x))
    1484             :   {
    1485             :     case t_REAL:
    1486             :     case t_COMPLEX:
    1487           0 :       x = gprec_w(x, newprec);
    1488             :   }
    1489          28 :   lA = lg(A); X = cgetg(lA, t_VEC);
    1490          42 :   for (i = 1; i < l; i++)
    1491          35 :     if (!odd(ind[i])) { allodd = 0; break; }
    1492             :   /* if only odd derivation orders, the value at 0 (A[1]) is not needed */
    1493          28 :   gel(X, 1) = gen_0;
    1494         707 :   for (i = allodd? 2: 1; i < lA; i++)
    1495             :   {
    1496         679 :     GEN t = eval(E, gadd(x, gmul2n(gel(A,i), -e)), newprec);
    1497         679 :     if (!gprecision(t)) t = gtofp(t, newprec);
    1498         679 :     gel(X, i) = t;
    1499             :   }
    1500             : 
    1501         147 :   for (i = 1; i < l; i++)
    1502             :   {
    1503             :     GEN t;
    1504         126 :     long m = ind[i]; chk_ord(m);
    1505         119 :     t = gmul2n(RgV_dotproduct(gel(D,m+1), X), e*m);
    1506         119 :     if (m < M) t = gdiv(t, mulu_interval(m+1,M));
    1507         119 :     gel(F,i) = t;
    1508             :   }
    1509          21 :   if (typ(ind0) == t_INT) F = gel(F,1);
    1510          21 :   return gerepileupto(av, gprec_w(F, nbits2prec(fpr)));
    1511             : }
    1512             : /* v(t') */
    1513             : static long
    1514          14 : rfrac_val_deriv(GEN t)
    1515             : {
    1516          14 :   long v = varn(gel(t,2));
    1517          14 :   return gvaluation(deriv(t, v), pol_x(v));
    1518             : }
    1519             : 
    1520             : GEN
    1521        1043 : derivfunk(void *E, GEN (*eval)(void *, GEN, long), GEN x, GEN ind0, long prec)
    1522             : {
    1523             :   pari_sp av;
    1524             :   GEN ind, xp, ixp, F, G;
    1525             :   long i, l, vx, M;
    1526        1043 :   if (!ind0) return derivfun(E, eval, x, prec);
    1527          63 :   switch(typ(x))
    1528             :   {
    1529             :   case t_REAL: case t_INT: case t_FRAC: case t_COMPLEX:
    1530          42 :     return derivnumk(E,eval, x, ind0, prec);
    1531             :   case t_POL:
    1532          14 :     ind = gtovecsmall(ind0);
    1533          14 :     M = vecsmall_max(ind);
    1534          14 :     xp = RgX_deriv(x);
    1535          14 :     x = RgX_to_ser(x, precdl+2 + M * (1+RgX_val(xp)));
    1536          14 :     break;
    1537             :   case t_RFRAC:
    1538           7 :     ind = gtovecsmall(ind0);
    1539           7 :     M = vecsmall_max(ind);
    1540           7 :     x = rfrac_to_ser(x, precdl+2 + M * (1+rfrac_val_deriv(x)));
    1541           7 :     xp = derivser(x);
    1542           7 :     break;
    1543             :   case t_SER:
    1544           0 :     ind = gtovecsmall(ind0);
    1545           0 :     M = vecsmall_max(ind);
    1546           0 :     xp = derivser(x);
    1547           0 :     break;
    1548           0 :   default: pari_err_TYPE("numerical derivation",x);
    1549             :     return NULL; /*LCOV_EXCL_LINE*/
    1550             :   }
    1551          21 :   av = avma; chk_ord(M);
    1552          21 :   vx = varn(x);
    1553          21 :   ixp = M? ginv(xp): NULL;
    1554          21 :   F = cgetg(M+2, t_VEC);
    1555          21 :   gel(F,1) = eval(E, x, prec);
    1556          21 :   for (i = 1; i <= M; i++) gel(F,i+1) = gmul(deriv(gel(F,i),vx), ixp);
    1557          21 :   l = lg(ind); G = cgetg(l, t_VEC);
    1558          42 :   for (i = 1; i < l; i++)
    1559             :   {
    1560          21 :     long m = ind[i]; chk_ord(m);
    1561          21 :     gel(G,i) = gel(F,m+1);
    1562             :   }
    1563          21 :   if (typ(ind0) == t_INT) G = gel(G,1);
    1564          21 :   return gerepilecopy(av, G);
    1565             : }
    1566             : 
    1567             : GEN
    1568        1043 : derivfun(void *E, GEN (*eval)(void *, GEN, long), GEN x, long prec)
    1569             : {
    1570        1043 :   pari_sp av = avma;
    1571             :   GEN xp;
    1572             :   long vx;
    1573        1043 :   switch(typ(x))
    1574             :   {
    1575             :   case t_REAL: case t_INT: case t_FRAC: case t_COMPLEX:
    1576        1008 :     return derivnum(E,eval, x, prec);
    1577             :   case t_POL:
    1578          14 :     xp = RgX_deriv(x);
    1579          14 :     x = RgX_to_ser(x, precdl+2+ (1 + RgX_val(xp)));
    1580          14 :     break;
    1581             :   case t_RFRAC:
    1582           7 :     x = rfrac_to_ser(x, precdl+2+ (1 + rfrac_val_deriv(x)));
    1583             :     /* fall through */
    1584             :   case t_SER:
    1585          14 :     xp = derivser(x);
    1586          14 :     break;
    1587           7 :   default: pari_err_TYPE("formal derivation",x);
    1588             :     return NULL; /*LCOV_EXCL_LINE*/
    1589             :   }
    1590          28 :   vx = varn(x);
    1591          28 :   return gerepileupto(av, gdiv(deriv(eval(E, x, prec),vx), xp));
    1592             : }
    1593             : 
    1594             : GEN
    1595        1043 : derivnum0(GEN a, GEN code, GEN ind, long prec)
    1596        1043 : { EXPR_WRAP(code, derivfunk(EXPR_ARGPREC,a,ind,prec)); }
    1597             : 
    1598             : GEN
    1599          63 : derivfun0(GEN code, GEN args, long prec)
    1600             : {
    1601             :   struct deriv_data E;
    1602          63 :   E.code=code; E.args=args;
    1603          63 :   return derivfun((void*)&E, deriv_eval, gel(args,1), prec);
    1604             : }
    1605             : 
    1606             : /********************************************************************/
    1607             : /**                   Numerical extrapolation                      **/
    1608             : /********************************************************************/
    1609             : /* for t_CLOSURE */
    1610             : static double
    1611          77 : fun_getmf(long mul)
    1612             : {
    1613          77 :   const double A[] = {0,
    1614             :     0.331,0.260,0.228,0.210,0.198,
    1615             :     0.188,0.181,0.175,0.170,0.166 };
    1616          77 :   const double B[] = {0,
    1617             :     0.166,0.142,0.132,0.125,0.120,
    1618             :     0.116,0.113,0.111,0.109,0.107 };
    1619          77 :   if (mul <=  10) return A[mul];
    1620          70 :   if (mul <= 109) return B[mul/10]; else return 0.105;
    1621             : }
    1622             : /* for t_VEC */
    1623             : static double
    1624          14 : vec_getmf(long mul)
    1625             : {
    1626          14 :   const double A[] = {0,
    1627             :     0.2062, 0.1760, 0.1613, 0.1519, 0.1459,
    1628             :     0.1401, 0.1360, 0.1327, 0.1299, 0.1280 };
    1629          14 :   const double B[] = {0,
    1630             :     0.1280, 0.1133, 0.1064, 0.1019, 0.0987,
    1631             :     0.0962, 0.0942, 0.0925, 0.0911, 0.0899 };
    1632          14 :   if (mul <=  10) return A[mul];
    1633          14 :   if (mul <= 109) return B[mul/10]; else return 0.0899;
    1634             : }
    1635             : 
    1636             : /* [u(n*muli), u <= N], muli = 1 unless f!=NULL */
    1637             : static GEN
    1638          91 : get_u(void *E, GEN (*f)(void *, GEN, long), long N, long muli, long prec)
    1639             : {
    1640             :   long n;
    1641          91 :   GEN u = cgetg(N+1, t_VEC);
    1642          91 :   if (f)
    1643             :   {
    1644          77 :     for (n = 1; n <= N; n++) gel(u,n) = f(E, stoi(muli*n), prec);
    1645             :   }
    1646             :   else
    1647             :   {
    1648          14 :     GEN v = (GEN)E;
    1649          14 :     long t = lg(v)-1;
    1650          14 :     if (t < N*muli) pari_err_COMPONENT("limitnum","<",stoi(N), stoi(t));
    1651          14 :     for (n = 1; n <= N; n++) gel(u,n) = gel(v, muli*n);
    1652             :   }
    1653        3997 :   for (n = 1; n <= N; n++)
    1654             :   {
    1655        3906 :     GEN un = gel(u,n);
    1656        3906 :     if (is_rational_t(typ(un))) gel(u,n) = gtofp(un, prec);
    1657             :   }
    1658          91 :   return u;
    1659             : }
    1660             : 
    1661             : struct limit
    1662             : {
    1663             :   long prec0; /* target accuracy */
    1664             :   long prec; /* working accuracy */
    1665             :   long N; /* number of terms */
    1666             :   GEN u; /* sequence to extrapolate */
    1667             :   GEN na; /* [n^alpha, n <= N] */
    1668             :   GEN nma; /* [n^-alpha, n <= N] or NULL (alpha = 1) */
    1669             :   GEN coef; /* or NULL (alpha != 1) */
    1670             : };
    1671             : 
    1672             : static void
    1673          91 : limit_init(struct limit *L, void *E, GEN (*f)(void*,GEN,long),
    1674             :            long muli, GEN alpha, long prec)
    1675             : {
    1676          91 :   long bitprec = prec2nbits(prec), n, N;
    1677             :   GEN na;
    1678             : 
    1679          91 :   if (muli <= 0) muli = 20;
    1680          91 :   L->N = N = (long)ceil((f? fun_getmf(muli): vec_getmf(muli)) * bitprec);
    1681          91 :   L->prec = nbits2prec(bitprec + (long)ceil(1.844*N));
    1682          91 :   L->prec0 = prec;
    1683          91 :   L->u = get_u(E, f, N, muli, L->prec);
    1684          91 :   if (alpha && !gequal1(alpha))
    1685          28 :   {
    1686          28 :     long prec2 = gprecision(alpha);
    1687             :     GEN nma;
    1688          28 :     if (!prec2) prec2 = L->prec;
    1689          28 :     na = vecpowug(N, alpha, prec2);
    1690          28 :     L->coef = NULL;
    1691          28 :     L->nma = nma = cgetg(N+1, t_VEC);
    1692          28 :     for (n = 1; n <= N; n++) gel(nma, n) = ginv(gel(na, n));
    1693          28 :     if (muli != 1) na = gmul(na, gpow(utor(muli,prec2), alpha, prec2));
    1694             :   }
    1695             :   else
    1696             :   {
    1697          63 :     GEN coef, C = vecbinomial(N), T = vecpowuu(N, N);
    1698          63 :     na = cgetg(N+1, t_VEC);
    1699          63 :     L->coef = coef = cgetg(N+1, t_VEC);
    1700          63 :     L->nma = NULL;
    1701        2366 :     for (n = 1; n <= N; n++)
    1702             :     {
    1703        2303 :       GEN c = mulii(gel(C,n+1), gel(T,n));
    1704        2303 :       if (odd(N-n)) togglesign_safe(&c);
    1705        2303 :       gel(coef, n) = c;
    1706        2303 :       gel(na, n) = utoipos(n*muli);
    1707             :     }
    1708             :   }
    1709          91 :   L->na = na;
    1710          91 : }
    1711             : 
    1712             : /* Zagier/Lagrange extrapolation */
    1713             : static GEN
    1714         777 : limitnum_i(struct limit *L)
    1715             : {
    1716         777 :   pari_sp av = avma;
    1717             :   GEN S;
    1718         777 :   if (L->nma)
    1719         343 :     S = polint(L->nma, L->u,gen_0,NULL);
    1720             :   else
    1721         434 :     S = gdiv(RgV_dotproduct(L->u,L->coef), mpfact(L->N));
    1722         777 :   return gerepilecopy(av, gprec_w(S, L->prec0));
    1723             : }
    1724             : GEN
    1725          49 : limitnum(void *E, GEN (*f)(void *, GEN, long), long muli, GEN alpha, long prec)
    1726             : {
    1727             :   struct limit L;
    1728          49 :   limit_init(&L, E,f, muli, alpha, prec);
    1729          49 :   return limitnum_i(&L);
    1730             : }
    1731             : GEN
    1732          49 : limitnum0(GEN u, long muli, GEN alpha, long prec)
    1733             : {
    1734          49 :   void *E = (void*)u;
    1735          49 :   GEN (*f)(void*,GEN,long) = NULL;
    1736          49 :   switch(typ(u))
    1737             :   {
    1738             :     case t_COL:
    1739           7 :     case t_VEC: break;
    1740          42 :     case t_CLOSURE: f = gp_callprec; break;
    1741           0 :     default: pari_err_TYPE("limitnum", u);
    1742             :   }
    1743          49 :   return limitnum(E,f, muli,alpha, prec);
    1744             : }
    1745             : 
    1746             : GEN
    1747          42 : asympnum(void *E, GEN (*f)(void *, GEN, long), long muli, GEN alpha, long prec)
    1748             : {
    1749          42 :   const long MAX = 100;
    1750          42 :   pari_sp av = avma;
    1751          42 :   GEN u, vres = vectrunc_init(MAX);
    1752          42 :   long i, B = prec2nbits(prec);
    1753          42 :   double LB = 0.9*expu(B); /* 0.9 and 0.95 below are heuristic */
    1754             :   struct limit L;
    1755          42 :   limit_init(&L, E,f, muli, alpha, prec);
    1756          42 :   if (alpha) LB *= gtodouble(alpha);
    1757          42 :   u = L.u;
    1758         728 :   for(i = 1; i <= MAX; i++)
    1759             :   {
    1760             :     GEN a, s, v, p, q;
    1761             :     long n;
    1762         728 :     s = limitnum_i(&L);
    1763             :     /* NOT bestappr: lindep properly ignores the lower bits */
    1764         728 :     v = lindep_bit(mkvec2(gen_1, s), maxss((long)(0.95*floor(B - i*LB)), 32));
    1765         728 :     p = negi(gel(v,1));
    1766         728 :     q = gel(v,2);
    1767         728 :     if (!signe(q)) break;
    1768         728 :     a = gdiv(p,q);
    1769         728 :     s = gsub(s, a);
    1770             :     /* |s|q^2 > eps */
    1771         728 :     if (!gequal0(s) && gexpo(s) + 2*expi(q) > -17) break;
    1772         686 :     vectrunc_append(vres, a);
    1773         686 :     for (n = 1; n <= L.N; n++) gel(u,n) = gmul(gsub(gel(u,n), a), gel(L.na,n));
    1774             :   }
    1775          42 :   return gerepilecopy(av, vres);
    1776             : }
    1777             : GEN
    1778          42 : asympnum0(GEN u, long muli, GEN alpha, long prec)
    1779             : {
    1780          42 :   void *E = (void*)u;
    1781          42 :   GEN (*f)(void*,GEN,long) = NULL;
    1782          42 :   switch(typ(u))
    1783             :   {
    1784             :     case t_COL:
    1785           7 :     case t_VEC: break;
    1786          35 :     case t_CLOSURE: f = gp_callprec; break;
    1787           0 :     default: pari_err_TYPE("asympnum", u);
    1788             :   }
    1789          42 :   return asympnum(E,f, muli,alpha, prec);
    1790             : }

Generated by: LCOV version 1.11