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

Generated by: LCOV version 1.11