Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - language - sumiter.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19047-48ce0fc) Lines: 776 809 95.9 %
Date: 2016-06-28 Functions: 65 65 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 438 565 77.5 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : #include "pari.h"
      15                 :            : #include "paripriv.h"
      16                 :            : #include "anal.h"
      17                 :            : 
      18                 :            : GEN
      19                 :    1235987 : iferrpari(GEN a, GEN b, GEN c)
      20                 :            : {
      21                 :            :   GEN res;
      22                 :            :   struct pari_evalstate state;
      23                 :    1235987 :   evalstate_save(&state);
      24         [ +  + ]:    1235987 :   pari_CATCH(CATCH_ALL)
      25                 :            :   {
      26                 :            :     GEN E;
      27 [ -  + ][ #  # ]:      15012 :     if (!b&&!c) return gnil;
      28                 :      15012 :     E = evalstate_restore_err(&state);
      29         [ +  + ]:      15012 :     if (c)
      30                 :            :     {
      31                 :        144 :       push_lex(E,c);
      32                 :        144 :       res = closure_evalgen(c);
      33                 :        144 :       pop_lex(1);
      34         [ +  + ]:        144 :       if (gequal0(res))
      35                 :          7 :         pari_err(0, E);
      36                 :            :     }
      37         [ -  + ]:      15005 :     if (!b) return gnil;
      38                 :      15005 :     push_lex(E,b);
      39                 :      15005 :     res = closure_evalgen(b);
      40                 :      15005 :     pop_lex(1);
      41                 :      15005 :     return res;
      42                 :            :   } pari_TRY {
      43                 :    1235987 :     res = closure_evalgen(a);
      44                 :    1220976 :   } pari_ENDCATCH;
      45                 :    1235981 :   return res;
      46                 :            : }
      47                 :            : 
      48                 :            : /********************************************************************/
      49                 :            : /**                                                                **/
      50                 :            : /**                        ITERATIONS                              **/
      51                 :            : /**                                                                **/
      52                 :            : /********************************************************************/
      53                 :            : 
      54                 :            : static void
      55                 :    3799618 : forparii(GEN a, GEN b, GEN code)
      56                 :            : {
      57                 :    3799618 :   pari_sp av, av0 = avma;
      58                 :            :   GEN aa;
      59         [ +  + ]:    7544944 :   if (gcmp(b,a) < 0) return;
      60         [ +  + ]:    3745383 :   if (typ(b) != t_INFINITY) b = gfloor(b);
      61                 :    3745376 :   aa = a = setloop(a);
      62                 :    3745381 :   av=avma;
      63                 :    3745381 :   push_lex(a,code);
      64         [ +  + ]:   46932359 :   while (gcmp(a,b) <= 0)
      65                 :            :   {
      66         [ +  + ]:   43466317 :     closure_evalvoid(code); if (loop_break()) break;
      67                 :   43248770 :     a = get_lex(-1);
      68         [ +  + ]:   43246903 :     if (a == aa)
      69                 :            :     {
      70                 :   43246875 :       a = incloop(a);
      71         [ +  + ]:   43186944 :       if (a != aa) { set_lex(-1,a); aa = a; }
      72                 :            :     }
      73                 :            :     else
      74                 :            :     { /* 'code' modified a ! Be careful (and slow) from now on */
      75                 :         28 :       a = gaddgs(a,1);
      76         [ -  + ]:         28 :       if (gc_needed(av,1))
      77                 :            :       {
      78         [ #  # ]:          0 :         if (DEBUGMEM>1) pari_warn(warnmem,"forparii");
      79                 :          0 :         a = gerepileupto(av,a);
      80                 :            :       }
      81                 :         28 :       set_lex(-1,a);
      82                 :            :     }
      83                 :            :   }
      84                 :    3679653 :   pop_lex(1);  avma = av0;
      85                 :            : }
      86                 :            : 
      87                 :            : void
      88                 :    3799625 : forpari(GEN a, GEN b, GEN code)
      89                 :            : {
      90                 :    3799625 :   pari_sp ltop=avma, av;
      91         [ +  + ]:    3799632 :   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); avma = ltop;
     107                 :            : }
     108                 :            : 
     109                 :            : void
     110                 :    1169651 : whilepari(GEN a, GEN b)
     111                 :            : {
     112                 :    1169651 :   pari_sp av = avma;
     113                 :            :   for(;;)
     114                 :            :   {
     115                 :   16126236 :     GEN res = closure_evalnobrk(a);
     116         [ +  + ]:   16126236 :     if (gequal0(res)) break;
     117                 :   14956585 :     avma = av;
     118         [ -  + ]:   14956585 :     closure_evalvoid(b); if (loop_break()) break;
     119                 :   14956585 :   }
     120                 :    1169651 :   avma = av;
     121                 :    1169651 : }
     122                 :            : 
     123                 :            : void
     124                 :     222070 : untilpari(GEN a, GEN b)
     125                 :            : {
     126                 :     222070 :   pari_sp av = avma;
     127                 :            :   for(;;)
     128                 :            :   {
     129                 :            :     GEN res;
     130         [ -  + ]:    2194053 :     closure_evalvoid(b); if (loop_break()) break;
     131                 :    2194053 :     res = closure_evalnobrk(a);
     132         [ +  + ]:    2194053 :     if (!gequal0(res)) break;
     133                 :    1971983 :     avma = av;
     134                 :    1971983 :   }
     135                 :     222070 :   avma = av;
     136                 :     222070 : }
     137                 :            : 
     138                 :         28 : static int negcmp(GEN x, GEN y) { return gcmp(y,x); }
     139                 :            : 
     140                 :            : void
     141                 :       1484 : forstep(GEN a, GEN b, GEN s, GEN code)
     142                 :            : {
     143                 :            :   long ss, i;
     144                 :       1484 :   pari_sp av, av0 = avma;
     145                 :       1484 :   GEN v = NULL;
     146                 :            :   int (*cmp)(GEN,GEN);
     147                 :            : 
     148                 :       1484 :   b = gcopy(b); s = gcopy(s); av=avma;
     149                 :       1484 :   push_lex(a,code);
     150         [ +  + ]:       1484 :   if (is_vec_t(typ(s)))
     151                 :            :   {
     152                 :          7 :     v = s; s = gen_0;
     153         [ +  + ]:         21 :     for (i=lg(v)-1; i; i--) s = gadd(s,gel(v,i));
     154                 :            :   }
     155                 :       1484 :   ss = gsigne(s);
     156         [ +  + ]:       1484 :   if (!ss) pari_err_DOMAIN("forstep","step","=",gen_0,s);
     157         [ +  + ]:       1477 :   cmp = (ss > 0)? &gcmp: &negcmp;
     158                 :       1477 :   i = 0;
     159         [ +  + ]:      49007 :   while (cmp(a,b) <= 0)
     160                 :            :   {
     161         [ -  + ]:      47530 :     closure_evalvoid(code); if (loop_break()) break;
     162         [ +  + ]:      47530 :     if (v)
     163                 :            :     {
     164         [ +  + ]:         42 :       if (++i >= lg(v)) i = 1;
     165                 :         42 :       s = gel(v,i);
     166                 :            :     }
     167                 :      47530 :     a = get_lex(-1); a = gadd(a,s);
     168                 :            : 
     169         [ -  + ]:      47530 :     if (gc_needed(av,1))
     170                 :            :     {
     171         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"forstep");
     172                 :          0 :       a = gerepileupto(av,a);
     173                 :            :     }
     174                 :      47530 :     set_lex(-1,a);
     175                 :            :   }
     176                 :       1477 :   pop_lex(1); avma = av0;
     177                 :       1477 : }
     178                 :            : 
     179                 :            : void
     180                 :          7 : fordiv(GEN a, GEN code)
     181                 :            : {
     182                 :            :   long i, l;
     183                 :          7 :   pari_sp av2, av = avma;
     184                 :          7 :   GEN t = divisors(a);
     185                 :            : 
     186                 :          7 :   push_lex(gen_0,code); l=lg(t); av2 = avma;
     187         [ +  + ]:         35 :   for (i=1; i<l; i++)
     188                 :            :   {
     189                 :         28 :     set_lex(-1,gel(t,i));
     190         [ -  + ]:         28 :     closure_evalvoid(code); if (loop_break()) break;
     191                 :         28 :     avma = av2;
     192                 :            :   }
     193                 :          7 :   pop_lex(1); avma=av;
     194                 :          7 : }
     195                 :            : 
     196                 :            : /* Embedded for loops:
     197                 :            :  *   fl = 0: execute ch (a), where a = (ai) runs through all n-uplets in
     198                 :            :  *     [m1,M1] x ... x [mn,Mn]
     199                 :            :  *   fl = 1: impose a1 <= ... <= an
     200                 :            :  *   fl = 2:        a1 <  ... <  an
     201                 :            :  */
     202                 :            : /* increment and return d->a [over integers]*/
     203                 :            : static GEN
     204                 :     130584 : _next_i(forvec_t *d)
     205                 :            : {
     206                 :     130584 :   long i = d->n;
     207         [ +  + ]:     130584 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     208                 :            :   for (;;) {
     209         [ +  + ]:     157867 :     if (cmpii(d->a[i], d->M[i]) < 0) {
     210                 :     130342 :       d->a[i] = incloop(d->a[i]);
     211                 :     130342 :       return (GEN)d->a;
     212                 :            :     }
     213                 :      27525 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     214         [ +  + ]:      27525 :     if (--i <= 0) return NULL;
     215                 :     157995 :   }
     216                 :            : }
     217                 :            : /* increment and return d->a [generic]*/
     218                 :            : static GEN
     219                 :         63 : _next(forvec_t *d)
     220                 :            : {
     221                 :         63 :   long i = d->n;
     222         [ +  + ]:         63 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     223                 :            :   for (;;) {
     224                 :         98 :     d->a[i] = gaddgs(d->a[i], 1);
     225         [ +  + ]:         98 :     if (gcmp(d->a[i], d->M[i]) <= 0) return (GEN)d->a;
     226                 :         49 :     d->a[i] = d->m[i];
     227         [ +  + ]:         49 :     if (--i <= 0) return NULL;
     228                 :        105 :   }
     229                 :            : }
     230                 :            : 
     231                 :            : /* non-decreasing order [over integers] */
     232                 :            : static GEN
     233                 :        113 : _next_le_i(forvec_t *d)
     234                 :            : {
     235                 :        113 :   long i = d->n;
     236         [ +  + ]:        113 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     237                 :            :   for (;;) {
     238         [ +  + ]:        175 :     if (cmpii(d->a[i], d->M[i]) < 0)
     239                 :            :     {
     240                 :         81 :       d->a[i] = incloop(d->a[i]);
     241                 :            :       /* m[i] < a[i] <= M[i] <= M[i+1] */
     242         [ +  + ]:        136 :       while (i < d->n)
     243                 :            :       {
     244                 :            :         GEN t;
     245                 :         55 :         i++;
     246         [ -  + ]:         55 :         if (cmpii(d->a[i-1], d->a[i]) <= 0) continue;
     247                 :            :         /* a[i] < a[i-1] <= M[i-1] <= M[i] */
     248         [ -  + ]:         55 :         t = d->a[i-1]; if (cmpii(t, d->m[i]) < 0) t = d->m[i];
     249                 :         55 :         d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1],m[i])*/
     250                 :            :       }
     251                 :         81 :       return (GEN)d->a;
     252                 :            :     }
     253                 :         94 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     254         [ +  + ]:         94 :     if (--i <= 0) return NULL;
     255                 :        191 :   }
     256                 :            : }
     257                 :            : /* non-decreasing order [generic] */
     258                 :            : static GEN
     259                 :        154 : _next_le(forvec_t *d)
     260                 :            : {
     261                 :        154 :   long i = d->n;
     262         [ +  + ]:        154 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     263                 :            :   for (;;) {
     264                 :        266 :     d->a[i] = gaddgs(d->a[i], 1);
     265         [ +  + ]:        266 :     if (gcmp(d->a[i], d->M[i]) <= 0)
     266                 :            :     {
     267         [ +  + ]:        224 :       while (i < d->n)
     268                 :            :       {
     269                 :            :         GEN c;
     270                 :         98 :         i++;
     271         [ -  + ]:         98 :         if (gcmp(d->a[i-1], d->a[i]) <= 0) continue;
     272                 :            :         /* M[i] >= M[i-1] >= a[i-1] > a[i] */
     273                 :         98 :         c = gceil(gsub(d->a[i-1], d->a[i]));
     274                 :         98 :         d->a[i] = gadd(d->a[i], c);
     275                 :            :         /* a[i-1] <= a[i] < M[i-1] + 1 => a[i] < M[i]+1 => a[i] <= M[i] */
     276                 :            :       }
     277                 :        126 :       return (GEN)d->a;
     278                 :            :     }
     279                 :        140 :     d->a[i] = d->m[i];
     280         [ +  + ]:        140 :     if (--i <= 0) return NULL;
     281                 :        280 :   }
     282                 :            : }
     283                 :            : /* strictly increasing order [over integers] */
     284                 :            : static GEN
     285                 :    1173502 : _next_lt_i(forvec_t *d)
     286                 :            : {
     287                 :    1173502 :   long i = d->n;
     288         [ +  + ]:    1173502 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     289                 :            :   for (;;) {
     290         [ +  + ]:    1290009 :     if (cmpii(d->a[i], d->M[i]) < 0)
     291                 :            :     {
     292                 :    1159904 :       d->a[i] = incloop(d->a[i]);
     293                 :            :       /* m[i] < a[i] <= M[i] < M[i+1] */
     294         [ +  + ]:    1276397 :       while (i < d->n)
     295                 :            :       {
     296                 :            :         pari_sp av;
     297                 :            :         GEN t;
     298                 :     116493 :         i++;
     299         [ -  + ]:     116493 :         if (cmpii(d->a[i-1], d->a[i]) < 0) continue;
     300                 :     116493 :         av = avma;
     301                 :            :         /* M[i] > M[i-1] >= a[i-1] */
     302         [ -  + ]:     116493 :         t = addis(d->a[i-1],1); if (cmpii(t, d->m[i]) < 0) t = d->m[i];
     303                 :     116493 :         d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1]+1,m[i]) <= M[i]*/
     304                 :     116493 :         avma = av;
     305                 :            :       }
     306                 :    1159904 :       return (GEN)d->a;
     307                 :            :     }
     308                 :     130105 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     309         [ +  + ]:     130105 :     if (--i <= 0) return NULL;
     310                 :    1296808 :   }
     311                 :            : }
     312                 :            : /* strictly increasing order [generic] */
     313                 :            : static GEN
     314                 :         84 : _next_lt(forvec_t *d)
     315                 :            : {
     316                 :         84 :   long i = d->n;
     317         [ +  + ]:         84 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     318                 :            :   for (;;) {
     319                 :        133 :     d->a[i] = gaddgs(d->a[i], 1);
     320         [ +  + ]:        133 :     if (gcmp(d->a[i], d->M[i]) <= 0)
     321                 :            :     {
     322         [ +  + ]:         91 :       while (i < d->n)
     323                 :            :       {
     324                 :            :         GEN c;
     325                 :         35 :         i++;
     326         [ -  + ]:         35 :         if (gcmp(d->a[i-1], d->a[i]) < 0) continue;
     327                 :            :         /* M[i] > M[i-1] >= a[i-1] >= a[i] */
     328                 :         35 :         c = addis(gfloor(gsub(d->a[i-1], d->a[i])), 1); /* > a[i-1] - a[i] */
     329                 :         35 :         d->a[i] = gadd(d->a[i], c);
     330                 :            :         /* a[i-1] < a[i] <= M[i-1] + 1 => a[i] < M[i]+1 => a[i] <= M[i] */
     331                 :            :       }
     332                 :         56 :       return (GEN)d->a;
     333                 :            :     }
     334                 :         77 :     d->a[i] = d->m[i];
     335         [ +  + ]:         77 :     if (--i <= 0) return NULL;
     336                 :        147 :   }
     337                 :            : }
     338                 :            : /* for forvec(v=[],) */
     339                 :            : static GEN
     340                 :         14 : _next_void(forvec_t *d)
     341                 :            : {
     342         [ +  + ]:         14 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     343                 :         14 :   return NULL;
     344                 :            : }
     345                 :            : 
     346                 :            : /* Initialize minima (m) and maxima (M); guarantee M[i] - m[i] integer and
     347                 :            :  *   if flag = 1: m[i-1] <= m[i] <= M[i] <= M[i+1]
     348                 :            :  *   if flag = 2: m[i-1] <  m[i] <= M[i] <  M[i+1],
     349                 :            :  * for all i */
     350                 :            : int
     351                 :       7020 : forvec_init(forvec_t *d, GEN x, long flag)
     352                 :            : {
     353                 :       7020 :   long i, tx = typ(x), l = lg(x), t = t_INT;
     354         [ -  + ]:       7020 :   if (!is_vec_t(tx)) pari_err_TYPE("forvec [not a vector]", x);
     355                 :       7020 :   d->first = 1;
     356                 :       7020 :   d->n = l - 1;
     357                 :       7020 :   d->a = (GEN*)cgetg(l,tx);
     358                 :       7020 :   d->m = (GEN*)cgetg(l,tx);
     359                 :       7020 :   d->M = (GEN*)cgetg(l,tx);
     360         [ +  + ]:       7020 :   if (l == 1) { d->next = &_next_void; return 1; }
     361         [ +  + ]:      21179 :   for (i = 1; i < l; i++)
     362                 :            :   {
     363                 :      14194 :     GEN a, e = gel(x,i), m = gel(e,1), M = gel(e,2);
     364                 :      14194 :     tx = typ(e);
     365 [ +  + ][ +  + ]:      14194 :     if (! is_vec_t(tx) || lg(e)!=3)
     366                 :         14 :       pari_err_TYPE("forvec [expected vector not of type [min,MAX]]",e);
     367         [ +  + ]:      14180 :     if (typ(m) != t_INT) t = t_REAL;
     368         [ +  + ]:      14180 :     if (i > 1) switch(flag)
              [ +  +  + ]
     369                 :            :     {
     370                 :            :       case 1: /* a >= m[i-1] - m */
     371                 :         51 :         a = gceil(gsub(d->m[i-1], m));
     372         [ -  + ]:         51 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     373         [ +  + ]:         51 :         if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
     374                 :         51 :         break;
     375                 :            :       case 2: /* a > m[i-1] - m */
     376                 :       6848 :         a = gfloor(gsub(d->m[i-1], m));
     377         [ -  + ]:       6848 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     378                 :       6848 :         a = addis(a, 1);
     379         [ +  + ]:       6848 :         if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
     380                 :       6848 :         break;
     381                 :        282 :       default: m = gcopy(m);
     382                 :        282 :         break;
     383                 :            :     }
     384                 :      14180 :     M = gadd(m, gfloor(gsub(M,m))); /* ensure M-m is an integer */
     385         [ +  + ]:      14173 :     if (gcmp(m,M) > 0) { d->a = NULL; d->next = &_next; return 0; }
     386                 :      14166 :     d->m[i] = m;
     387                 :      14166 :     d->M[i] = M;
     388                 :            :   }
     389 [ +  + ][ +  + ]:       7036 :   if (flag == 1) for (i = l-2; i >= 1; i--)
     390                 :            :   {
     391                 :         51 :     GEN M = d->M[i], a = gfloor(gsub(d->M[i+1], M));
     392         [ -  + ]:         51 :     if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     393                 :            :     /* M[i]+a <= M[i+1] */
     394         [ +  + ]:         51 :     if (signe(a) < 0) d->M[i] = gadd(M, a);
     395                 :            :   }
     396 [ +  + ][ +  + ]:      13796 :   else if (flag == 2) for (i = l-2; i >= 1; i--)
     397                 :            :   {
     398                 :       6841 :     GEN M = d->M[i], a = gceil(gsub(d->M[i+1], M));
     399         [ -  + ]:       6841 :     if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     400                 :       6841 :     a = subiu(a, 1);
     401                 :            :     /* M[i]+a < M[i+1] */
     402         [ +  + ]:       6841 :     if (signe(a) < 0) d->M[i] = gadd(M, a);
     403                 :            :   }
     404         [ +  + ]:       6985 :   if (t == t_INT) {
     405         [ +  + ]:      21004 :     for (i = 1; i < l; i++) {
     406                 :      14054 :       d->a[i] = setloop(d->m[i]);
     407         [ -  + ]:      14054 :       if (typ(d->M[i]) != t_INT) d->M[i] = gfloor(d->M[i]);
     408                 :            :     }
     409                 :            :   } else {
     410         [ +  + ]:        140 :     for (i = 1; i < l; i++) d->a[i] = d->m[i];
     411                 :            :   }
     412   [ +  +  +  + ]:       6985 :   switch(flag)
     413                 :            :   {
     414         [ +  + ]:        135 :     case 0: d->next = t==t_INT? &_next_i:    &_next; break;
     415         [ +  + ]:         30 :     case 1: d->next = t==t_INT? &_next_le_i: &_next_le; break;
     416         [ +  + ]:       6813 :     case 2: d->next = t==t_INT? &_next_lt_i: &_next_lt; break;
     417                 :          7 :     default: pari_err_FLAG("forvec");
     418                 :            :   }
     419                 :       6992 :   return 1;
     420                 :            : }
     421                 :            : GEN
     422                 :    1304514 : forvec_next(forvec_t *d) { return d->next(d); }
     423                 :            : 
     424                 :            : void
     425                 :       6972 : forvec(GEN x, GEN code, long flag)
     426                 :            : {
     427                 :       6972 :   pari_sp av = avma;
     428                 :            :   forvec_t T;
     429                 :            :   GEN v;
     430         [ +  + ]:      13909 :   if (!forvec_init(&T, x, flag)) { avma = av; return; }
     431                 :       6937 :   push_lex((GEN)T.a, code);
     432         [ +  + ]:    1304345 :   while ((v = forvec_next(&T)))
     433                 :            :   {
     434                 :    1297408 :     closure_evalvoid(code);
     435         [ -  + ]:    1297408 :     if (loop_break()) break;
     436                 :            :   }
     437                 :       6937 :   pop_lex(1); avma = av;
     438                 :            : }
     439                 :            : 
     440                 :            : /********************************************************************/
     441                 :            : /**                                                                **/
     442                 :            : /**                              SUMS                              **/
     443                 :            : /**                                                                **/
     444                 :            : /********************************************************************/
     445                 :            : 
     446                 :            : GEN
     447                 :       1113 : somme(GEN a, GEN b, GEN code, GEN x)
     448                 :            : {
     449                 :       1113 :   pari_sp av, av0 = avma;
     450                 :            :   GEN p1;
     451                 :            : 
     452         [ -  + ]:       1113 :   if (typ(a) != t_INT) pari_err_TYPE("sum",a);
     453         [ +  + ]:       1113 :   if (!x) x = gen_0;
     454         [ -  + ]:       1113 :   if (gcmp(b,a) < 0) return gcopy(x);
     455                 :            : 
     456                 :       1113 :   b = gfloor(b);
     457                 :       1113 :   a = setloop(a);
     458                 :       1113 :   av=avma;
     459                 :       1113 :   push_lex(a,code);
     460                 :            :   for(;;)
     461                 :            :   {
     462                 :     867153 :     p1 = closure_evalnobrk(code);
     463         [ +  + ]:     867153 :     x=gadd(x,p1); if (cmpii(a,b) >= 0) break;
     464                 :     866040 :     a = incloop(a);
     465         [ -  + ]:     866040 :     if (gc_needed(av,1))
     466                 :            :     {
     467         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sum");
     468                 :          0 :       x = gerepileupto(av,x);
     469                 :            :     }
     470                 :     866040 :     set_lex(-1,a);
     471                 :     866040 :   }
     472                 :       1113 :   pop_lex(1); return gerepileupto(av0,x);
     473                 :            : }
     474                 :            : 
     475                 :            : GEN
     476                 :         14 : suminf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     477                 :            : {
     478                 :            :   long fl, G;
     479                 :         14 :   pari_sp av0 = avma, av;
     480                 :         14 :   GEN p1,x = real_1(prec);
     481                 :            : 
     482         [ -  + ]:         14 :   if (typ(a) != t_INT) pari_err_TYPE("suminf",a);
     483                 :         14 :   a = setloop(a);
     484                 :         14 :   av = avma;
     485                 :         14 :   fl=0; G = prec2nbits(prec) + 5;
     486                 :            :   for(;;)
     487                 :            :   {
     488                 :       8162 :     p1 = eval(E, a); x = gadd(x,p1); a = incloop(a);
     489 [ +  - ][ +  + ]:       8162 :     if (gequal0(p1) || gexpo(p1) <= gexpo(x)-G)
     490         [ +  + ]:         42 :       { if (++fl==3) break; }
     491                 :            :     else
     492                 :       8120 :       fl=0;
     493         [ -  + ]:       8148 :     if (gc_needed(av,1))
     494                 :            :     {
     495         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"suminf");
     496                 :          0 :       x = gerepileupto(av,x);
     497                 :            :     }
     498                 :       8148 :   }
     499                 :         14 :   return gerepileupto(av0, gaddgs(x,-1));
     500                 :            : }
     501                 :            : GEN
     502                 :         14 : suminf0(GEN a, GEN code, long prec)
     503                 :         14 : { EXPR_WRAP(code, suminf(EXPR_ARG, a, prec)); }
     504                 :            : 
     505                 :            : GEN
     506                 :         49 : sumdivexpr(GEN num, GEN code)
     507                 :            : {
     508                 :         49 :   pari_sp av = avma;
     509                 :         49 :   GEN y = gen_0, t = divisors(num);
     510                 :         49 :   long i, l = lg(t);
     511                 :            : 
     512                 :         49 :   push_lex(gen_0, code);
     513         [ +  + ]:       9289 :   for (i=1; i<l; i++)
     514                 :            :   {
     515                 :       9240 :     set_lex(-1,gel(t,i));
     516                 :       9240 :     y = gadd(y, closure_evalnobrk(code));
     517                 :            :   }
     518                 :         49 :   pop_lex(1); return gerepileupto(av,y);
     519                 :            : }
     520                 :            : GEN
     521                 :         42 : sumdivmultexpr(GEN num, GEN code)
     522                 :            : {
     523                 :         42 :   pari_sp av = avma;
     524                 :         42 :   GEN y = gen_1, P,E;
     525                 :         42 :   int isint = divisors_init(num, &P,&E);
     526                 :         42 :   long i, l = lg(P);
     527                 :            : 
     528         [ -  + ]:         42 :   if (l == 1) { avma = av; return gen_1; }
     529                 :         42 :   push_lex(gen_0, code);
     530         [ +  + ]:        196 :   for (i=1; i<l; i++)
     531                 :            :   {
     532                 :        154 :     GEN p = gel(P,i), q = p, z = gen_1;
     533                 :        154 :     long j, e = E[i];
     534 [ +  + ][ +  - ]:        560 :     for (j = 1; j <= e; j++, q = isint?mulii(q, p): gmul(q,p))
     535                 :            :     {
     536                 :        560 :       set_lex(-1, q);
     537                 :        560 :       z = gadd(z, closure_evalnobrk(code));
     538         [ +  + ]:        560 :       if (j == e) break;
     539                 :            :     }
     540                 :        154 :     y = gmul(y, z);
     541                 :            :   }
     542                 :         42 :   pop_lex(1); return gerepileupto(av,y);
     543                 :            : }
     544                 :            : 
     545                 :            : /********************************************************************/
     546                 :            : /**                                                                **/
     547                 :            : /**                           PRODUCTS                             **/
     548                 :            : /**                                                                **/
     549                 :            : /********************************************************************/
     550                 :            : 
     551                 :            : GEN
     552                 :     148288 : produit(GEN a, GEN b, GEN code, GEN x)
     553                 :            : {
     554                 :     148288 :   pari_sp av, av0 = avma;
     555                 :            :   GEN p1;
     556                 :            : 
     557         [ -  + ]:     148288 :   if (typ(a) != t_INT) pari_err_TYPE("prod",a);
     558         [ +  - ]:     148288 :   if (!x) x = gen_1;
     559         [ -  + ]:     148288 :   if (gcmp(b,a) < 0) return gcopy(x);
     560                 :            : 
     561                 :     148288 :   b = gfloor(b);
     562                 :     148288 :   a = setloop(a);
     563                 :     148288 :   av=avma;
     564                 :     148288 :   push_lex(a,code);
     565                 :            :   for(;;)
     566                 :            :   {
     567                 :     393981 :     p1 = closure_evalnobrk(code);
     568         [ +  + ]:     393981 :     x = gmul(x,p1); if (cmpii(a,b) >= 0) break;
     569                 :     245693 :     a = incloop(a);
     570         [ -  + ]:     245693 :     if (gc_needed(av,1))
     571                 :            :     {
     572         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prod");
     573                 :          0 :       x = gerepileupto(av,x);
     574                 :            :     }
     575                 :     245693 :     set_lex(-1,a);
     576                 :     245693 :   }
     577                 :     148288 :   pop_lex(1); return gerepileupto(av0,x);
     578                 :            : }
     579                 :            : 
     580                 :            : GEN
     581                 :          7 : prodinf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     582                 :            : {
     583                 :          7 :   pari_sp av0 = avma, av;
     584                 :            :   long fl,G;
     585                 :          7 :   GEN p1,x = real_1(prec);
     586                 :            : 
     587         [ -  + ]:          7 :   if (typ(a) != t_INT) pari_err_TYPE("prodinf",a);
     588                 :          7 :   a = setloop(a);
     589                 :          7 :   av = avma;
     590                 :          7 :   fl=0; G = -prec2nbits(prec)-5;
     591                 :            :   for(;;)
     592                 :            :   {
     593         [ -  + ]:        952 :     p1 = eval(E, a); if (gequal0(p1)) { x = p1; break; }
     594                 :        952 :     x = gmul(x,p1); a = incloop(a);
     595                 :        952 :     p1 = gsubgs(p1, 1);
     596 [ +  - ][ +  + ]:        952 :     if (gequal0(p1) || gexpo(p1) <= G) { if (++fl==3) break; } else fl=0;
                 [ +  + ]
     597         [ -  + ]:        945 :     if (gc_needed(av,1))
     598                 :            :     {
     599         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodinf");
     600                 :          0 :       x = gerepileupto(av,x);
     601                 :            :     }
     602                 :        945 :   }
     603                 :          7 :   return gerepilecopy(av0,x);
     604                 :            : }
     605                 :            : GEN
     606                 :          7 : prodinf1(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     607                 :            : {
     608                 :          7 :   pari_sp av0 = avma, av;
     609                 :            :   long fl,G;
     610                 :          7 :   GEN p1,p2,x = real_1(prec);
     611                 :            : 
     612         [ -  + ]:          7 :   if (typ(a) != t_INT) pari_err_TYPE("prodinf1",a);
     613                 :          7 :   a = setloop(a);
     614                 :          7 :   av = avma;
     615                 :          7 :   fl=0; G = -prec2nbits(prec)-5;
     616                 :            :   for(;;)
     617                 :            :   {
     618                 :        952 :     p2 = eval(E, a); p1 = gaddgs(p2,1);
     619         [ -  + ]:        952 :     if (gequal0(p1)) { x = p1; break; }
     620                 :        952 :     x = gmul(x,p1); a = incloop(a);
     621 [ +  - ][ +  + ]:        952 :     if (gequal0(p2) || gexpo(p2) <= G) { if (++fl==3) break; } else fl=0;
                 [ +  + ]
     622         [ -  + ]:        945 :     if (gc_needed(av,1))
     623                 :            :     {
     624         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodinf1");
     625                 :          0 :       x = gerepileupto(av,x);
     626                 :            :     }
     627                 :        945 :   }
     628                 :          7 :   return gerepilecopy(av0,x);
     629                 :            : }
     630                 :            : GEN
     631                 :         14 : prodinf0(GEN a, GEN code, long flag, long prec)
     632                 :            : {
     633      [ +  +  - ]:         14 :   switch(flag)
     634                 :            :   {
     635                 :          7 :     case 0: EXPR_WRAP(code, prodinf (EXPR_ARG, a, prec));
     636                 :          7 :     case 1: EXPR_WRAP(code, prodinf1(EXPR_ARG, a, prec));
     637                 :            :   }
     638                 :          0 :   pari_err_FLAG("prodinf");
     639                 :         14 :   return NULL; /* not reached */
     640                 :            : }
     641                 :            : 
     642                 :            : GEN
     643                 :          7 : prodeuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
     644                 :            : {
     645                 :          7 :   pari_sp av, av0 = avma;
     646                 :          7 :   GEN x = real_1(prec), prime;
     647                 :            :   forprime_t T;
     648                 :            : 
     649                 :          7 :   av = avma;
     650         [ -  + ]:          7 :   if (!forprime_init(&T, a,b)) { avma = av; return x; }
     651                 :            : 
     652                 :          7 :   av = avma;
     653         [ +  + ]:       8610 :   while ( (prime = forprime_next(&T)) )
     654                 :            :   {
     655                 :       8603 :     x = gmul(x, eval(E, prime));
     656         [ -  + ]:       8603 :     if (gc_needed(av,1))
     657                 :            :     {
     658         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodeuler");
     659                 :          0 :       x = gerepilecopy(av, x);
     660                 :            :     }
     661                 :            :   }
     662                 :          7 :   return gerepilecopy(av0,x);
     663                 :            : }
     664                 :            : GEN
     665                 :          7 : prodeuler0(GEN a, GEN b, GEN code, long prec)
     666                 :          7 : { EXPR_WRAP(code, prodeuler(EXPR_ARG, a, b, prec)); }
     667                 :            : GEN
     668                 :        133 : direuler0(GEN a, GEN b, GEN code, GEN c)
     669                 :        133 : { EXPR_WRAP(code, direuler(EXPR_ARG, a, b, c)); }
     670                 :            : 
     671                 :            : /********************************************************************/
     672                 :            : /**                                                                **/
     673                 :            : /**                       VECTORS & MATRICES                       **/
     674                 :            : /**                                                                **/
     675                 :            : /********************************************************************/
     676                 :            : 
     677                 :            : INLINE GEN
     678                 :    5687159 : copyupto(GEN z, GEN t)
     679                 :            : {
     680 [ +  + ][ +  + ]:    5687159 :   if (is_universal_constant(z) || (z>(GEN)pari_mainstack->bot && z<=t))
                 [ +  + ]
     681                 :    5682315 :     return z;
     682                 :            :   else
     683                 :    5687159 :     return gcopy(z);
     684                 :            : }
     685                 :            : 
     686                 :            : GEN
     687                 :     110824 : vecexpr0(GEN vec, GEN code, GEN pred)
     688                 :            : {
     689      [ +  +  - ]:     110824 :   switch(typ(vec))
     690                 :            :   {
     691                 :            :     case t_LIST:
     692                 :            :     {
     693         [ -  + ]:         14 :       if (list_typ(vec)==t_LIST_MAP)
     694                 :          0 :         vec = mapdomain_shallow(vec);
     695                 :            :       else
     696                 :         14 :         vec = list_data(vec);
     697         [ +  + ]:         14 :       if (!vec) return cgetg(1, t_VEC);
     698                 :          7 :       break;
     699                 :            :     }
     700                 :     110810 :     case t_VEC: case t_COL: case t_MAT: break;
     701                 :          0 :     default: pari_err_TYPE("[_|_<-_,_]",vec);
     702                 :            :   }
     703 [ +  + ][ +  - ]:     110817 :   if (pred && code)
     704                 :        252 :     EXPR_WRAP(code,vecselapply((void*)pred,&gp_evalbool,EXPR_ARGUPTO,vec))
     705         [ +  - ]:     110565 :   else if (code)
     706                 :     110565 :     EXPR_WRAP(code,vecapply(EXPR_ARGUPTO,vec))
     707                 :            :   else
     708                 :     110824 :     EXPR_WRAP(pred,vecselect(EXPR_ARGBOOL,vec))
     709                 :            : }
     710                 :            : 
     711                 :            : GEN
     712                 :        175 : vecexpr1(GEN vec, GEN code, GEN pred)
     713                 :            : {
     714                 :        175 :   GEN v = vecexpr0(vec, code, pred);
     715         [ +  + ]:        175 :   return lg(v) == 1? v: shallowconcat1(v);
     716                 :            : }
     717                 :            : 
     718                 :            : GEN
     719                 :    2338628 : vecteur(GEN nmax, GEN code)
     720                 :            : {
     721                 :            :   GEN y,p1;
     722                 :            :   long i,m;
     723                 :    2338628 :   GEN c=utoipos(1);
     724                 :            : 
     725                 :    2338628 :   m = gtos(nmax);
     726         [ +  + ]:    2338628 :   if (m < 0)  pari_err_DOMAIN("vector", "dimension", "<", gen_0, stoi(m));
     727         [ +  + ]:    2338614 :   if (!code) return zerovec(m);
     728                 :      18639 :   y = cgetg(m+1,t_VEC); push_lex(c, code);
     729         [ +  + ]:    5095503 :   for (i=1; i<=m; i++)
     730                 :            :   {
     731                 :    5076871 :     c[2] = i; p1 = closure_evalnobrk(code);
     732                 :    5076864 :     gel(y,i) = copyupto(p1, y);
     733                 :    5076864 :     set_lex(-1,c);
     734                 :            :   }
     735                 :    2338607 :   pop_lex(1); return y;
     736                 :            : }
     737                 :            : 
     738                 :            : GEN
     739                 :         42 : vecteursmall(GEN nmax, GEN code)
     740                 :            : {
     741                 :            :   GEN y;
     742                 :            :   long i,m;
     743                 :         42 :   GEN c=utoipos(1);
     744                 :            : 
     745                 :         42 :   m = gtos(nmax);
     746         [ +  + ]:         42 :   if (m < 0)  pari_err_DOMAIN("vectorsmall", "dimension", "<", gen_0, stoi(m));
     747         [ +  + ]:         35 :   if (!code) return zero_zv(m);
     748                 :         21 :   y = cgetg(m+1,t_VECSMALL); push_lex(c,code);
     749         [ +  + ]:         84 :   for (i=1; i<=m; i++)
     750                 :            :   {
     751                 :         70 :     c[2] = i;
     752                 :         70 :     y[i] = gtos(closure_evalnobrk(code));
     753                 :         63 :     set_lex(-1,c);
     754                 :            :   }
     755                 :         28 :   pop_lex(1); return y;
     756                 :            : }
     757                 :            : 
     758                 :            : GEN
     759                 :        483 : vvecteur(GEN nmax, GEN n)
     760                 :            : {
     761                 :        483 :   GEN y = vecteur(nmax,n);
     762                 :        476 :   settyp(y,t_COL); return y;
     763                 :            : }
     764                 :            : 
     765                 :            : GEN
     766                 :        511 : matrice(GEN nlig, GEN ncol, GEN code)
     767                 :            : {
     768                 :            :   GEN y, z, p1;
     769                 :            :   long i, j, m, n;
     770                 :        511 :   GEN c1 = utoipos(1), c2 = utoipos(1);
     771                 :            : 
     772                 :        511 :   m = gtos(ncol);
     773                 :        511 :   n = gtos(nlig);
     774         [ +  + ]:        511 :   if (m < 0)  pari_err_DOMAIN("matrix", "nbcols", "<", gen_0, stoi(m));
     775         [ +  + ]:        504 :   if (n < 0)  pari_err_DOMAIN("matrix", "nbrows", "<", gen_0, stoi(n));
     776         [ +  + ]:        497 :   if (!m) return cgetg(1,t_MAT);
     777 [ +  + ][ -  + ]:        455 :   if (!code || !n) return zeromatcopy(n, m);
     778                 :        371 :   push_lex(c1,code);
     779                 :        371 :   push_lex(c2,NULL); y = cgetg(m+1,t_MAT);
     780         [ +  + ]:       5572 :   for (i=1; i<=m; i++)
     781                 :            :   {
     782                 :       5201 :     c2[2] = i; z = cgetg(n+1,t_COL); gel(y,i) = z;
     783         [ +  + ]:     615496 :     for (j=1; j<=n; j++)
     784                 :            :     {
     785                 :     610295 :       c1[2] = j; p1 = closure_evalnobrk(code);
     786                 :     610295 :       gel(z,j) = copyupto(p1, y);
     787                 :     610295 :       set_lex(-2,c1);
     788                 :     610295 :       set_lex(-1,c2);
     789                 :            :     }
     790                 :            :   }
     791                 :        497 :   pop_lex(2); return y;
     792                 :            : }
     793                 :            : 
     794                 :            : /********************************************************************/
     795                 :            : /**                                                                **/
     796                 :            : /**                         SUMMING SERIES                         **/
     797                 :            : /**                                                                **/
     798                 :            : /********************************************************************/
     799                 :            : /* h = (2+2x)g'- g; g has t_INT coeffs */
     800                 :            : static GEN
     801                 :       1246 : delt(GEN g, long n)
     802                 :            : {
     803                 :       1246 :   GEN h = cgetg(n+3,t_POL);
     804                 :            :   long k;
     805                 :       1246 :   h[1] = g[1];
     806                 :       1246 :   gel(h,2) = gel(g,2);
     807         [ +  + ]:     359597 :   for (k=1; k<n; k++)
     808                 :     358351 :     gel(h,k+2) = addii(mului(k+k+1,gel(g,k+2)), mului(k<<1,gel(g,k+1)));
     809                 :       1246 :   gel(h,n+2) = mului(n<<1, gel(g,n+1)); return h;
     810                 :            : }
     811                 :            : 
     812                 :            : #ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */
     813                 :            : #pragma optimize("g",off)
     814                 :            : #endif
     815                 :            : /* P = polzagier(n,m)(-X), unnormalized (P(0) != 1) */
     816                 :            : static GEN
     817                 :         49 : polzag1(long n, long m)
     818                 :            : {
     819                 :         49 :   const long d = n - m, d2 = d<<1, r = (m+1)>>1, D = (d+1)>>1;
     820                 :            :   long i, k;
     821                 :         49 :   pari_sp av = avma;
     822                 :            :   GEN g, T;
     823                 :            : 
     824 [ +  - ][ -  + ]:         49 :   if (d <= 0 || m < 0) return pol_0(0);
     825                 :         49 :   g = cgetg(d+2, t_POL);
     826                 :         49 :   g[1] = evalsigne(1)|evalvarn(0);
     827                 :         49 :   T = cgetg(d+1,t_VEC);
     828                 :            :   /* T[k+1] = binomial(2d,2k+1), 0 <= k < d */
     829                 :         49 :   gel(T,1) = utoipos(d2);
     830         [ +  + ]:       1267 :   for (k = 1; k < D; k++)
     831                 :            :   {
     832                 :       1218 :     long k2 = k<<1;
     833                 :       1218 :     gel(T,k+1) = diviiexact(mulii(gel(T,k), muluu(d2-k2+1, d2-k2)),
     834                 :       1218 :                             muluu(k2,k2+1));
     835                 :            :   }
     836         [ +  + ]:       1281 :   for (; k < d; k++) gel(T,k+1) = gel(T,d-k);
     837                 :         49 :   gel(g,2) = gel(T,d); /* binomial(2d, 2(d-1)+1) */
     838         [ +  + ]:       2499 :   for (i = 1; i < d; i++)
     839                 :            :   {
     840                 :       2450 :     pari_sp av2 = avma;
     841                 :       2450 :     GEN s, t = gel(T,d-i); /* binomial(2d, 2(d-1-i)+1) */
     842                 :       2450 :     s = t;
     843         [ +  + ]:     180082 :     for (k = d-i; k < d; k++)
     844                 :            :     {
     845                 :     177632 :       long k2 = k<<1;
     846                 :     177632 :       t = diviiexact(mulii(t, muluu(d2-k2+1, d-k)), muluu(k2+1,k-(d-i)+1));
     847                 :     177632 :       s = addii(s, t);
     848                 :            :     }
     849                 :            :     /* g_i = sum_{d-1-i <= k < d}, binomial(2*d, 2*k+1)*binomial(k,d-1-i) */
     850                 :       2450 :     gel(g,i+2) = gerepileuptoint(av2, s);
     851                 :            :   }
     852                 :            :   /* sum_{0 <= i < d} g_i x^i * (x+x^2)^r */
     853                 :         49 :   g = RgX_mulXn(gmul(g, gpowgs(deg1pol(gen_1,gen_1,0),r)), r);
     854         [ +  + ]:         49 :   if (!odd(m)) g = delt(g, n);
     855         [ +  + ]:       1274 :   for (i=1; i<=r; i++)
     856                 :            :   {
     857                 :       1225 :     g = delt(ZX_deriv(g), n);
     858         [ -  + ]:       1225 :     if (gc_needed(av,4))
     859                 :            :     {
     860         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"polzag, i = %ld/%ld", i,r);
     861                 :          0 :       g = gerepilecopy(av, g);
     862                 :            :     }
     863                 :            :   }
     864                 :         49 :   return g;
     865                 :            : }
     866                 :            : GEN
     867                 :         28 : polzag(long n, long m)
     868                 :            : {
     869                 :         28 :   pari_sp av = avma;
     870                 :         28 :   GEN g = ZX_unscale(polzag1(n,m), gen_m1);
     871                 :         28 :   return gerepileupto(av, RgX_Rg_div(g,gel(g,2)));
     872                 :            : }
     873                 :            : 
     874                 :            : GEN
     875                 :          7 : sumalt(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     876                 :            : {
     877                 :            :   ulong k, N;
     878                 :          7 :   pari_sp av = avma, av2;
     879                 :            :   GEN s, az, c, d;
     880                 :            : 
     881         [ -  + ]:          7 :   if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
     882                 :          7 :   N = (ulong)(0.39322*(prec2nbits(prec) + 7)); /*0.39322 > 1/log_2(3+sqrt(8))*/
     883                 :          7 :   d = powru(addsr(3, sqrtr(stor(8,prec))), N);
     884                 :          7 :   d = shiftr(addrr(d, invr(d)),-1);
     885                 :          7 :   a = setloop(a);
     886                 :          7 :   az = gen_m1; c = d;
     887                 :          7 :   s = gen_0;
     888                 :          7 :   av2 = avma;
     889                 :          7 :   for (k=0; ; k++) /* k < N */
     890                 :            :   {
     891                 :        371 :     c = addir(az,c); s = gadd(s, gmul(c, eval(E, a)));
     892         [ +  + ]:        371 :     if (k==N-1) break;
     893                 :        364 :     az = diviuuexact(muluui((N-k)<<1,N+k,az), k+1, (k<<1)+1);
     894                 :        364 :     a = incloop(a); /* in place! */
     895         [ -  + ]:        364 :     if (gc_needed(av,4))
     896                 :            :     {
     897         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumalt, k = %ld/%ld", k,N-1);
     898                 :          0 :       gerepileall(av2, 3, &az,&c,&s);
     899                 :            :     }
     900                 :        364 :   }
     901                 :          7 :   return gerepileupto(av, gdiv(s,d));
     902                 :            : }
     903                 :            : 
     904                 :            : GEN
     905                 :          7 : sumalt2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     906                 :            : {
     907                 :            :   long k, N;
     908                 :          7 :   pari_sp av = avma, av2;
     909                 :            :   GEN s, dn, pol;
     910                 :            : 
     911         [ -  + ]:          7 :   if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
     912                 :          7 :   N = (long)(0.307073*(prec2nbits(prec) + 5)); /*0.307073 > 1/log_2(\beta_B)*/
     913                 :          7 :   pol = ZX_div_by_X_1(polzag1(N,N>>1), &dn);
     914                 :          7 :   a = setloop(a);
     915                 :          7 :   N = degpol(pol);
     916                 :          7 :   s = gen_0;
     917                 :          7 :   av2 = avma;
     918         [ +  - ]:        280 :   for (k=0; k<=N; k++)
     919                 :            :   {
     920                 :        280 :     GEN t = itor(gel(pol,k+2), prec+EXTRAPRECWORD);
     921                 :        280 :     s = gadd(s, gmul(t, eval(E, a)));
     922         [ +  + ]:        280 :     if (k == N) break;
     923                 :        273 :     a = incloop(a); /* in place! */
     924         [ -  + ]:        273 :     if (gc_needed(av,4))
     925                 :            :     {
     926         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumalt2, k = %ld/%ld", k,N-1);
     927                 :          0 :       s = gerepileupto(av2, s);
     928                 :            :     }
     929                 :            :   }
     930                 :          7 :   return gerepileupto(av, gdiv(s,dn));
     931                 :            : }
     932                 :            : 
     933                 :            : GEN
     934                 :         14 : sumalt0(GEN a, GEN code, long flag, long prec)
     935                 :            : {
     936      [ +  +  - ]:         14 :   switch(flag)
     937                 :            :   {
     938                 :          7 :     case 0: EXPR_WRAP(code, sumalt (EXPR_ARG,a,prec));
     939                 :          7 :     case 1: EXPR_WRAP(code, sumalt2(EXPR_ARG,a,prec));
     940                 :          0 :     default: pari_err_FLAG("sumalt");
     941                 :            :   }
     942                 :         14 :   return NULL; /* not reached */
     943                 :            : }
     944                 :            : 
     945                 :            : /* For k > 0, set S[k*2^i] <- g(k*2^i), k*2^i <= N = #S.
     946                 :            :  * Only needed with k odd (but also works for g even). */
     947                 :            : static void
     948                 :       7406 : binsum(GEN S, ulong k, void *E, GEN (*f)(void *, GEN), GEN a,
     949                 :            :         long G, long prec)
     950                 :            : {
     951                 :       7406 :   long e, i, N = lg(S)-1, l = expu(N / k); /* k 2^l <= N < k 2^(l+1) */
     952                 :            :   pari_sp av;
     953                 :       7406 :   GEN r, t = gen_0;
     954                 :            : 
     955                 :       7406 :   gel(S, k << l) = cgetr(prec); av = avma;
     956                 :       7406 :   G -= l;
     957                 :       7406 :   r = utoipos(k<<l);
     958                 :       7406 :   for(e=0;;e++) /* compute g(k 2^l) with absolute error ~ 2^(G-l) */
     959                 :            :   {
     960                 :    3964471 :     GEN u = gtofp(f(E, addii(a,r)), prec);
     961         [ -  + ]:    3964471 :     if (typ(u) != t_REAL) pari_err_TYPE("sumpos",u);
     962         [ +  + ]:    3964471 :     if (!signe(u)) break;
     963         [ +  + ]:    3964282 :     if (!e)
     964                 :       7217 :       t = u;
     965                 :            :     else {
     966                 :    3957065 :       shiftr_inplace(u, e);
     967                 :    3957065 :       t = addrr(t,u);
     968         [ +  + ]:    3957065 :       if (expo(u) < G) break;
     969                 :            :     }
     970                 :    3957065 :     r = shifti(r,1);
     971                 :    3957065 :   }
     972                 :       7406 :   gel(S, k << l) = t = gerepileuptoleaf(av, t);
     973                 :            :   /* g(j) = 2g(2j) + f(a+j) for all j > 0 */
     974         [ +  + ]:      14812 :   for(i = l-1; i >= 0; i--)
     975                 :            :   { /* t ~ g(2 * k*2^i) with error ~ 2^(G-i-1) */
     976                 :            :     GEN u;
     977                 :       7406 :     av = avma; u = gtofp(f(E, addiu(a, k << i)), prec);
     978         [ -  + ]:       7406 :     if (typ(u) != t_REAL) pari_err_TYPE("sumpos",u);
     979                 :       7406 :     t = addrr(gtofp(u,prec), shiftr(t,1)); /* ~ g(k*2^i) */
     980                 :       7406 :     gel(S, k << i) = t = gerepileuptoleaf(av, t);
     981                 :            :   }
     982                 :       7406 : }
     983                 :            : /* For k > 0, let g(k) := \sum_{e >= 0} 2^e f(a + k*2^e).
     984                 :            :  * Return [g(k), 1 <= k <= N] */
     985                 :            : static GEN
     986                 :         70 : sumpos_init(void *E, GEN (*f)(void *, GEN), GEN a, long N, long prec)
     987                 :            : {
     988                 :         70 :   GEN S = cgetg(N+1,t_VEC);
     989                 :         70 :   long k, G = -prec2nbits(prec) - 5;
     990         [ +  + ]:       7476 :   for (k=1; k<=N; k+=2) binsum(S,k, E,f, a,G,prec);
     991                 :         70 :   return S;
     992                 :            : }
     993                 :            : 
     994                 :            : GEN
     995                 :         56 : sumpos(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     996                 :            : {
     997                 :            :   ulong k, N;
     998                 :         56 :   pari_sp av = avma;
     999                 :            :   GEN s, az, c, d, S;
    1000                 :            : 
    1001         [ -  + ]:         56 :   if (typ(a) != t_INT) pari_err_TYPE("sumpos",a);
    1002                 :         56 :   a = subiu(a,1);
    1003                 :         56 :   N = (ulong)(0.4*(prec2nbits(prec) + 7));
    1004         [ -  + ]:         56 :   if (odd(N)) N++; /* extra precision for free */
    1005                 :         56 :   d = powru(addsr(3, sqrtr(stor(8,prec))), N);
    1006                 :         56 :   d = shiftr(addrr(d, invr(d)),-1);
    1007                 :         56 :   az = gen_m1; c = d;
    1008                 :            : 
    1009                 :         56 :   S = sumpos_init(E, eval, a, N, prec);
    1010                 :         56 :   s = gen_0;
    1011         [ +  - ]:      10360 :   for (k=0; k<N; k++)
    1012                 :            :   {
    1013                 :            :     GEN t;
    1014                 :      10360 :     c = addir(az,c);
    1015                 :      10360 :     t = mulrr(gel(S,k+1), c);
    1016         [ +  + ]:      10360 :     s = odd(k)? mpsub(s, t): mpadd(s, t);
    1017         [ +  + ]:      10360 :     if (k == N-1) break;
    1018                 :      10304 :     az = diviuuexact(muluui((N-k)<<1,N+k,az), k+1, (k<<1)+1);
    1019                 :            :   }
    1020                 :         56 :   return gerepileupto(av, gdiv(s,d));
    1021                 :            : }
    1022                 :            : 
    1023                 :            : GEN
    1024                 :         14 : sumpos2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1025                 :            : {
    1026                 :            :   ulong k, N;
    1027                 :         14 :   pari_sp av = avma;
    1028                 :            :   GEN s, pol, dn, S;
    1029                 :            : 
    1030         [ -  + ]:         14 :   if (typ(a) != t_INT) pari_err_TYPE("sumpos2",a);
    1031                 :         14 :   a = subiu(a,1);
    1032                 :         14 :   N = (ulong)(0.31*(prec2nbits(prec) + 5));
    1033                 :            : 
    1034         [ -  + ]:         14 :   if (odd(N)) N++; /* extra precision for free */
    1035                 :         14 :   S = sumpos_init(E, eval, a, N, prec);
    1036                 :         14 :   pol = ZX_div_by_X_1(polzag1(N,N>>1), &dn);
    1037                 :         14 :   s = gen_0;
    1038         [ +  + ]:       4466 :   for (k=0; k<N; k++)
    1039                 :            :   {
    1040                 :       4452 :     GEN t = mulri(gel(S,k+1), gel(pol,k+2));
    1041         [ +  + ]:       4452 :     s = odd(k)? mpsub(s,t): mpadd(s,t);
    1042                 :            :   }
    1043                 :         14 :   return gerepileupto(av, gdiv(s,dn));
    1044                 :            : }
    1045                 :            : 
    1046                 :            : GEN
    1047                 :         70 : sumpos0(GEN a, GEN code, long flag, long prec)
    1048                 :            : {
    1049      [ +  +  - ]:         70 :   switch(flag)
    1050                 :            :   {
    1051                 :         56 :     case 0: EXPR_WRAP(code, sumpos (EXPR_ARG,a,prec));
    1052                 :         14 :     case 1: EXPR_WRAP(code, sumpos2(EXPR_ARG,a,prec));
    1053                 :          0 :     default: pari_err_FLAG("sumpos");
    1054                 :            :   }
    1055                 :         70 :   return NULL; /* not reached */
    1056                 :            : }
    1057                 :            : 
    1058                 :            : /********************************************************************/
    1059                 :            : /**                                                                **/
    1060                 :            : /**            SEARCH FOR REAL ZEROS of an expression              **/
    1061                 :            : /**                                                                **/
    1062                 :            : /********************************************************************/
    1063                 :            : /* Brent's method, [a,b] bracketing interval */
    1064                 :            : GEN
    1065                 :       1043 : zbrent(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
    1066                 :            : {
    1067                 :            :   long sig, iter, itmax;
    1068                 :       1043 :   pari_sp av = avma;
    1069                 :            :   GEN c, d, e, tol, fa, fb, fc;
    1070                 :            : 
    1071 [ +  + ][ -  + ]:       1043 :   if (typ(a) != t_REAL || realprec(a) < prec) a = gtofp(a, prec);
    1072 [ +  + ][ -  + ]:       1043 :   if (typ(b) != t_REAL || realprec(b) < prec) b = gtofp(b, prec);
    1073                 :       1043 :   sig = cmprr(b, a);
    1074         [ -  + ]:       1043 :   if (!sig) return gerepileupto(av, a);
    1075         [ -  + ]:       1043 :   if (sig < 0) {c = a; a = b; b = c;} else c = b;
    1076                 :       1043 :   fa = eval(E, a);
    1077                 :       1043 :   fb = eval(E, b);
    1078         [ +  + ]:       1043 :   if (gsigne(fa)*gsigne(fb) > 0)
    1079                 :          7 :     pari_err_DOMAIN("solve", "f(a)f(b)", ">", gen_0, mkvec2(fa, fb));
    1080                 :       1036 :   itmax = prec2nbits(prec) * 2 + 1;
    1081                 :       1036 :   tol = real2n(5-prec2nbits(prec), LOWDEFAULTPREC);
    1082                 :       1036 :   fc = fb;
    1083                 :       1036 :   e = d = NULL; /* gcc -Wall */
    1084         [ +  - ]:      15810 :   for (iter = 1; iter <= itmax; ++iter)
    1085                 :            :   {
    1086                 :            :     GEN xm, tol1;
    1087         [ +  + ]:      15810 :     if (gsigne(fb)*gsigne(fc) > 0)
    1088                 :            :     {
    1089                 :       8675 :       c = a; fc = fa; e = d = subrr(b, a);
    1090                 :            :     }
    1091         [ +  + ]:      15810 :     if (gcmp(gabs(fc, 0), gabs(fb, 0)) < 0)
    1092                 :            :     {
    1093                 :       4147 :       a = b; b = c; c = a; fa = fb; fb = fc; fc = fa;
    1094                 :            :     }
    1095         [ -  + ]:      15810 :     tol1 = absr_cmp(tol, b) > 0? sqrr(tol): mulrr(tol, absr(b));
    1096                 :      15810 :     xm = shiftr(subrr(c, b), -1);
    1097 [ +  + ][ +  + ]:      15810 :     if (absr_cmp(xm, tol1) <= 0 || gequal0(fb)) break; /* SUCCESS */
    1098                 :            : 
    1099 [ +  - ][ +  + ]:      14774 :     if (absr_cmp(e, tol1) >= 0 && gcmp(gabs(fa, 0), gabs(fb, 0)) > 0)
    1100                 :      13543 :     { /* attempt interpolation */
    1101                 :      13543 :       GEN min1, min2, p, q, s = gdiv(fb, fa);
    1102         [ +  + ]:      13543 :       if (cmprr(a, c) == 0)
    1103                 :            :       {
    1104                 :       8285 :         p = gmul2n(gmul(xm, s), 1);
    1105                 :       8285 :         q = gsubsg(1, s);
    1106                 :            :       }
    1107                 :            :       else
    1108                 :            :       {
    1109                 :       5258 :         GEN r = gdiv(fb, fc);
    1110                 :       5258 :         q = gdiv(fa, fc);
    1111                 :       5258 :         p = gmul2n(gmul(gsub(q, r), gmul(xm, q)), 1);
    1112                 :       5258 :         p = gmul(s, gsub(p, gmul(gsub(b, a), gsubgs(r, 1))));
    1113                 :       5258 :         q = gmul(gmul(gsubgs(q, 1), gsubgs(r, 1)), gsubgs(s, 1));
    1114                 :            :       }
    1115         [ +  + ]:      13543 :       if (gsigne(p) > 0) q = gneg_i(q); else p = gneg_i(p);
    1116                 :      13543 :       min1 = gsub(gmulsg(3, gmul(xm,q)), gabs(gmul(q, tol1), 0));
    1117                 :      13543 :       min2 = gabs(gmul(e, q), 0);
    1118         [ +  + ]:      13543 :       if (gcmp(gmul2n(p, 1), gmin(min1, min2)) < 0)
    1119                 :      11698 :         { e = d; d = gdiv(p, q); } /* interpolation OK */
    1120                 :            :       else
    1121                 :       1845 :         { d = xm; e = d; } /* failed, use bisection */
    1122                 :            :     }
    1123                 :       1231 :     else { d = xm; e = d; } /* bound decreasing too slowly, use bisection */
    1124                 :      14774 :     a = b; fa = fb;
    1125         [ +  + ]:      14774 :     if (gcmp(gabs(d, 0), tol1) > 0) b = gadd(b, d);
    1126         [ +  + ]:        998 :     else if (gsigne(xm) > 0)      b = addrr(b, tol1);
    1127                 :        573 :     else                          b = subrr(b, tol1);
    1128         [ +  + ]:      14774 :     if (realprec(b) < prec) b = rtor(b, prec);
    1129                 :      14774 :     fb = eval(E, b);
    1130                 :            :   }
    1131         [ -  + ]:       1036 :   if (iter > itmax) pari_err_IMPL("solve recovery [too many iterations]");
    1132                 :       1036 :   return gerepileuptoleaf(av, rcopy(b));
    1133                 :            : }
    1134                 :            : 
    1135                 :            : GEN
    1136                 :         21 : zbrent0(GEN a, GEN b, GEN code, long prec)
    1137                 :         21 : { EXPR_WRAP(code, zbrent(EXPR_ARG, a, b, prec)); }
    1138                 :            : 
    1139                 :            : /* x = solve_start(&D, a, b, prec)
    1140                 :            :  * while (x) {
    1141                 :            :  *   y = ...(x);
    1142                 :            :  *   x = solve_next(&D, y);
    1143                 :            :  * }
    1144                 :            :  * return D.res; */
    1145                 :            : 
    1146                 :            : /* Find zeros of a function in the real interval [a,b] by interval splitting */
    1147                 :            : GEN
    1148                 :         84 : solvestep(void *E, GEN (*f)(void *,GEN), GEN a, GEN b, GEN step, long flag, long prec)
    1149                 :            : {
    1150                 :         84 :   const long ITMAX = 10;
    1151                 :         84 :   pari_sp av = avma;
    1152                 :            :   GEN fa, ainit, binit;
    1153                 :         84 :   long sainit, it, bit = prec2nbits(prec) / 2, ct = 0, s = gcmp(a,b);
    1154                 :            : 
    1155 [ -  + ][ #  # ]:         84 :   if (!s) return gequal0(f(E, a)) ? gcopy(mkvec(a)): cgetg(1,t_VEC);
    1156         [ -  + ]:         84 :   if (s > 0) swap(a, b);
    1157         [ +  + ]:         84 :   if (flag&4)
    1158                 :            :   {
    1159         [ -  + ]:         77 :     if (gcmpgs(step,1)<=0) pari_err_DOMAIN("solvestep","step","<=",gen_1,step);
    1160         [ -  + ]:         77 :     if (gsigne(a) <= 0) pari_err_DOMAIN("solvestep","a","<=",gen_0,a);
    1161                 :            :   }
    1162         [ -  + ]:          7 :   else if (gsigne(step) <= 0)
    1163                 :          0 :     pari_err_DOMAIN("solvestep","step","<=",gen_0,step);
    1164                 :         84 :   ainit = a = gtofp(a, prec); fa = f(E, a);
    1165                 :         84 :   binit = b = gtofp(b, prec); step = gtofp(step, prec);
    1166                 :         84 :   sainit = gsigne(fa);
    1167         [ -  + ]:         84 :   if (gexpo(fa) < -bit) sainit = 0;
    1168         [ +  - ]:         91 :   for (it = 0; it < ITMAX; it++)
    1169                 :            :   {
    1170                 :         91 :     pari_sp av2 = avma;
    1171                 :         91 :     GEN v = cgetg(1, t_VEC);
    1172                 :            :     long sa;
    1173                 :         91 :     a = ainit;
    1174                 :         91 :     b = binit;
    1175                 :         91 :     sa = sainit;
    1176         [ +  + ]:       2317 :     while (gcmp(a,b) < 0)
    1177                 :            :     {
    1178         [ +  + ]:       2226 :       GEN fc, c = (flag&4)? gmul(a, step): gadd(a, step);
    1179                 :            :       long sc;
    1180         [ +  + ]:       2226 :       if (gcmp(c,b) > 0) c = b;
    1181                 :       2226 :       fc = f(E, c);
    1182                 :       2226 :       sc = gsigne(fc);
    1183         [ +  + ]:       2226 :       if (gexpo(fc) < -bit) sc = 0;
    1184 [ +  + ][ +  + ]:       2226 :       if (!sc || sa*sc < 0)
    1185                 :            :       {
    1186                 :            :         long e;
    1187                 :            :         GEN z;
    1188         [ +  + ]:        434 :         z = sc? zbrent(E, f, a, c, prec): c;
    1189                 :        434 :         (void)grndtoi(z, &e);
    1190         [ +  + ]:        434 :         if (e  <= -bit) ct = 1;
    1191 [ -  + ][ #  # ]:        434 :         if ((flag&1) && ((!(flag&8)) || ct)) return gerepileupto(av, z);
                 [ #  # ]
    1192                 :        434 :         v = gconcat(v, z);
    1193                 :            :       }
    1194                 :       2226 :       a = c; fa = fc; sa = sc;
    1195                 :            :     }
    1196 [ +  + ][ +  - ]:         91 :     if ((!(flag&2) || lg(v) > 1) && (!(flag&8) || ct))
         [ +  + ][ +  + ]
    1197                 :         84 :       return gerepilecopy(av, v);
    1198         [ +  - ]:          7 :     step = (flag&4)? sqrtr(sqrtr(step)): gmul2n(step, -2);
    1199                 :          7 :     gerepileall(av2, 2, &fa, &step);
    1200                 :            :   }
    1201         [ #  # ]:          0 :   if (it == ITMAX) pari_err_IMPL("solvestep recovery [too many iterations]");
    1202                 :         84 :   return NULL;
    1203                 :            : }
    1204                 :            : 
    1205                 :            : GEN
    1206                 :          7 : solvestep0(GEN a, GEN b, GEN step, GEN code, long flag, long prec)
    1207                 :          7 : { EXPR_WRAP(code, solvestep(EXPR_ARG, a,b, step, flag, prec)); }
    1208                 :            : 
    1209                 :            : /********************************************************************/
    1210                 :            : /**                     Numerical derivation                       **/
    1211                 :            : /********************************************************************/
    1212                 :            : 
    1213                 :            : struct deriv_data
    1214                 :            : {
    1215                 :            :   GEN code;
    1216                 :            :   GEN args;
    1217                 :            : };
    1218                 :            : 
    1219                 :         91 : static GEN deriv_eval(void *E, GEN x, long prec)
    1220                 :            : {
    1221                 :         91 :  struct deriv_data *data=(struct deriv_data *)E;
    1222                 :         91 :  gel(data->args,1)=x;
    1223                 :         91 :  return closure_callgenvecprec(data->code, data->args, prec);
    1224                 :            : }
    1225                 :            : 
    1226                 :            : /* Rationale: (f(2^-e) - f(-2^-e) + O(2^-pr)) / (2 * 2^-e) = f'(0) + O(2^-2e)
    1227                 :            :  * since 2nd derivatives cancel.
    1228                 :            :  *   prec(LHS) = pr - e
    1229                 :            :  *   prec(RHS) = 2e, equal when  pr = 3e = 3/2 fpr (fpr = required final prec)
    1230                 :            :  *
    1231                 :            :  * For f'(x), x far from 0: prec(LHS) = pr - e - expo(x)
    1232                 :            :  * --> pr = 3/2 fpr + expo(x) */
    1233                 :            : GEN
    1234                 :       1379 : derivnum(void *E, GEN (*eval)(void *, GEN, long), GEN x, long prec)
    1235                 :            : {
    1236                 :            :   GEN eps,a,b, y;
    1237                 :            :   long pr, l, e, ex, newprec;
    1238                 :       1379 :   pari_sp av = avma;
    1239                 :       1379 :   long p = precision(x);
    1240         [ +  + ]:       1379 :   long fpr = p ? prec2nbits(p): prec2nbits(prec);
    1241                 :       1379 :   ex = gexpo(x);
    1242         [ +  + ]:       1379 :   if (ex < 0) ex = 0; /* near 0 */
    1243                 :       1379 :   pr = (long)ceil(fpr * 1.5 + ex);
    1244                 :       1379 :   l = nbits2prec(pr);
    1245                 :       1379 :   newprec = nbits2prec(pr + ex + BITS_IN_LONG);
    1246         [ +  + ]:       1379 :   switch(typ(x))
    1247                 :            :   {
    1248                 :            :     case t_REAL:
    1249                 :            :     case t_COMPLEX:
    1250                 :        784 :       x = gprec_w(x, newprec);
    1251                 :            :   }
    1252                 :            : 
    1253                 :       1379 :   e = fpr/2; /* 1/2 required prec (in sig. bits) */
    1254                 :       1379 :   eps = real2n(-e, l);
    1255                 :       1379 :   a = eval(E, gsub(x, eps), newprec);
    1256                 :       1379 :   b = eval(E, gadd(x, eps), newprec);
    1257                 :       1379 :   y = gmul2n(gsub(b,a), e-1);
    1258                 :       1379 :   return gerepileupto(av, gprec_w(y, nbits2prec(fpr)));
    1259                 :            : }
    1260                 :            : 
    1261                 :            : GEN
    1262                 :       1022 : derivfun(void *E, GEN (*eval)(void *, GEN, long), GEN x, long prec)
    1263                 :            : {
    1264                 :       1022 :   pari_sp av = avma;
    1265                 :            :   long vx;
    1266   [ +  +  +  + ]:       1022 :   switch(typ(x))
    1267                 :            :   {
    1268                 :            :   case t_REAL: case t_INT: case t_FRAC: case t_COMPLEX:
    1269                 :       1001 :     return derivnum(E,eval, x, prec);
    1270                 :            :   case t_POL:
    1271                 :          7 :     x = RgX_to_ser(x, precdl+2+1); /* +1 because deriv reduce the precision by 1 */
    1272                 :            :   case t_SER: /* FALL THROUGH */
    1273                 :         14 :     vx = varn(x);
    1274                 :         14 :     return gerepileupto(av, gdiv(deriv(eval(E, x, prec),vx), deriv(x,vx)));
    1275                 :          7 :   default: pari_err_TYPE("formal derivation",x);
    1276                 :       1015 :     return NULL; /*NOT REACHED*/
    1277                 :            :   }
    1278                 :            : }
    1279                 :            : 
    1280                 :            : GEN
    1281                 :        966 : derivnum0(GEN a, GEN code, long prec)
    1282                 :            : {
    1283                 :        966 :   EXPR_WRAP(code, derivfun (EXPR_ARGPREC,a,prec));
    1284                 :            : }
    1285                 :            : 
    1286                 :            : GEN
    1287                 :         56 : derivfun0(GEN code, GEN args, long prec)
    1288                 :            : {
    1289                 :            :   struct deriv_data E;
    1290                 :         56 :   E.code=code; E.args=args;
    1291                 :         56 :   return derivfun((void*)&E, deriv_eval, gel(args,1), prec);
    1292                 :            : }
    1293                 :            : 
    1294                 :            : /********************************************************************/
    1295                 :            : /**                   Numerical extrapolation                      **/
    1296                 :            : /********************************************************************/
    1297                 :            : 
    1298                 :            : static double
    1299                 :         91 : extgetmf(long muli)
    1300                 :            : {
    1301                 :         91 :   double mulfact[] = {0.5,0.5,0.48,0.43,0.41,0.39,0.38,0.37,0.36,0.36,0.35};
    1302         [ -  + ]:         91 :   if (muli > 100) return 0.35*LOG10_2;
    1303                 :         91 :   return mulfact[muli/10]*LOG10_2;
    1304                 :            : }
    1305                 :            : 
    1306                 :            : /* [u(n*muli), u <= N], muli = 1 unless f!=NULL */
    1307                 :            : static GEN
    1308                 :         91 : get_u(void *E, GEN (*f)(void *, GEN, long), long N, long muli, long prec)
    1309                 :            : {
    1310                 :            :   long n;
    1311                 :            :   GEN u;
    1312         [ +  + ]:         91 :   if (f)
    1313                 :            :   {
    1314                 :         77 :     u = cgetg(N+1, t_VEC);
    1315         [ +  + ]:       3367 :     for (n = 1; n <= N; n++) gel(u,n) = f(E, stoi(muli*n), prec);
    1316                 :            :   }
    1317                 :            :   else
    1318                 :            :   {
    1319                 :         14 :     u = (GEN)E;
    1320                 :         14 :     n = lg(u)-1;
    1321         [ -  + ]:         14 :     if (n < N) pari_err_COMPONENT("limitnum","<",stoi(N), stoi(n));
    1322                 :         14 :     u = vecslice(u, 1, N);
    1323                 :            :   }
    1324         [ +  + ]:       4193 :   for (n = 1; n <= N; n++)
    1325                 :            :   {
    1326                 :       4102 :     GEN un = gel(u,n);
    1327         [ +  + ]:       4102 :     if (is_rational_t(typ(un))) gel(u,n) = gtofp(un, prec);
    1328                 :            :   }
    1329                 :         91 :   return u;
    1330                 :            : }
    1331                 :            : 
    1332                 :            : struct limit
    1333                 :            : {
    1334                 :            :   long prec0; /* target accuracy */
    1335                 :            :   long prec; /* working accuracy */
    1336                 :            :   long N; /* number of terms */
    1337                 :            :   GEN u; /* sequence to extrapolate */
    1338                 :            :   GEN na; /* [n^alpha, n <= N] */
    1339                 :            :   GEN nma; /* [n^-alpha, n <= N] or NULL (alpha = 1) */
    1340                 :            :   GEN coef; /* or NULL (alpha != 1) */
    1341                 :            : };
    1342                 :            : 
    1343                 :            : static void
    1344                 :         91 : limit_init(struct limit *L, void *E, GEN (*f)(void*,GEN,long),
    1345                 :            :            long muli, GEN alpha, long prec)
    1346                 :            : {
    1347                 :         91 :   long bitprec = prec2nbits(prec), N;
    1348                 :            :   GEN na;
    1349                 :            :   long n;
    1350                 :            : 
    1351         [ +  + ]:         91 :   if (muli <= 0) muli = 20;
    1352         [ +  + ]:         91 :   if (!f) muli = 1;
    1353                 :         91 :   L->N = N = (long)ceil(extgetmf(muli)*bitprec);
    1354                 :         91 :   L->prec = nbits2prec((long)ceil(1.25*bitprec) + 32);
    1355                 :         91 :   L->prec0 = prec;
    1356                 :         91 :   L->u = get_u(E, f, N, muli, L->prec);
    1357 [ +  + ][ -  + ]:         91 :   if (alpha && gequal1(alpha)) alpha = NULL;
    1358                 :         91 :   L->na = na  = cgetg(N+1, t_VEC);
    1359         [ +  + ]:       4193 :   for (n = 1; n <= N; n++)
    1360                 :            :   {
    1361                 :       4102 :     GEN c = utoipos(n*muli);
    1362         [ +  + ]:       4102 :     if (alpha) c = gpow(c, alpha, L->prec);
    1363                 :       4102 :     gel(na,n) = c;
    1364                 :            :   }
    1365         [ +  + ]:         91 :   if (alpha)
    1366                 :            :   {
    1367                 :         28 :     GEN nma, malpha = gneg(alpha);
    1368                 :         28 :     L->coef = NULL;
    1369                 :         28 :     L->nma= nma = cgetg(N+1, t_VEC);
    1370         [ +  + ]:       1610 :     for (n = 1; n <= N; n++)
    1371                 :            :     {
    1372                 :       1582 :       GEN c = gpow(utoipos(n),malpha,L->prec);
    1373         [ +  + ]:       1582 :       if (typ(c) != t_REAL) c = gtofp(c, L->prec);
    1374                 :       1582 :       gel(nma, n) = c;
    1375                 :            :     }
    1376                 :            :   }
    1377                 :            :   else
    1378                 :            :   {
    1379                 :         63 :     GEN coef, C = vecbinome(N);
    1380                 :         63 :     L->coef = coef = cgetg(N+1, t_VEC);
    1381                 :         63 :     L->nma = NULL;
    1382         [ +  + ]:       2583 :     for (n = 1; n <= N; n++)
    1383                 :            :     {
    1384                 :       2520 :       GEN c = mulii(gel(C,n+1), powuu(n, N));
    1385         [ +  + ]:       2520 :       if (odd(N-n)) togglesign_safe(&c);
    1386                 :       2520 :       gel(coef, n) = c;
    1387                 :            :     }
    1388                 :            :   }
    1389                 :         91 : }
    1390                 :            : 
    1391                 :            : /* Zagier/Lagrange extrapolation */
    1392                 :            : static GEN
    1393                 :        676 : limitnum_i(struct limit *L)
    1394                 :            : {
    1395                 :        676 :   pari_sp av = avma;
    1396                 :            :   GEN S;
    1397         [ +  + ]:        676 :   if (L->nma)
    1398                 :        284 :     S = polint(L->nma, L->u,gen_0,NULL);
    1399                 :            :   else
    1400                 :        392 :     S = gdiv(RgV_dotproduct(L->u,L->coef), mpfact(L->N));
    1401                 :        676 :   return gerepilecopy(av, gprec_w(S, L->prec0));
    1402                 :            : }
    1403                 :            : GEN
    1404                 :         49 : limitnum(void *E, GEN (*f)(void *, GEN, long), long muli, GEN alpha, long prec)
    1405                 :            : {
    1406                 :            :   struct limit L;
    1407                 :         49 :   limit_init(&L, E,f, muli, alpha, prec);
    1408                 :         49 :   return limitnum_i(&L);
    1409                 :            : }
    1410                 :            : GEN
    1411                 :         49 : limitnum0(GEN u, long muli, GEN alpha, long prec)
    1412                 :            : {
    1413                 :         49 :   void *E = (void*)u;
    1414                 :         49 :   GEN (*f)(void*,GEN,long) = NULL;
    1415      [ +  +  - ]:         49 :   switch(typ(u))
    1416                 :            :   {
    1417                 :            :     case t_COL:
    1418                 :          7 :     case t_VEC: break;
    1419                 :         42 :     case t_CLOSURE: f = gp_callprec; break;
    1420                 :          0 :     default: pari_err_TYPE("limitnum", u);
    1421                 :            :   }
    1422                 :         49 :   return limitnum(E,f, muli,alpha, prec);
    1423                 :            : }
    1424                 :            : 
    1425                 :            : GEN
    1426                 :         42 : asympnum(void *E, GEN (*f)(void *, GEN, long), long muli, GEN alpha, long prec)
    1427                 :            : {
    1428                 :         42 :   const long MAX = 100;
    1429                 :         42 :   pari_sp av = avma;
    1430                 :         42 :   GEN u, vres = vectrunc_init(MAX);
    1431                 :            :   long i;
    1432                 :            :   struct limit L;
    1433                 :         42 :   limit_init(&L, E,f, muli, alpha, prec);
    1434                 :         42 :   u = L.u;
    1435         [ +  - ]:        627 :   for(i = 1; i <= MAX; i++)
    1436                 :            :   {
    1437                 :            :     GEN a, s, v, p, q;
    1438                 :            :     long n;
    1439                 :        627 :     s = limitnum_i(&L);
    1440                 :            :     /* NOT bestappr: lindep will "properly" ignore the lower bits */
    1441                 :        627 :     v = lindep(mkvec2(gen_1, s));
    1442                 :        627 :     p = negi(gel(v,1));
    1443                 :        627 :     q = gel(v,2);
    1444         [ -  + ]:        627 :     if (!signe(q)) break;
    1445                 :        627 :     a = gdiv(p,q);
    1446                 :        627 :     s = gsub(s, a);
    1447                 :            :     /* |s|q^2 > eps */
    1448 [ +  + ][ +  + ]:        627 :     if (!gequal0(s) && gexpo(s) + 2*expi(q) > -17) break;
    1449                 :        585 :     vectrunc_append(vres, a);
    1450         [ +  + ]:      29684 :     for (n = 1; n <= L.N; n++) gel(u,n) = gmul(gsub(gel(u,n), a), gel(L.na,n));
    1451                 :            :   }
    1452                 :         42 :   return gerepilecopy(av, vres);
    1453                 :            : }
    1454                 :            : GEN
    1455                 :         42 : asympnum0(GEN u, long muli, GEN alpha, long prec)
    1456                 :            : {
    1457                 :         42 :   void *E = (void*)u;
    1458                 :         42 :   GEN (*f)(void*,GEN,long) = NULL;
    1459      [ +  +  - ]:         42 :   switch(typ(u))
    1460                 :            :   {
    1461                 :            :     case t_COL:
    1462                 :          7 :     case t_VEC: break;
    1463                 :         35 :     case t_CLOSURE: f = gp_callprec; break;
    1464                 :          0 :     default: pari_err_TYPE("asympnum", u);
    1465                 :            :   }
    1466                 :         42 :   return asympnum(E,f, muli,alpha, prec);
    1467                 :            : }

Generated by: LCOV version 1.9