Code coverage tests

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

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

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

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

Generated by: LCOV version 1.14