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-bordeaux1.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 16375-9f41ae0) Lines: 801 995 80.5 %
Date: 2014-04-19 Functions: 63 70 90.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 486 766 63.4 %

           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                 :    1231706 : iferrpari(GEN a, GEN b, GEN c)
      20                 :            : {
      21                 :            :   GEN res;
      22                 :            :   struct pari_evalstate state;
      23                 :    1231706 :   evalstate_save(&state);
      24         [ +  + ]:    1231706 :   pari_CATCH(CATCH_ALL)
      25                 :            :   {
      26                 :            :     GEN E;
      27 [ -  + ][ #  # ]:      14679 :     if (!b&&!c) return gnil;
      28                 :      14679 :     evalstate_restore(&state);
      29                 :      14679 :     E = gerepilecopy(avma, pari_err_last());
      30         [ +  + ]:      14679 :     if (c)
      31                 :            :     {
      32                 :         63 :       push_lex(E,c);
      33                 :         63 :       res = closure_evalgen(c);
      34                 :         63 :       pop_lex(1);
      35         [ -  + ]:         63 :       if (gequal0(res))
      36                 :          0 :         pari_err(0, E);
      37                 :            :     }
      38         [ -  + ]:      14679 :     if (!b) return gnil;
      39                 :      14679 :     push_lex(E,b);
      40                 :      14679 :     res = closure_evalgen(b);
      41                 :      14679 :     pop_lex(1);
      42                 :      14679 :     return res;
      43                 :            :   } pari_TRY {
      44                 :    1231706 :     res = closure_evalgen(a);
      45                 :    1217027 :   } pari_ENDCATCH;
      46                 :    1231706 :   return res;
      47                 :            : }
      48                 :            : 
      49                 :            : /********************************************************************/
      50                 :            : /**                                                                **/
      51                 :            : /**                        ITERATIONS                              **/
      52                 :            : /**                                                                **/
      53                 :            : /********************************************************************/
      54                 :            : 
      55                 :            : static void
      56                 :    2480814 : forparii(GEN a, GEN b, GEN code)
      57                 :            : {
      58                 :    2480814 :   pari_sp av, av0 = avma, lim;
      59                 :            :   GEN aa;
      60         [ +  + ]:    4945738 :   if (gcmp(b,a) < 0) return;
      61                 :    2464959 :   b = gfloor(b);
      62                 :    2464959 :   aa = a = setloop(a);
      63                 :    2464959 :   av=avma; lim = stack_lim(av,1);
      64                 :    2464959 :   push_lex(a,code);
      65         [ +  + ]:   24444841 :   while (gcmp(a,b) <= 0)
      66                 :            :   {
      67         [ +  + ]:   22050064 :     closure_evalvoid(code); if (loop_break()) break;
      68                 :   21979882 :     a = get_lex(-1);
      69         [ +  - ]:   21979882 :     if (a == aa)
      70                 :            :     {
      71                 :   21979882 :       a = incloop(a);
      72         [ +  + ]:   21979882 :       if (a != aa) { set_lex(-1,a); aa = a; }
      73                 :            :     }
      74                 :            :     else
      75                 :            :     { /* 'code' modified a ! Be careful (and slow) from now on */
      76                 :          0 :       a = gaddgs(a,1);
      77         [ #  # ]:          0 :       if (low_stack(lim, stack_lim(av,1)))
      78                 :            :       {
      79         [ #  # ]:          0 :         if (DEBUGMEM>1) pari_warn(warnmem,"forparii");
      80                 :          0 :         a = gerepileupto(av,a);
      81                 :            :       }
      82                 :          0 :       set_lex(-1,a);
      83                 :            :     }
      84                 :            :   }
      85                 :    2464924 :   pop_lex(1);  avma = av0;
      86                 :            : }
      87                 :            : 
      88                 :            : void
      89                 :    2480814 : forpari(GEN a, GEN b, GEN code)
      90                 :            : {
      91                 :    2480814 :   pari_sp ltop=avma, av, lim;
      92         [ +  - ]:    2480814 :   if (typ(a) == t_INT) { forparii(a,b,code); return; }
      93                 :          0 :   b = gcopy(b); /* Kludge to work-around the a+(a=2) bug */
      94                 :          0 :   av=avma; lim = stack_lim(av,1);
      95                 :          0 :   push_lex(a,code);
      96         [ #  # ]:          0 :   while (gcmp(a,b) <= 0)
      97                 :            :   {
      98         [ #  # ]:          0 :     closure_evalvoid(code); if (loop_break()) break;
      99                 :          0 :     a = get_lex(-1); a = gaddgs(a,1);
     100         [ #  # ]:          0 :     if (low_stack(lim, stack_lim(av,1)))
     101                 :            :     {
     102         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"forpari");
     103                 :          0 :       a = gerepileupto(av,a);
     104                 :            :     }
     105                 :          0 :     set_lex(-1, a);
     106                 :            :   }
     107                 :          0 :   pop_lex(1); avma = ltop;
     108                 :            : }
     109                 :            : 
     110                 :            : void
     111                 :       1197 : whilepari(GEN a, GEN b)
     112                 :            : {
     113                 :       1197 :   pari_sp av = avma;
     114                 :            :   for(;;)
     115                 :            :   {
     116                 :       8939 :     GEN res = closure_evalnobrk(a);
     117         [ +  + ]:       8939 :     if (gequal0(res)) break;
     118                 :       7742 :     avma = av;
     119         [ -  + ]:       7742 :     closure_evalvoid(b); if (loop_break()) break;
     120                 :       7742 :   }
     121                 :       1197 :   avma = av;
     122                 :       1197 : }
     123                 :            : 
     124                 :            : void
     125                 :     222075 : untilpari(GEN a, GEN b)
     126                 :            : {
     127                 :     222075 :   pari_sp av = avma;
     128                 :            :   for(;;)
     129                 :            :   {
     130                 :            :     GEN res;
     131         [ -  + ]:    2112845 :     closure_evalvoid(b); if (loop_break()) break;
     132                 :    2112845 :     res = closure_evalnobrk(a);
     133         [ +  + ]:    2112845 :     if (!gequal0(res)) break;
     134                 :    1890770 :     avma = av;
     135                 :    1890770 :   }
     136                 :     222075 :   avma = av;
     137                 :     222075 : }
     138                 :            : 
     139                 :          0 : static int negcmp(GEN x, GEN y) { return gcmp(y,x); }
     140                 :            : 
     141                 :            : void
     142                 :       1449 : forstep(GEN a, GEN b, GEN s, GEN code)
     143                 :            : {
     144                 :            :   long ss, i;
     145                 :       1449 :   pari_sp av, av0 = avma, lim;
     146                 :       1449 :   GEN v = NULL;
     147                 :            :   int (*cmp)(GEN,GEN);
     148                 :            : 
     149                 :       1449 :   b = gcopy(b); av=avma; lim = stack_lim(av,1);
     150                 :       1449 :   push_lex(a,code);
     151         [ -  + ]:       1449 :   if (is_vec_t(typ(s)))
     152                 :            :   {
     153                 :          0 :     v = s; s = gen_0;
     154         [ #  # ]:          0 :     for (i=lg(v)-1; i; i--) s = gadd(s,gel(v,i));
     155                 :            :   }
     156                 :       1449 :   ss = gsigne(s);
     157         [ +  + ]:       1449 :   if (!ss) pari_err_DOMAIN("forstep","step","=",gen_0,s);
     158         [ +  - ]:       1442 :   cmp = (ss > 0)? &gcmp: &negcmp;
     159                 :       1442 :   i = 0;
     160         [ +  + ]:      48741 :   while (cmp(a,b) <= 0)
     161                 :            :   {
     162         [ -  + ]:      47299 :     closure_evalvoid(code); if (loop_break()) break;
     163         [ -  + ]:      47299 :     if (v)
     164                 :            :     {
     165         [ #  # ]:          0 :       if (++i >= lg(v)) i = 1;
     166                 :          0 :       s = gel(v,i);
     167                 :            :     }
     168                 :      47299 :     a = get_lex(-1); a = gadd(a,s);
     169                 :            : 
     170         [ -  + ]:      47299 :     if (low_stack(lim, stack_lim(av,1)))
     171                 :            :     {
     172         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"forstep");
     173                 :          0 :       a = gerepileupto(av,a);
     174                 :            :     }
     175                 :      47299 :     set_lex(-1,a);
     176                 :            :   }
     177                 :       1442 :   pop_lex(1); avma = av0;
     178                 :       1442 : }
     179                 :            : 
     180                 :            : /* return good chunk size for sieve, 16 | chunk + 2 */
     181                 :            : static ulong
     182                 :    5174407 : optimize_chunk(ulong a, ulong b)
     183                 :            : {
     184                 :            :   /* TODO: Optimize size (surely < 512k to stay in L1 cache, but not so large */
     185                 :            :   /* as to force recalculating too often). */
     186                 :            :   /* Guesstimate: greater of sqrt(n) * lg(n) or 1M */
     187                 :    5174407 :   ulong chunk = maxuu(0x100000, usqrt(b) * expu(b));
     188                 :    5174407 :   ulong tmp = (b - a) / chunk + 1;
     189                 :            : 
     190         [ +  + ]:    5174407 :   if (tmp == 1)
     191                 :         56 :     chunk = b - a + 16;
     192                 :            :   else
     193                 :    5174351 :     chunk = (b - a) / tmp + 15;
     194                 :            :   /* Don't take up more than 2/3 of the stack */
     195                 :    5174407 :   chunk = minuu(chunk, avma - stack_lim(avma, 2));
     196                 :            :   /* ensure 16 | chunk + 2 */
     197                 :    5174407 :   return (((chunk + 2)>>4)<<4) - 2;
     198                 :            : }
     199                 :            : static void
     200                 :    5174407 : sieve_init(forprime_t *T, ulong a, ulong b)
     201                 :            : {
     202                 :    5174407 :   T->sieveb = b;
     203                 :    5174407 :   T->chunk = optimize_chunk(a, b);
     204                 :    5174407 :   T->sieve = (unsigned char*)stack_malloc(((T->chunk+2) >> 4) + 1);
     205                 :    5174407 :   T->cache[0] = 0;
     206                 :            :   /* >> 1 [only odds] + 3 [convert from bits to bytes] */
     207                 :    5174407 :   T->a = a;
     208                 :    5174407 :   T->end = minuu(a + T->chunk, b);
     209                 :    5174407 :   T->pos = T->maxpos = 0;
     210                 :    5174407 : }
     211                 :            : 
     212                 :            : static void
     213                 :    6984611 : u_forprime_set_prime_table(forprime_t *T, ulong a)
     214                 :            : {
     215                 :    6984611 :   T->strategy = 1;
     216         [ +  + ]:    6984611 :   if (a < 3)
     217                 :            :   {
     218                 :     821498 :     T->p = 0;
     219                 :     821498 :     T->d = diffptr;
     220                 :            :   }
     221                 :            :   else
     222                 :    6163113 :     T->p = init_primepointer_lt(a, &T->d);
     223                 :    6984611 : }
     224                 :            : 
     225                 :            : /* Set p so that p + q the smallest integer = c (mod q) and > original p.
     226                 :            :  * Assume 0 < c < q. Set p = 0 on overflow */
     227                 :            : static void
     228                 :        504 : arith_set(forprime_t *T)
     229                 :            : {
     230                 :        504 :   ulong r = T->p % T->q; /* 0 <= r <= min(p, q-1) */
     231                 :        504 :   pari_sp av = avma;
     232                 :        504 :   GEN d = adduu(T->p - r, T->c);
     233         [ -  + ]:        504 :   if (T->c > r) d = subiu(d, T->q);
     234                 :            :   /* d = c mod q,  d = c > r? p-r+c-q: p-r+c, so that
     235                 :            :    *  d <= p  and  d+q = c>r? p-r+c  : p-r+c+q > p */
     236                 :        504 :   T->p = itou_or_0(d); avma = av; /* d = 0 is impossible */
     237                 :        504 : }
     238                 :            : 
     239                 :            : /* run through primes in arithmetic progression = c (mod q).
     240                 :            :  * Assume (c,q)=1, 0 <= c < q */
     241                 :            : int
     242                 :    7833097 : u_forprime_arith_init(forprime_t *T, ulong a, ulong b, ulong c, ulong q)
     243                 :            : {
     244                 :            :   ulong maxp, maxp2;
     245 [ +  + ][ -  + ]:    7833097 :   if (a > b || b < 2)
     246                 :            :   {
     247                 :     847513 :     T->strategy = 1; /* paranoia */
     248                 :     847513 :     T->p = 0; /* empty */
     249                 :     847513 :     T->b = 0; /* empty */
     250                 :     847513 :     T->d = diffptr;
     251                 :     847513 :     return 0;
     252                 :            :   }
     253                 :    6985584 :   maxp = maxprime();
     254 [ +  + ][ +  - ]:    6985584 :   if (q != 1 && c != 2 && odd(q)) {
                 [ +  + ]
     255                 :            :     /* only allow *odd* primes. If c = 2, then p = 2 must be included :-( */
     256         [ -  + ]:      30226 :     if (!odd(c)) c += q;
     257                 :      30226 :     q <<= 1;
     258                 :            :   }
     259                 :    6985584 :   T->q = q;
     260                 :    6985584 :   T->c = c;
     261                 :    6985584 :   T->strategy = 0; /* unknown */
     262                 :    6985584 :   T->sieve = NULL; /* unused for now */
     263                 :    6985584 :   T->b = b;
     264         [ +  + ]:    6985584 :   if (maxp >= b) { /* [a,b] \subset prime table */
     265                 :    1780048 :     u_forprime_set_prime_table(T, a);
     266                 :    1780048 :     return 1;
     267                 :            :   }
     268                 :            :   /* b > maxp */
     269         [ +  + ]:    5205536 :   if (a >= maxp)
     270                 :            :   {
     271         [ +  - ]:        973 :     if (T->q == 1)
     272                 :        973 :       T->p = a - 1;
     273                 :            :     else
     274                 :          0 :       arith_set(T);
     275                 :            :   }
     276                 :            :   else
     277                 :    5204563 :     u_forprime_set_prime_table(T, a);
     278                 :            : 
     279         [ +  + ]:    5205536 :   maxp2 = (maxp & HIGHMASK)? 0 : maxp*maxp;
     280                 :            :   /* FIXME: should sieve as well if q != 1, adapt sieve code */
     281 [ +  + ][ +  + ]:    5205536 :   if (q != 1 || (maxp2 && maxp2 <= a)
                 [ +  + ]
     282         [ +  + ]:    5175228 :              || T->b - maxuu(a,maxp) < maxp / expu(b))
     283         [ +  + ]:      31129 :   { if (!T->strategy) T->strategy = 3; }
     284                 :            :   else
     285                 :            :   { /* worth sieving */
     286                 :            : #ifdef LONG_IS_64BIT
     287                 :    4438734 :     const ulong UPRIME_MAX = 18446744073709551557UL;
     288                 :            : #else
     289                 :     735673 :     const ulong UPRIME_MAX = 4294967291UL;
     290                 :            : #endif
     291                 :            :     ulong sieveb;
     292         [ +  + ]:    5174407 :     if (b > UPRIME_MAX) b = UPRIME_MAX;
     293                 :    5174407 :     sieveb = b;
     294 [ +  + ][ +  + ]:    5174407 :     if (maxp2 && maxp2 < b) sieveb = maxp2;
     295         [ +  + ]:    5174407 :     if (!T->strategy) T->strategy = 2;
     296         [ +  + ]:    5174407 :     if (!odd(sieveb)) sieveb--;
     297                 :    5174407 :     sieve_init(T, maxuu(maxp+2, a), sieveb);
     298                 :            :   }
     299                 :    7833097 :   return 1;
     300                 :            : }
     301                 :            : 
     302                 :            : /* will run through primes in [a,b] */
     303                 :            : int
     304                 :    7802801 : u_forprime_init(forprime_t *T, ulong a, ulong b)
     305                 :    7802801 : { return u_forprime_arith_init(T, a,b, 0,1); }
     306                 :            : /* now only run through primes <= c; assume c <= b above */
     307                 :            : void
     308                 :         42 : u_forprime_restrict(forprime_t *T, ulong c) { T->b = c; }
     309                 :            : 
     310                 :            : /* b = NULL: loop forever */
     311                 :            : int
     312                 :        569 : forprime_init(forprime_t *T, GEN a, GEN b)
     313                 :            : {
     314                 :            :   long lb;
     315         [ -  + ]:        569 :   a = gceil(a); if (typ(a) != t_INT) pari_err_TYPE("forprime_init",a);
     316         [ -  + ]:        569 :   if (signe(a) <= 0) a = gen_1;
     317         [ +  + ]:        569 :   if (b)
     318                 :            :   {
     319                 :        492 :     b = gfloor(b);
     320         [ -  + ]:        492 :     if (typ(b) != t_INT) pari_err_TYPE("forprime_init",b);
     321 [ +  - ][ -  + ]:        492 :     if (signe(b) < 0 || cmpii(a,b) > 0)
     322                 :            :     {
     323                 :          0 :       T->strategy = 4; /* paranoia */
     324                 :          0 :       T->bb = gen_0;
     325                 :          0 :       T->pp = gen_0;
     326                 :          0 :       return 0;
     327                 :            :     }
     328                 :        492 :     lb = lgefint(b);
     329                 :            :   }
     330                 :            :   else
     331                 :         77 :     lb = lgefint(a) + 4;
     332                 :        569 :   T->bb = b;
     333                 :        569 :   T->pp = cgeti(lb);
     334                 :            :   /* a, b are positive integers, a <= b */
     335         [ +  + ]:        569 :   if (lgefint(a) == 3) /* lb == 3 implies b != NULL */
     336         [ +  + ]:        560 :     return u_forprime_init(T, a[2], lb == 3? (ulong)b[2]: ULONG_MAX);
     337                 :          9 :   T->strategy = 4;
     338                 :          9 :   affii(subiu(a,1), T->pp);
     339                 :        569 :   return 1;
     340                 :            : }
     341                 :            : 
     342                 :            : /* assume a <= b <= maxprime()^2, a,b odd, sieve[n] corresponds to
     343                 :            :  *   a+16*n, a+16*n+2, ..., a+16*n+14 (bits 0 to 7)
     344                 :            :  * maxpos = index of last sieve cell. */
     345                 :            : static void
     346                 :        903 : sieve_block(ulong a, ulong b, ulong maxpos, unsigned char* sieve)
     347                 :            : {
     348                 :        903 :   ulong p = 2, lim = usqrt(b), sz = (b-a) >> 1;
     349                 :        903 :   byteptr d = diffptr+1;
     350                 :        903 :   (void)memset(sieve, 0, maxpos+1);
     351                 :            :   for (;;)
     352                 :            :   { /* p is odd */
     353                 :            :     ulong k, r;
     354                 :    2154673 :     NEXT_PRIME_VIADIFF(p, d); /* starts at p = 3 */
     355         [ +  + ]:    2154673 :     if (p > lim) break;
     356                 :            : 
     357                 :            :     /* solve a + 2k = 0 (mod p) */
     358                 :    2153770 :     r = a % p;
     359         [ +  + ]:    2153770 :     if (r == 0)
     360                 :       1785 :       k = 0;
     361                 :            :     else
     362                 :            :     {
     363                 :    2151985 :       k = p - r;
     364         [ +  + ]:    2151985 :       if (odd(k)) k += p;
     365                 :    2151985 :       k >>= 1;
     366                 :            :     }
     367                 :            :     /* m = a + 2k is the smallest odd m >= a, p | m */
     368                 :            :     /* position n (corresponds to a+2n) is sieve[n>>3], bit n&7 */
     369         [ +  + ]: 1445494720 :     while (k <= sz) { sieve[k>>3] |= 1 << (k&7); k += p; /* 2k += 2p */ }
     370                 :    2153770 :   }
     371                 :        903 : }
     372                 :            : 
     373                 :            : /* T->cache is a 0-terminated list of primes, return the first one and
     374                 :            :  * remove it from list. Most of the time the list contains a single prime */
     375                 :            : static ulong
     376                 :   57544767 : shift_cache(forprime_t *T)
     377                 :            : {
     378                 :            :   long i;
     379                 :   57544767 :   T->p = T->cache[0];
     380                 :   57544767 :   for (i = 1;; i++)  /* remove one prime from cache */
     381         [ +  + ]:   75756709 :     if (! (T->cache[i-1] = T->cache[i]) ) break;
     382                 :   75756709 :   return T->p;
     383                 :            : }
     384                 :            : 
     385                 :            : ulong
     386                 :   78042325 : u_forprime_next(forprime_t *T)
     387                 :            : {
     388         [ +  + ]:   78042325 :   if (T->strategy == 1)
     389                 :            :   {
     390                 :            :     for(;;)
     391                 :            :     {
     392         [ +  + ]:  161958112 :       if (!*(T->d))
     393                 :            :       {
     394         [ +  + ]:       5378 :         T->strategy = T->sieve? 2: 3;
     395 [ +  + ][ -  + ]:       5378 :         if (T->q != 1) { arith_set(T); if (!T->p) return 0; }
     396                 :            :         /* T->p possibly not a prime if q != 1 */
     397                 :       5378 :         break;
     398                 :            :       }
     399                 :            :       else
     400                 :            :       {
     401                 :  161952734 :         NEXT_PRIME_VIADIFF(T->p, T->d);
     402         [ +  + ]:  161952734 :         if (T->p > T->b) return 0;
     403 [ +  + ][ +  + ]:  161912267 :         if (T->q == 1 || T->p % T->q == T->c) return T->p;
     404                 :            :       }
     405                 :  143936541 :     }
     406                 :            :   }
     407         [ +  + ]:   60026132 :   if (T->strategy == 2)
     408                 :            :   {
     409                 :            :     ulong n;
     410         [ +  + ]:   57544893 :     if (T->cache[0]) return shift_cache(T);
     411                 :            : NEXT_CHUNK:
     412         [ +  + ]:   69110013 :     for (n = T->pos; n < T->maxpos; n++)
     413         [ +  + ]:   69108431 :       if (T->sieve[n] != 0xFF)
     414                 :            :       {
     415                 :   41857627 :         unsigned char mask = T->sieve[n];
     416                 :   41857627 :         ulong p = T->a + (n<<4);
     417                 :   41857627 :         long i = 0;
     418                 :   41857627 :         T->pos = n;
     419         [ +  + ]:   41857627 :         if (!(mask &  1)) T->cache[i++] = p;
     420         [ +  + ]:   41857627 :         if (!(mask &  2)) T->cache[i++] = p+2;
     421         [ +  + ]:   41857627 :         if (!(mask &  4)) T->cache[i++] = p+4;
     422         [ +  + ]:   41857627 :         if (!(mask &  8)) T->cache[i++] = p+6;
     423         [ +  + ]:   41857627 :         if (!(mask & 16)) T->cache[i++] = p+8;
     424         [ +  + ]:   41857627 :         if (!(mask & 32)) T->cache[i++] = p+10;
     425         [ +  + ]:   41857627 :         if (!(mask & 64)) T->cache[i++] = p+12;
     426         [ +  + ]:   41857627 :         if (!(mask &128)) T->cache[i++] = p+14;
     427                 :   41857627 :         T->cache[i] = 0;
     428                 :   41857627 :         T->pos = n+1;
     429                 :   41857627 :         return shift_cache(T);
     430                 :            :       }
     431                 :            :     /* n = T->maxpos, last cell: check p <= b */
     432 [ +  + ][ +  + ]:       1582 :     if (T->maxpos && n == T->maxpos && T->sieve[n] != 0xFF)
                 [ +  + ]
     433                 :            :     {
     434                 :        581 :       unsigned char mask = T->sieve[n];
     435                 :        581 :       ulong p = T->a + (n<<4);
     436                 :        581 :       long i = 0;
     437                 :        581 :       T->pos = n;
     438 [ +  + ][ +  - ]:        581 :       if (!(mask &  1) && p <= T->sieveb) T->cache[i++] = p;
     439 [ +  + ][ +  - ]:        581 :       if (!(mask &  2) && p <= T->sieveb-2) T->cache[i++] = p+2;
     440 [ +  + ][ +  + ]:        581 :       if (!(mask &  4) && p <= T->sieveb-4) T->cache[i++] = p+4;
     441 [ +  + ][ +  + ]:        581 :       if (!(mask &  8) && p <= T->sieveb-6) T->cache[i++] = p+6;
     442 [ +  + ][ +  + ]:        581 :       if (!(mask & 16) && p <= T->sieveb-8) T->cache[i++] = p+8;
     443 [ +  + ][ +  + ]:        581 :       if (!(mask & 32) && p <= T->sieveb-10) T->cache[i++] = p+10;
     444 [ +  + ][ +  + ]:        581 :       if (!(mask & 64) && p <= T->sieveb-12) T->cache[i++] = p+12;
     445 [ +  + ][ +  + ]:        581 :       if (!(mask &128) && p <= T->sieveb-14) T->cache[i++] = p+14;
     446         [ +  + ]:        581 :       if (i)
     447                 :            :       {
     448                 :        553 :         T->cache[i] = 0;
     449                 :        553 :         T->pos = n+1;
     450                 :        553 :         return shift_cache(T);
     451                 :            :       }
     452                 :            :     }
     453                 :            : 
     454         [ +  + ]:       1029 :     if (T->end >= T->sieveb) /* done */
     455                 :        126 :       T->strategy = 3;
     456                 :            :     else
     457                 :            :     { /* initialize next chunk */
     458         [ +  + ]:        903 :       if (T->maxpos == 0)
     459                 :        140 :         T->a |= 1; /* first time; ensure odd */
     460                 :            :       else
     461                 :        763 :         T->a = (T->end + 2) | 1;
     462                 :        903 :       T->end = T->a + T->chunk; /* may overflow */
     463 [ +  + ][ +  + ]:        903 :       if (T->end < T->a || T->end > T->sieveb) T->end = T->sieveb;
     464                 :            :       /* end and a are odd; sieve[k] contains the a + 8*2k + (0,2,...,14).
     465                 :            :        * The largest k is (end-a) >> 4 */
     466                 :        903 :       T->pos = 0;
     467                 :        903 :       T->maxpos = (T->end - T->a) >> 4;
     468                 :        903 :       sieve_block(T->a, T->end, T->maxpos, T->sieve);
     469                 :        903 :       goto NEXT_CHUNK;
     470                 :            :     }
     471                 :            :   }
     472         [ +  - ]:    2481365 :   if (T->strategy == 3)
     473                 :            :   {
     474         [ +  + ]:    2481365 :     if (T->q == 1)
     475                 :    2480861 :       T->p = unextprime(T->p + 1);
     476                 :            :     else
     477                 :            :     {
     478                 :            :       do {
     479                 :       2499 :         T->p += T->q;
     480         [ -  + ]:       2499 :         if (T->p < T->q) return 0; /* overflow */
     481         [ +  + ]:       2499 :       } while (!uisprime(T->p));
     482                 :            :     }
     483         [ +  + ]:    2481365 :     if (!T->p) /* overflow ulong, switch to GEN */
     484                 :         15 :       T->strategy = 4;
     485                 :            :     else
     486                 :            :     {
     487         [ +  + ]:    2481350 :       if (T->p > T->b) return 0;
     488                 :    2475595 :       return T->p;
     489                 :            :     }
     490                 :            :   }
     491                 :   78042325 :   return 0; /* overflow */
     492                 :            : }
     493                 :            : 
     494                 :            : GEN
     495                 :    8457930 : forprime_next(forprime_t *T)
     496                 :            : {
     497                 :            :   pari_sp av;
     498                 :            :   GEN p;
     499         [ +  + ]:    8457930 :   if (T->strategy <= 3)
     500                 :            :   {
     501                 :    8457847 :     ulong q = u_forprime_next(T);
     502         [ +  + ]:    8457847 :     if (q) { affui(q, T->pp); return T->pp; }
     503                 :            :     /* failure */
     504         [ +  + ]:        469 :     if (T->strategy <= 3) return NULL; /* we're done */
     505                 :            :     /* overflow ulong, switch to GEN */
     506                 :         15 :     affui(ULONG_MAX, T->pp); /* initialize */
     507                 :            :   }
     508                 :         98 :   av = avma;
     509                 :         98 :   p = nextprime(addiu(T->pp, 1));
     510 [ +  - ][ +  + ]:         98 :   if (T->bb && absi_cmp(p, T->bb) > 0) { avma = av; return NULL; }
     511                 :    8457930 :   affii(p, T->pp); avma = av; return T->pp;
     512                 :            : }
     513                 :            : 
     514                 :            : void
     515                 :        147 : forprime(GEN a, GEN b, GEN code)
     516                 :            : {
     517                 :        147 :   pari_sp av = avma;
     518                 :            :   forprime_t T;
     519                 :            : 
     520         [ -  + ]:        294 :   if (!forprime_init(&T, a,b)) { avma = av; return; }
     521                 :            : 
     522                 :        147 :   push_lex(T.pp,code);
     523         [ +  + ]:       2254 :   while(forprime_next(&T))
     524                 :            :   {
     525         [ -  + ]:       2107 :     closure_evalvoid(code); if (loop_break()) break;
     526                 :            :     /* p changed in 'code', complain */
     527         [ -  + ]:       2107 :     if (get_lex(-1) != T.pp)
     528                 :          0 :       pari_err(e_MISC,"prime index read-only: was changed to %Ps", get_lex(-1));
     529                 :            :   }
     530                 :        147 :   pop_lex(1); avma = av;
     531                 :            : }
     532                 :            : 
     533                 :            : int
     534                 :         35 : forcomposite_init(forcomposite_t *C, GEN a, GEN b)
     535                 :            : {
     536                 :         35 :   pari_sp av = avma;
     537         [ -  + ]:         35 :   a = gceil(a); if (typ(a)!=t_INT) pari_err_TYPE("forcomposite",a);
     538         [ +  + ]:         35 :   if (b) {
     539         [ -  + ]:         28 :     b = gfloor(b);if (typ(b)!=t_INT) pari_err_TYPE("forcomposite",b);
     540                 :            :   }
     541         [ -  + ]:         35 :   if (signe(a) < 0) pari_err_DOMAIN("forcomposite", "a", "<", gen_0, a);
     542         [ +  + ]:         35 :   if (cmpiu(a, 4) < 0) a = utoipos(4);
     543                 :         35 :   C->first = 1;
     544         [ -  + ]:         35 :   if (!forprime_init(&C->T, a,b))
     545                 :            :   {
     546                 :          0 :     C->n = gen_1; /* in case caller forgets to check the return value */
     547                 :          0 :     C->b = gen_0;
     548                 :          0 :     avma = av; return 0;
     549                 :            :   }
     550                 :         35 :   C->n = setloop(a);
     551                 :         35 :   C->b = b;
     552                 :         35 :   C->p = NULL; return 1;
     553                 :            : }
     554                 :            : 
     555                 :            : GEN
     556                 :        147 : forcomposite_next(forcomposite_t *C)
     557                 :            : {
     558         [ +  + ]:        147 :   if (C->first) /* first call ever */
     559                 :            :   {
     560                 :         35 :     C->first = 0;
     561                 :         35 :     C->p = forprime_next(&C->T);
     562                 :            :   }
     563                 :            :   else
     564                 :        112 :     incloop(C->n);
     565         [ +  + ]:        147 :   if (C->p)
     566                 :            :   {
     567         [ +  + ]:        119 :     if (cmpii(C->n, C->p) < 0) return C->n;
     568                 :         56 :     incloop(C->n);
     569                 :            :     /* n = p+1 */
     570                 :         56 :     C->p = forprime_next(&C->T); /* nextprime(p) > n */
     571         [ +  + ]:         56 :     if (C->p) return C->n;
     572                 :            :   }
     573 [ +  - ][ +  + ]:         49 :   if (!C->b || cmpii(C->n, C->b) <= 0) return C->n;
     574                 :        147 :   return NULL;
     575                 :            : }
     576                 :            : 
     577                 :            : void
     578                 :         35 : forcomposite(GEN a, GEN b, GEN code)
     579                 :            : {
     580                 :         35 :   pari_sp av = avma;
     581                 :            :   forcomposite_t T;
     582                 :            :   GEN n;
     583         [ +  - ]:         63 :   if (!forcomposite_init(&T,a,b)) return;
     584                 :         35 :   push_lex(T.n,code);
     585         [ +  + ]:        147 :   while((n = forcomposite_next(&T)))
     586                 :            :   {
     587         [ +  + ]:        126 :     closure_evalvoid(code); if (loop_break()) break;
     588                 :            :     /* n changed in 'code', complain */
     589         [ +  + ]:        119 :     if (get_lex(-1) != n)
     590                 :          7 :       pari_err(e_MISC,"index read-only: was changed to %Ps", get_lex(-1));
     591                 :            :   }
     592                 :         28 :   pop_lex(1); avma = av;
     593                 :            : }
     594                 :            : 
     595                 :            : void
     596                 :         14 : fordiv(GEN a, GEN code)
     597                 :            : {
     598                 :            :   long i, l;
     599                 :         14 :   pari_sp av2, av = avma;
     600                 :         14 :   GEN t = divisors(a);
     601                 :            : 
     602                 :         14 :   push_lex(gen_0,code); l=lg(t); av2 = avma;
     603         [ +  + ]:         70 :   for (i=1; i<l; i++)
     604                 :            :   {
     605                 :         56 :     set_lex(-1,gel(t,i));
     606         [ -  + ]:         56 :     closure_evalvoid(code); if (loop_break()) break;
     607                 :         56 :     avma = av2;
     608                 :            :   }
     609                 :         14 :   pop_lex(1); avma=av;
     610                 :         14 : }
     611                 :            : 
     612                 :            : /* Embedded for loops:
     613                 :            :  *   fl = 0: execute ch (a), where a = (ai) runs through all n-uplets in
     614                 :            :  *     [m1,M1] x ... x [mn,Mn]
     615                 :            :  *   fl = 1: impose a1 <= ... <= an
     616                 :            :  *   fl = 2:        a1 <  ... <  an
     617                 :            :  */
     618                 :            : /* increment and return d->a [over integers]*/
     619                 :            : static GEN
     620                 :      27083 : _next_i(forvec_t *d)
     621                 :            : {
     622                 :      27083 :   long i = d->n;
     623         [ +  + ]:      27083 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     624                 :            :   for (;;) {
     625         [ +  + ]:      34993 :     if (cmpii(d->a[i], d->M[i]) < 0) {
     626                 :      26957 :       d->a[i] = incloop(d->a[i]);
     627                 :      26957 :       return (GEN)d->a;
     628                 :            :     }
     629                 :       8036 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     630         [ +  + ]:       8036 :     if (--i <= 0) return NULL;
     631                 :      35056 :   }
     632                 :            : }
     633                 :            : /* increment and return d->a [generic]*/
     634                 :            : static GEN
     635                 :          0 : _next(forvec_t *d)
     636                 :            : {
     637                 :          0 :   long i = d->n;
     638         [ #  # ]:          0 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     639                 :            :   for (;;) {
     640                 :          0 :     d->a[i] = gaddgs(d->a[i], 1);
     641         [ #  # ]:          0 :     if (gcmp(d->a[i], d->M[i]) <= 0) return (GEN)d->a;
     642                 :          0 :     d->a[i] = d->m[i];
     643         [ #  # ]:          0 :     if (--i <= 0) return NULL;
     644                 :          0 :   }
     645                 :            : }
     646                 :            : 
     647                 :            : /* non-decreasing order [over integers] */
     648                 :            : static GEN
     649                 :          0 : _next_le_i(forvec_t *d)
     650                 :            : {
     651                 :          0 :   long i = d->n;
     652         [ #  # ]:          0 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     653                 :            :   for (;;) {
     654         [ #  # ]:          0 :     if (cmpii(d->a[i], d->M[i]) < 0)
     655                 :            :     {
     656                 :          0 :       d->a[i] = incloop(d->a[i]);
     657                 :            :       /* m[i] < a[i] <= M[i] < M[i+1] */
     658         [ #  # ]:          0 :       while (i < d->n)
     659                 :            :       {
     660                 :            :         GEN t;
     661                 :          0 :         i++;
     662         [ #  # ]:          0 :         if (cmpii(d->a[i-1], d->a[i]) <= 0) continue;
     663                 :            :         /* a[i-1] <= M[i-1] <= M[i] */
     664         [ #  # ]:          0 :         t = d->a[i-1]; if (cmpii(t, d->m[i]) < 0) t = d->m[i];
     665                 :          0 :         d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1],m[i])*/
     666                 :            :       }
     667                 :          0 :       return (GEN)d->a;
     668                 :            :     }
     669                 :          0 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     670         [ #  # ]:          0 :     if (--i <= 0) return NULL;
     671                 :          0 :   }
     672                 :            : }
     673                 :            : /* non-decreasing order [generic] */
     674                 :            : static GEN
     675                 :          0 : _next_le(forvec_t *d)
     676                 :            : {
     677                 :          0 :   long i = d->n, imin = d->n;
     678         [ #  # ]:          0 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     679                 :            :   for (;;) {
     680                 :          0 :     d->a[i] = gaddgs(d->a[i], 1);
     681         [ #  # ]:          0 :     if (gcmp(d->a[i], d->M[i]) <= 0)
     682                 :            :     {
     683         [ #  # ]:          0 :       while (i < d->n)
     684                 :            :       {
     685                 :          0 :         i++;
     686         [ #  # ]:          0 :         if (gcmp(d->a[i-1], d->a[i]) <= 0) continue;
     687         [ #  # ]:          0 :         while (gcmp(d->a[i-1], d->M[i]) > 0)
     688                 :            :         {
     689         [ #  # ]:          0 :           i = imin - 1; if (!i) return NULL;
     690                 :          0 :           imin = i;
     691                 :          0 :           d->a[i] = gaddgs(d->a[i], 1);
     692         [ #  # ]:          0 :           if (gcmp(d->a[i], d->M[i]) <= 0) break;
     693                 :            :         }
     694         [ #  # ]:          0 :         if (i > 1) { /* a >= a[i-1] - a[i] */
     695                 :          0 :           GEN a = gceil(gsub(d->a[i-1], d->a[i]));
     696                 :          0 :           d->a[i] = gadd(d->a[i], a);
     697                 :            :         }
     698                 :            :       }
     699                 :          0 :       return (GEN)d->a;
     700                 :            :     }
     701                 :          0 :     d->a[i] = d->m[i];
     702         [ #  # ]:          0 :     if (--i <= 0) return NULL;
     703         [ #  # ]:          0 :     if (i < imin) imin = i;
     704                 :          0 :   }
     705                 :            : }
     706                 :            : /* strictly increasing order [over integers] */
     707                 :            : static GEN
     708                 :    1173466 : _next_lt_i(forvec_t *d)
     709                 :            : {
     710                 :    1173466 :   long i = d->n;
     711         [ +  + ]:    1173466 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     712                 :            :   for (;;) {
     713         [ +  + ]:    1289960 :     if (cmpii(d->a[i], d->M[i]) < 0)
     714                 :            :     {
     715                 :    1159886 :       d->a[i] = incloop(d->a[i]);
     716                 :            :       /* m[i] < a[i] <= M[i] < M[i+1] */
     717         [ +  + ]:    1276373 :       while (i < d->n)
     718                 :            :       {
     719                 :            :         pari_sp av;
     720                 :            :         GEN t;
     721                 :     116487 :         i++;
     722         [ -  + ]:     116487 :         if (cmpii(d->a[i-1], d->a[i]) < 0) continue;
     723                 :     116487 :         av = avma;
     724                 :            :         /* a[i-1] <= M[i-1] < M[i] */
     725         [ -  + ]:     116487 :         t = addis(d->a[i-1],1); if (cmpii(t, d->m[i]) < 0) t = d->m[i];
     726                 :     116487 :         d->a[i] = resetloop(d->a[i], t);/*a[i]:=max(a[i-1]+1,m[i]) <= M[i]*/
     727                 :     116487 :         avma = av;
     728                 :            :       }
     729                 :    1159886 :       return (GEN)d->a;
     730                 :            :     }
     731                 :     130074 :     d->a[i] = resetloop(d->a[i], d->m[i]);
     732         [ +  + ]:     130074 :     if (--i <= 0) return NULL;
     733                 :    1296750 :   }
     734                 :            : }
     735                 :            : /* strictly increasing order [generic] */
     736                 :            : static GEN
     737                 :          0 : _next_lt(forvec_t *d)
     738                 :            : {
     739                 :          0 :   long i = d->n, imin = d->n;
     740         [ #  # ]:          0 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     741                 :            :   for (;;) {
     742                 :          0 :     d->a[i] = gaddgs(d->a[i], 1);
     743         [ #  # ]:          0 :     if (gcmp(d->a[i], d->M[i]) <= 0)
     744                 :            :     {
     745         [ #  # ]:          0 :       while (i < d->n)
     746                 :            :       {
     747                 :          0 :         i++;
     748         [ #  # ]:          0 :         if (gcmp(d->a[i-1], d->a[i]) < 0) continue;
     749                 :            :         for(;;)
     750                 :            :         {
     751                 :            :           GEN a, b;
     752                 :          0 :           a = addis(gfloor(gsub(d->a[i-1], d->a[i])), 1); /* a> v[i-1]-v[i] */
     753                 :          0 :           b = gadd(d->a[i], a);
     754                 :            :           /* v[i-1] < b <= v[i-1] + 1 */
     755         [ #  # ]:          0 :           if (gcmp(b, d->M[i]) <= 0) { d->a[i] = b; break; }
     756                 :            : 
     757         [ #  # ]:          0 :           for (; i >= imin; i--) d->a[i] = d->m[i];
     758         [ #  # ]:          0 :           if (!i) return NULL;
     759                 :          0 :           imin = i;
     760                 :          0 :           d->a[i] = gaddgs(d->a[i], 1);
     761         [ #  # ]:          0 :           if (gcmp(d->a[i], d->M[i]) <= 0) break;
     762                 :          0 :         }
     763                 :            :       }
     764                 :          0 :       return (GEN)d->a;
     765                 :            :     }
     766                 :          0 :     d->a[i] = d->m[i];
     767         [ #  # ]:          0 :     if (--i <= 0) return NULL;
     768         [ #  # ]:          0 :     if (i < imin) imin = i;
     769                 :          0 :   }
     770                 :            : }
     771                 :            : /* for forvec(v=[],) */
     772                 :            : static GEN
     773                 :          0 : _next_void(forvec_t *d)
     774                 :            : {
     775         [ #  # ]:          0 :   if (d->first) { d->first = 0; return (GEN)d->a; }
     776                 :          0 :   return NULL;
     777                 :            : }
     778                 :            : 
     779                 :            : /* Initialize minima (m) and maxima (M); guarantee
     780                 :            :  *   if flag = 1: m[i-1] <= m[i] <= M[i] <= M[i+1]
     781                 :            :  *   if flag = 2: m[i-1] <  m[i] <= M[i] <  M[i+1] */
     782                 :            : int
     783                 :       6853 : forvec_init(forvec_t *d, GEN x, long flag)
     784                 :            : {
     785                 :       6853 :   long i, tx = typ(x), l = lg(x), t = t_INT;
     786         [ -  + ]:       6853 :   if (!is_vec_t(tx)) pari_err_TYPE("forvec [not a vector]", x);
     787                 :       6853 :   d->first = 1;
     788                 :       6853 :   d->n = l - 1;
     789                 :       6853 :   d->a = (GEN*)cgetg(l,tx);
     790                 :       6853 :   d->m = (GEN*)cgetg(l,tx);
     791                 :       6853 :   d->M = (GEN*)cgetg(l,tx);
     792         [ -  + ]:       6853 :   if (l == 1) { d->next = &_next_void; return 1; }
     793         [ +  + ]:      20678 :   for (i = 1; i < l; i++)
     794                 :            :   {
     795                 :      13825 :     GEN a, e = gel(x,i), m = gel(e,1), M = gel(e,2);
     796                 :      13825 :     tx = typ(e);
     797 [ +  - ][ -  + ]:      13825 :     if (! is_vec_t(tx) || lg(e)!=3)
     798                 :          0 :       pari_err_TYPE("forvec [expected vector not of type [min,MAX]]",e);
     799         [ -  + ]:      13825 :     if (typ(m) != t_INT) t = t_REAL;
     800         [ +  + ]:      13825 :     if (i > 1) switch(flag)
              [ -  +  + ]
     801                 :            :     {
     802                 :            :       case 1: /* a >= m[i-1] - m */
     803                 :          0 :         a = gceil(gsub(d->m[i-1], m));
     804         [ #  # ]:          0 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     805         [ #  # ]:          0 :         if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
     806                 :          0 :         break;
     807                 :            :       case 2: /* a > m[i-1] - m */
     808                 :       6797 :         a = gfloor(gsub(d->m[i-1], m));
     809         [ -  + ]:       6797 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     810                 :       6797 :         a = addis(a, 1);
     811         [ +  - ]:       6797 :         if (signe(a) > 0) m = gadd(m, a); else m = gcopy(m);
     812                 :       6797 :         break;
     813                 :        175 :       default: m = gcopy(m);
     814                 :        175 :         break;
     815                 :            :     }
     816         [ -  + ]:      13825 :     if (gcmp(m,M) > 0) { d->a = NULL; d->next = &_next; return 0; }
     817                 :      13825 :     d->m[i] = m;
     818                 :      13825 :     d->M[i] = M;
     819                 :            :   }
     820         [ +  + ]:      13825 :   for (i = l-2; i >= 1; i--)
     821                 :            :   {
     822                 :       6972 :     GEN a, M = d->M[i];
     823      [ -  +  + ]:       6972 :     switch(flag) {
     824                 :            :       case 1:/* a >= M - M[i] */
     825                 :          0 :         a = gfloor(gsub(d->M[i+1], M));
     826         [ #  # ]:          0 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     827         [ #  # ]:          0 :         if (signe(a) < 0) M = gadd(M, a); else M = gcopy(M);
     828                 :            :         /* M <= M[i+1] */
     829                 :          0 :         break;
     830                 :            :       case 2:
     831                 :       6797 :         a = gceil(gsub(d->M[i+1], M));
     832         [ -  + ]:       6797 :         if (typ(a) != t_INT) pari_err_TYPE("forvec",a);
     833                 :       6797 :         a = subis(a, 1);
     834         [ +  - ]:       6797 :         if (signe(a) < 0) M = gadd(M, a); else M = gcopy(M);
     835                 :            :         /* M < M[i+1] */
     836                 :       6797 :         break;
     837                 :            :       default:
     838                 :        175 :         M = gcopy(M);
     839                 :        175 :         break;
     840                 :            :     }
     841                 :       6972 :     d->M[i] = M;
     842                 :            :   }
     843         [ +  - ]:       6853 :   if (t == t_INT) {
     844         [ +  + ]:      20678 :     for (i = 1; i < l; i++) {
     845                 :      13825 :       d->a[i] = setloop(d->m[i]);
     846         [ -  + ]:      13825 :       if (typ(d->M[i]) != t_INT) d->M[i] = gfloor(d->M[i]);
     847                 :            :     }
     848                 :            :   } else {
     849         [ #  # ]:          0 :     for (i = 1; i < l; i++) d->a[i] = d->m[i];
     850                 :            :   }
     851   [ +  -  +  - ]:       6853 :   switch(flag)
     852                 :            :   {
     853         [ +  - ]:         63 :     case 0: d->next = t==t_INT? &_next_i:    &_next; break;
     854         [ #  # ]:          0 :     case 1: d->next = t==t_INT? &_next_le_i: &_next_le; break;
     855         [ +  - ]:       6790 :     case 2: d->next = t==t_INT? &_next_lt_i: &_next_lt; break;
     856                 :          0 :     default: pari_err_FLAG("forvec");
     857                 :            :   }
     858                 :       6853 :   return 1;
     859                 :            : }
     860                 :            : GEN
     861                 :    1200549 : forvec_next(forvec_t *d) { return d->next(d); }
     862                 :            : 
     863                 :            : void
     864                 :       6853 : forvec(GEN x, GEN code, long flag)
     865                 :            : {
     866                 :       6853 :   pari_sp av = avma;
     867                 :            :   forvec_t T;
     868                 :            :   GEN v;
     869         [ -  + ]:      13706 :   if (!forvec_init(&T, x, flag)) { avma = av; return; }
     870                 :       6853 :   push_lex((GEN)T.a, code);
     871         [ +  + ]:    1200549 :   while ((v = forvec_next(&T)))
     872                 :            :   {
     873                 :    1193696 :     closure_evalvoid(code);
     874         [ -  + ]:    1193696 :     if (loop_break()) break;
     875                 :            :   }
     876                 :       6853 :   pop_lex(1); avma = av;
     877                 :            : }
     878                 :            : 
     879                 :            : /********************************************************************/
     880                 :            : /**                                                                **/
     881                 :            : /**                              SUMS                              **/
     882                 :            : /**                                                                **/
     883                 :            : /********************************************************************/
     884                 :            : 
     885                 :            : GEN
     886                 :        980 : somme(GEN a, GEN b, GEN code, GEN x)
     887                 :            : {
     888                 :        980 :   pari_sp av, av0 = avma, lim;
     889                 :            :   GEN p1;
     890                 :            : 
     891         [ -  + ]:        980 :   if (typ(a) != t_INT) pari_err_TYPE("sum",a);
     892         [ +  + ]:        980 :   if (!x) x = gen_0;
     893         [ -  + ]:        980 :   if (gcmp(b,a) < 0) return gcopy(x);
     894                 :            : 
     895                 :        980 :   b = gfloor(b);
     896                 :        980 :   a = setloop(a);
     897                 :        980 :   av=avma; lim = stack_lim(av,1);
     898                 :        980 :   push_lex(a,code);
     899                 :            :   for(;;)
     900                 :            :   {
     901                 :     865809 :     p1 = closure_evalnobrk(code);
     902         [ +  + ]:     865809 :     x=gadd(x,p1); if (cmpii(a,b) >= 0) break;
     903                 :     864829 :     a = incloop(a);
     904         [ -  + ]:     864829 :     if (low_stack(lim, stack_lim(av,1)))
     905                 :            :     {
     906         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sum");
     907                 :          0 :       x = gerepileupto(av,x);
     908                 :            :     }
     909                 :     864829 :     set_lex(-1,a);
     910                 :     864829 :   }
     911                 :        980 :   pop_lex(1); return gerepileupto(av0,x);
     912                 :            : }
     913                 :            : 
     914                 :            : GEN
     915                 :         21 : suminf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
     916                 :            : {
     917                 :            :   long fl, G;
     918                 :         21 :   pari_sp av0 = avma, av, lim;
     919                 :         21 :   GEN p1,x = real_1(prec);
     920                 :            : 
     921         [ -  + ]:         21 :   if (typ(a) != t_INT) pari_err_TYPE("suminf",a);
     922                 :         21 :   a = setloop(a);
     923                 :         21 :   av = avma; lim = stack_lim(av,1);
     924                 :         21 :   fl=0; G = prec2nbits(prec) + 5;
     925                 :            :   for(;;)
     926                 :            :   {
     927                 :       9107 :     p1 = eval(E, a); x = gadd(x,p1); a = incloop(a);
     928 [ +  - ][ +  + ]:       9107 :     if (gequal0(p1) || gexpo(p1) <= gexpo(x)-G)
     929         [ +  + ]:         63 :       { if (++fl==3) break; }
     930                 :            :     else
     931                 :       9044 :       fl=0;
     932         [ -  + ]:       9086 :     if (low_stack(lim, stack_lim(av,1)))
     933                 :            :     {
     934         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"suminf");
     935                 :          0 :       x = gerepileupto(av,x);
     936                 :            :     }
     937                 :       9086 :   }
     938                 :         21 :   return gerepileupto(av0, gaddgs(x,-1));
     939                 :            : }
     940                 :            : GEN
     941                 :         21 : suminf0(GEN a, GEN code, long prec)
     942                 :         21 : { EXPR_WRAP(code, suminf(EXPR_ARG, a, prec)); }
     943                 :            : 
     944                 :            : GEN
     945                 :         56 : sumdivexpr(GEN num, GEN code)
     946                 :            : {
     947                 :         56 :   pari_sp av = avma;
     948                 :         56 :   GEN y = gen_0, t = divisors(num);
     949                 :         56 :   long i, l = lg(t);
     950                 :            : 
     951                 :         56 :   push_lex(gen_0, code);
     952         [ +  + ]:       9968 :   for (i=1; i<l; i++)
     953                 :            :   {
     954                 :       9912 :     set_lex(-1,gel(t,i));
     955                 :       9912 :     y = gadd(y, closure_evalnobrk(code));
     956                 :            :   }
     957                 :         56 :   pop_lex(1); return gerepileupto(av,y);
     958                 :            : }
     959                 :            : GEN
     960                 :         42 : sumdivmultexpr(GEN num, GEN code)
     961                 :            : {
     962                 :         42 :   pari_sp av = avma;
     963                 :         42 :   GEN y = gen_1, P,E;
     964                 :         42 :   int isint = divisors_init(num, &P,&E);
     965                 :         42 :   long i, l = lg(P);
     966                 :            : 
     967         [ -  + ]:         42 :   if (l == 1) { avma = av; return gen_1; }
     968                 :         42 :   push_lex(gen_0, code);
     969         [ +  + ]:        196 :   for (i=1; i<l; i++)
     970                 :            :   {
     971                 :        154 :     GEN p = gel(P,i), q = p, z = gen_1;
     972                 :        154 :     long j, e = E[i];
     973 [ +  + ][ +  - ]:        560 :     for (j = 1; j <= e; j++, q = isint?mulii(q, p): gmul(q,p))
     974                 :            :     {
     975                 :        560 :       set_lex(-1, q);
     976                 :        560 :       z = gadd(z, closure_evalnobrk(code));
     977         [ +  + ]:        560 :       if (j == e) break;
     978                 :            :     }
     979                 :        154 :     y = gmul(y, z);
     980                 :            :   }
     981                 :         42 :   pop_lex(1); return gerepileupto(av,y);
     982                 :            : }
     983                 :            : 
     984                 :            : /********************************************************************/
     985                 :            : /**                                                                **/
     986                 :            : /**                           PRODUCTS                             **/
     987                 :            : /**                                                                **/
     988                 :            : /********************************************************************/
     989                 :            : 
     990                 :            : GEN
     991                 :      47124 : produit(GEN a, GEN b, GEN code, GEN x)
     992                 :            : {
     993                 :      47124 :   pari_sp av, av0 = avma, lim;
     994                 :            :   GEN p1;
     995                 :            : 
     996         [ -  + ]:      47124 :   if (typ(a) != t_INT) pari_err_TYPE("prod",a);
     997         [ +  + ]:      47124 :   if (!x) x = gen_1;
     998         [ -  + ]:      47124 :   if (gcmp(b,a) < 0) return gcopy(x);
     999                 :            : 
    1000                 :      47124 :   b = gfloor(b);
    1001                 :      47124 :   a = setloop(a);
    1002                 :      47124 :   av=avma; lim = stack_lim(av,1);
    1003                 :      47124 :   push_lex(a,code);
    1004                 :            :   for(;;)
    1005                 :            :   {
    1006                 :     189147 :     p1 = closure_evalnobrk(code);
    1007         [ +  + ]:     189147 :     x = gmul(x,p1); if (cmpii(a,b) >= 0) break;
    1008                 :     142023 :     a = incloop(a);
    1009         [ -  + ]:     142023 :     if (low_stack(lim, stack_lim(av,1)))
    1010                 :            :     {
    1011         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prod");
    1012                 :          0 :       x = gerepileupto(av,x);
    1013                 :            :     }
    1014                 :     142023 :     set_lex(-1,a);
    1015                 :     142023 :   }
    1016                 :      47124 :   pop_lex(1); return gerepileupto(av0,x);
    1017                 :            : }
    1018                 :            : 
    1019                 :            : GEN
    1020                 :         14 : prodinf(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1021                 :            : {
    1022                 :         14 :   pari_sp av0 = avma, av, lim;
    1023                 :            :   long fl,G;
    1024                 :         14 :   GEN p1,x = real_1(prec);
    1025                 :            : 
    1026         [ -  + ]:         14 :   if (typ(a) != t_INT) pari_err_TYPE("prodinf",a);
    1027                 :         14 :   a = setloop(a);
    1028                 :         14 :   av = avma; lim = stack_lim(av,1);
    1029                 :         14 :   fl=0; G = -prec2nbits(prec)-5;
    1030                 :            :   for(;;)
    1031                 :            :   {
    1032         [ -  + ]:       1904 :     p1 = eval(E, a); if (gequal0(p1)) { x = p1; break; }
    1033                 :       1904 :     x = gmul(x,p1); a = incloop(a);
    1034                 :       1904 :     p1 = gsubgs(p1, 1);
    1035 [ +  - ][ +  + ]:       1904 :     if (gequal0(p1) || gexpo(p1) <= G) { if (++fl==3) break; } else fl=0;
                 [ +  + ]
    1036         [ -  + ]:       1890 :     if (low_stack(lim, stack_lim(av,1)))
    1037                 :            :     {
    1038         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodinf");
    1039                 :          0 :       x = gerepileupto(av,x);
    1040                 :            :     }
    1041                 :       1890 :   }
    1042                 :         14 :   return gerepilecopy(av0,x);
    1043                 :            : }
    1044                 :            : GEN
    1045                 :         14 : prodinf1(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1046                 :            : {
    1047                 :         14 :   pari_sp av0 = avma, av, lim;
    1048                 :            :   long fl,G;
    1049                 :         14 :   GEN p1,p2,x = real_1(prec);
    1050                 :            : 
    1051         [ -  + ]:         14 :   if (typ(a) != t_INT) pari_err_TYPE("prodinf1",a);
    1052                 :         14 :   a = setloop(a);
    1053                 :         14 :   av = avma; lim = stack_lim(av,1);
    1054                 :         14 :   fl=0; G = -prec2nbits(prec)-5;
    1055                 :            :   for(;;)
    1056                 :            :   {
    1057                 :       1904 :     p2 = eval(E, a); p1 = gaddgs(p2,1);
    1058         [ -  + ]:       1904 :     if (gequal0(p1)) { x = p1; break; }
    1059                 :       1904 :     x = gmul(x,p1); a = incloop(a);
    1060 [ +  - ][ +  + ]:       1904 :     if (gequal0(p2) || gexpo(p2) <= G) { if (++fl==3) break; } else fl=0;
                 [ +  + ]
    1061         [ -  + ]:       1890 :     if (low_stack(lim, stack_lim(av,1)))
    1062                 :            :     {
    1063         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodinf1");
    1064                 :          0 :       x = gerepileupto(av,x);
    1065                 :            :     }
    1066                 :       1890 :   }
    1067                 :         14 :   return gerepilecopy(av0,x);
    1068                 :            : }
    1069                 :            : GEN
    1070                 :         28 : prodinf0(GEN a, GEN code, long flag, long prec)
    1071                 :            : {
    1072      [ +  +  - ]:         28 :   switch(flag)
    1073                 :            :   {
    1074                 :         14 :     case 0: EXPR_WRAP(code, prodinf (EXPR_ARG, a, prec));
    1075                 :         14 :     case 1: EXPR_WRAP(code, prodinf1(EXPR_ARG, a, prec));
    1076                 :            :   }
    1077                 :          0 :   pari_err_FLAG("prodinf");
    1078                 :         28 :   return NULL; /* not reached */
    1079                 :            : }
    1080                 :            : 
    1081                 :            : GEN
    1082                 :         14 : prodeuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
    1083                 :            : {
    1084                 :         14 :   pari_sp av, av0 = avma, lim;
    1085                 :         14 :   GEN x = real_1(prec), prime;
    1086                 :            :   forprime_t T;
    1087                 :            : 
    1088                 :         14 :   av = avma;
    1089         [ -  + ]:         14 :   if (!forprime_init(&T, a,b)) { avma = av; return x; }
    1090                 :            : 
    1091                 :         14 :   av = avma;
    1092                 :         14 :   lim = stack_lim(avma,1);
    1093         [ +  + ]:      17220 :   while ( (prime = forprime_next(&T)) )
    1094                 :            :   {
    1095                 :      17206 :     x = gmul(x, eval(E, prime));
    1096         [ -  + ]:      17206 :     if (low_stack(lim, stack_lim(av,1)))
    1097                 :            :     {
    1098         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"prodeuler");
    1099                 :          0 :       x = gerepilecopy(av, x);
    1100                 :            :     }
    1101                 :            :   }
    1102                 :         14 :   return gerepilecopy(av0,x);
    1103                 :            : }
    1104                 :            : GEN
    1105                 :         14 : prodeuler0(GEN a, GEN b, GEN code, long prec)
    1106                 :         14 : { EXPR_WRAP(code, prodeuler(EXPR_ARG, a, b, prec)); }
    1107                 :            : 
    1108                 :            : static void
    1109                 :          7 : err_direuler(GEN x)
    1110                 :          7 : { pari_err_DOMAIN("direuler","constant term","!=", gen_1,x); }
    1111                 :            : 
    1112                 :            : GEN
    1113                 :         28 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
    1114                 :            : {
    1115                 :            :   ulong i, k, n;
    1116                 :         28 :   pari_sp av0 = avma, av, lim = stack_lim(av0, 1);
    1117                 :            :   long j, tx, lx;
    1118                 :            :   GEN x, y, s, polnum, polden, prime;
    1119                 :            :   forprime_t T;
    1120                 :            : 
    1121         [ -  + ]:         28 :   if (c)
    1122                 :            :   {
    1123         [ #  # ]:          0 :     if (typ(c) != t_INT)
    1124                 :            :     {
    1125                 :          0 :       c = gfloor(c);
    1126         [ #  # ]:          0 :       if (typ(c) != t_INT) pari_err_TYPE("direuler", c);
    1127                 :            :     }
    1128         [ #  # ]:          0 :     if (signe(c) <= 0) { avma = av0; return mkvec(gen_1); }
    1129                 :          0 :     n = itou(c);
    1130         [ #  # ]:          0 :     if (cmpui(n, b) < 0) b = c;
    1131                 :            :   }
    1132         [ -  + ]:         28 :   if (!forprime_init(&T, a,b)) { avma = av0; return mkvec(gen_1); }
    1133                 :            : 
    1134         [ -  + ]:         28 :   if (c)
    1135                 :            :   {
    1136                 :          0 :     n = itou(c);
    1137         [ #  # ]:          0 :     if (cmpui(n, b) < 0) b = c;
    1138                 :            :   }
    1139                 :            :   else
    1140                 :            :   {
    1141         [ -  + ]:         28 :     if (lgefint(b) > 3) pari_err_OVERFLOW("direuler");
    1142                 :         28 :     n = itou(b);
    1143                 :            :   }
    1144                 :            : 
    1145                 :         28 :   y = cgetg(n+1,t_VEC); av = avma;
    1146                 :         28 :   x = zerovec(n); gel(x,1) = gen_1;
    1147         [ +  + ]:        462 :   while ( (prime = forprime_next(&T)) )
    1148                 :            :   {
    1149                 :        441 :     ulong p = prime[2];
    1150                 :        441 :     s = eval(E, prime);
    1151                 :        441 :     polnum = numer(s);
    1152                 :        441 :     polden = denom(s);
    1153                 :        441 :     tx = typ(polnum);
    1154         [ +  + ]:        441 :     if (is_scalar_t(tx))
    1155                 :            :     {
    1156         [ +  + ]:        357 :       if (!gequal1(polnum))
    1157                 :            :       {
    1158         [ +  - ]:          7 :         if (!gequalm1(polnum)) err_direuler(polnum);
    1159                 :          0 :         polden = gneg(polden);
    1160                 :            :       }
    1161                 :            :     }
    1162                 :            :     else
    1163                 :            :     {
    1164                 :            :       ulong k1, q, qlim;
    1165         [ -  + ]:         84 :       if (tx != t_POL) pari_err_TYPE("direuler",polnum);
    1166                 :         84 :       lx = degpol(polnum);
    1167         [ -  + ]:         84 :       if (lx < 0) err_direuler(polnum);
    1168                 :         84 :       c = gel(polnum,2);
    1169         [ -  + ]:         84 :       if (!gequal1(c))
    1170                 :            :       {
    1171         [ #  # ]:          0 :         if (!gequalm1(c)) err_direuler(polnum);
    1172                 :          0 :         polnum = gneg(polnum);
    1173                 :          0 :         polden = gneg(polden);
    1174                 :            :       }
    1175         [ +  + ]:       3444 :       for (i=1; i<=n; i++) gel(y,i) = gel(x,i);
    1176                 :         84 :       q = p; qlim = n/p;
    1177 [ +  - ][ +  + ]:        105 :       for (j = 1; q<=n && j<=lx; j++)
    1178                 :            :       {
    1179                 :         84 :         c = gel(polnum,j+2);
    1180         [ +  - ]:         84 :         if (!gequal0(c))
    1181         [ +  + ]:        504 :           for (k=1,k1=q; k1<=n; k1+=q,k++)
    1182                 :        420 :             gel(x,k1) = gadd(gel(x,k1), gmul(c,gel(y,k)));
    1183         [ +  + ]:         84 :         if (q > qlim) break;
    1184                 :         21 :         q *= p;
    1185                 :            :       }
    1186                 :            :     }
    1187                 :        434 :     tx = typ(polden);
    1188         [ -  + ]:        434 :     if (is_scalar_t(tx))
    1189                 :            :     {
    1190         [ #  # ]:          0 :       if (!gequal1(polden))
    1191                 :          0 :         pari_err_DOMAIN("direuler","constant term", "!=", gen_1,polden);
    1192                 :            :     }
    1193                 :            :     else
    1194                 :            :     {
    1195         [ -  + ]:        434 :       if (tx != t_POL) pari_err_TYPE("direuler",polden);
    1196                 :        434 :       c = gel(polden,2);
    1197         [ -  + ]:        434 :       if (!gequal1(c))
    1198                 :          0 :         pari_err_DOMAIN("direuler","constant term", "!=", gen_1,polden);
    1199                 :        434 :       lx = degpol(polden);
    1200         [ +  + ]:       3248 :       for (i=p; i<=n; i+=p)
    1201                 :            :       {
    1202                 :       2814 :         s = gen_0; k = i;
    1203 [ +  + ][ +  + ]:       5796 :         for (j = 1; !(k%p) && j<=lx; j++)
    1204                 :            :         {
    1205                 :       2982 :           c = gel(polden,j+2); k /= p;
    1206         [ +  + ]:       2982 :           if (!gequal0(c)) s = gadd(s, gmul(c,gel(x,k)));
    1207                 :            :         }
    1208                 :       2814 :         gel(x,i) = gsub(gel(x,i),s);
    1209                 :            :       }
    1210                 :            :     }
    1211         [ -  + ]:        434 :     if (low_stack(lim, stack_lim(av,1)))
    1212                 :            :     {
    1213         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"direuler");
    1214                 :          0 :       x = gerepilecopy(av, x);
    1215                 :            :     }
    1216                 :            :   }
    1217                 :         21 :   return gerepilecopy(av0,x);
    1218                 :            : }
    1219                 :            : GEN
    1220                 :         28 : direuler0(GEN a, GEN b, GEN code, GEN c)
    1221                 :         28 : { EXPR_WRAP(code, direuler(EXPR_ARG, a, b, c)); }
    1222                 :            : 
    1223                 :            : /********************************************************************/
    1224                 :            : /**                                                                **/
    1225                 :            : /**                       VECTORS & MATRICES                       **/
    1226                 :            : /**                                                                **/
    1227                 :            : /********************************************************************/
    1228                 :            : 
    1229                 :            : INLINE GEN
    1230                 :     678234 : copyupto(GEN z, GEN t)
    1231                 :            : {
    1232 [ +  + ][ +  - ]:     678234 :   if (is_universal_constant(z) || (z>(GEN)pari_mainstack->bot && z<=t))
                 [ +  + ]
    1233                 :     673306 :     return z;
    1234                 :            :   else
    1235                 :     678234 :     return gcopy(z);
    1236                 :            : }
    1237                 :            : 
    1238                 :            : GEN
    1239                 :        805 : vecexpr0(GEN vec, GEN code, GEN pred)
    1240                 :            : {
    1241      [ -  +  - ]:        805 :   switch(typ(vec))
    1242                 :            :   {
    1243                 :            :     case t_LIST:
    1244                 :            :     {
    1245                 :          0 :       vec = list_data(vec);
    1246         [ #  # ]:          0 :       if (!vec) return cgetg(1, t_VEC);
    1247                 :          0 :       break;
    1248                 :            :     }
    1249                 :        805 :     case t_VEC: case t_COL: case t_MAT: break;
    1250                 :          0 :     default: pari_err_TYPE("[_|_<-_,_]",vec);
    1251                 :            :   }
    1252 [ +  + ][ +  - ]:        805 :   if (pred && code)
    1253                 :        175 :     EXPR_WRAP(code,vecselapply((void*)pred,&gp_evalbool,EXPR_ARGUPTO,vec))
    1254         [ +  - ]:        630 :   else if (code)
    1255                 :        630 :     EXPR_WRAP(code,vecapply(EXPR_ARGUPTO,vec))
    1256                 :            :   else
    1257                 :        805 :     EXPR_WRAP(pred,vecselect(EXPR_ARGBOOL,vec))
    1258                 :            : }
    1259                 :            : 
    1260                 :            : GEN
    1261                 :        154 : vecexpr1(GEN vec, GEN code, GEN pred)
    1262                 :            : {
    1263                 :        154 :   GEN v = vecexpr0(vec, code, pred);
    1264         [ +  + ]:        154 :   return lg(v) == 1? v: shallowconcat1(v);
    1265                 :            : }
    1266                 :            : 
    1267                 :            : GEN
    1268                 :    2323302 : vecteur(GEN nmax, GEN code)
    1269                 :            : {
    1270                 :            :   GEN y,p1;
    1271                 :            :   long i,m;
    1272                 :    2323302 :   GEN c=utoipos(1);
    1273                 :            : 
    1274                 :    2323302 :   m = gtos(nmax);
    1275         [ +  + ]:    2323302 :   if (m < 0)  pari_err_DOMAIN("vector", "dimension", "<", gen_0, stoi(m));
    1276         [ +  + ]:    2323288 :   if (!code) return zerovec(m);
    1277                 :       3404 :   y = cgetg(m+1,t_VEC); push_lex(c, code);
    1278         [ +  + ]:     112003 :   for (i=1; i<=m; i++)
    1279                 :            :   {
    1280                 :     108606 :     c[2] = i; p1 = closure_evalnobrk(code);
    1281                 :     108599 :     gel(y,i) = copyupto(p1, y);
    1282                 :     108599 :     set_lex(-1,c);
    1283                 :            :   }
    1284                 :    2323281 :   pop_lex(1); return y;
    1285                 :            : }
    1286                 :            : 
    1287                 :            : GEN
    1288                 :          7 : vecteursmall(GEN nmax, GEN code)
    1289                 :            : {
    1290                 :            :   GEN y;
    1291                 :            :   long i,m;
    1292                 :          7 :   GEN c=utoipos(1);
    1293                 :            : 
    1294                 :          7 :   m = gtos(nmax);
    1295         [ +  - ]:          7 :   if (m < 0)  pari_err_DOMAIN("vectorsmall", "dimension", "<", gen_0, stoi(m));
    1296         [ #  # ]:          0 :   if (!code) return zero_zv(m);
    1297                 :          0 :   y = cgetg(m+1,t_VECSMALL); push_lex(c,code);
    1298         [ #  # ]:          0 :   for (i=1; i<=m; i++)
    1299                 :            :   {
    1300                 :          0 :     c[2] = i;
    1301                 :          0 :     y[i] = gtos(closure_evalnobrk(code));
    1302                 :          0 :     set_lex(-1,c);
    1303                 :            :   }
    1304                 :          0 :   pop_lex(1); return y;
    1305                 :            : }
    1306                 :            : 
    1307                 :            : GEN
    1308                 :        504 : vvecteur(GEN nmax, GEN n)
    1309                 :            : {
    1310                 :        504 :   GEN y = vecteur(nmax,n);
    1311                 :        497 :   settyp(y,t_COL); return y;
    1312                 :            : }
    1313                 :            : 
    1314                 :            : GEN
    1315                 :        478 : matrice(GEN nlig, GEN ncol, GEN code)
    1316                 :            : {
    1317                 :            :   GEN y, z, p1;
    1318                 :            :   long i, j, m, n;
    1319                 :        478 :   GEN c1 = utoipos(1), c2 = utoipos(1);
    1320                 :            : 
    1321                 :        478 :   m = gtos(ncol);
    1322                 :        478 :   n = gtos(nlig);
    1323         [ +  + ]:        478 :   if (m < 0)  pari_err_DOMAIN("matrix", "nbcols", "<", gen_0, stoi(m));
    1324         [ +  + ]:        471 :   if (n < 0)  pari_err_DOMAIN("matrix", "nbrows", "<", gen_0, stoi(n));
    1325         [ +  + ]:        464 :   if (!m) return cgetg(1,t_MAT);
    1326 [ +  + ][ -  + ]:        422 :   if (!code || !n) return zeromatcopy(n, m);
    1327                 :        373 :   push_lex(c1,code);
    1328                 :        373 :   push_lex(c2,NULL); y = cgetg(m+1,t_MAT);
    1329         [ +  + ]:       5408 :   for (i=1; i<=m; i++)
    1330                 :            :   {
    1331                 :       5035 :     c2[2] = i; z = cgetg(n+1,t_COL); gel(y,i) = z;
    1332         [ +  + ]:     574670 :     for (j=1; j<=n; j++)
    1333                 :            :     {
    1334                 :     569635 :       c1[2] = j; p1 = closure_evalnobrk(code);
    1335                 :     569635 :       gel(z,j) = copyupto(p1, y);
    1336                 :     569635 :       set_lex(-2,c1);
    1337                 :     569635 :       set_lex(-1,c2);
    1338                 :            :     }
    1339                 :            :   }
    1340                 :        464 :   pop_lex(2); return y;
    1341                 :            : }
    1342                 :            : 
    1343                 :            : /********************************************************************/
    1344                 :            : /**                                                                **/
    1345                 :            : /**                         SUMMING SERIES                         **/
    1346                 :            : /**                                                                **/
    1347                 :            : /********************************************************************/
    1348                 :            : GEN
    1349                 :         35 : polzag(long n, long m)
    1350                 :            : {
    1351                 :         35 :   pari_sp av = avma;
    1352                 :         35 :   long k, d = n - m;
    1353                 :            :   GEN A, Bx, g, s;
    1354                 :            : 
    1355 [ +  - ][ -  + ]:         35 :   if (d <= 0 || m < 0) return gen_0;
    1356                 :         35 :   A  = mkpoln(2, stoi(-2), gen_1); /* 1 - 2x */
    1357                 :         35 :   Bx = mkpoln(3, stoi(-2), gen_2, gen_0); /* 2x - 2x^2 */
    1358                 :         35 :   g = gmul(poleval(ZX_deriv(polchebyshev1(d,0)), A), gpowgs(Bx, (m+1)>>1));
    1359         [ +  + ]:        147 :   for (k = m; k >= 0; k--)
    1360         [ +  + ]:        112 :     g = (k&1)? ZX_deriv(g): gadd(gmul(A,g), gmul(Bx,ZX_deriv(g)));
    1361                 :         35 :   s = mulii(sqru(d), mpfact(m+1));
    1362                 :         35 :   return gerepileupto(av, gdiv(g,s));
    1363                 :            : }
    1364                 :            : 
    1365                 :            : #ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */
    1366                 :            : #pragma optimize("g",off)
    1367                 :            : #endif
    1368                 :            : static GEN
    1369                 :         14 : polzagreel(long n, long m, long prec)
    1370                 :            : {
    1371                 :         14 :   const long d = n - m, d2 = d<<1, r = (m+1)>>1;
    1372                 :            :   long j, k, k2;
    1373                 :         14 :   pari_sp av = avma;
    1374                 :            :   GEN Bx, g, h, v, b, s;
    1375                 :            : 
    1376 [ +  - ][ -  + ]:         14 :   if (d <= 0 || m < 0) return gen_0;
    1377                 :         14 :   Bx = mkpoln(3, gen_1, gen_1, gen_0); /* x + x^2 */
    1378                 :         14 :   v = cgetg(d+1,t_VEC);
    1379                 :         14 :   g = cgetg(d+1,t_VEC);
    1380                 :         14 :   gel(v,d) = gen_1; b = stor(d2, prec);
    1381                 :         14 :   gel(g,d) = b;
    1382         [ +  + ]:        294 :   for (k = 1; k < d; k++)
    1383                 :            :   {
    1384                 :        280 :     gel(v,d-k) = gen_1;
    1385         [ +  + ]:       2940 :     for (j=1; j<k; j++)
    1386                 :       2660 :       gel(v,d-k+j) = addii(gel(v,d-k+j), gel(v,d-k+j+1));
    1387                 :            :     /* v[d-k+j] = binom(k, j), j = 0..k */
    1388                 :        280 :     k2 = k+k; b = divri(mulri(b,mulss(d2-k2+1,d2-k2)), mulss(k2,k2+1));
    1389         [ +  + ]:       3220 :     for (j=1; j<=k; j++)
    1390                 :       2940 :       gel(g,d-k+j) = mpadd(gel(g,d-k+j), mulri(b,gel(v,d-k+j)));
    1391                 :        280 :     gel(g,d-k) = b;
    1392                 :            :   }
    1393                 :         14 :   g = gmul(RgV_to_RgX(g,0), gpowgs(Bx,r));
    1394         [ +  + ]:        168 :   for (j=0; j<=r; j++)
    1395                 :            :   {
    1396         [ +  + ]:        154 :     if (j) g = RgX_deriv(g);
    1397 [ +  + ][ +  - ]:        154 :     if (j || !(m&1))
    1398                 :            :     {
    1399                 :        154 :       h = cgetg(n+3,t_POL);
    1400                 :        154 :       h[1] = evalsigne(1);
    1401                 :        154 :       gel(h,2) = gel(g,2);
    1402         [ +  + ]:       6314 :       for (k=1; k<n; k++)
    1403                 :       6160 :         gel(h,k+2) = gadd(gmulsg(k+k+1,gel(g,k+2)), gmulsg(k<<1,gel(g,k+1)));
    1404                 :        154 :       gel(h,n+2) = gmulsg(n<<1, gel(g,n+1));
    1405                 :        154 :       g = h;
    1406                 :            :     }
    1407                 :            :   }
    1408                 :         14 :   g = gmul2n(g, r-1);
    1409                 :         14 :   s = mului(d, mpfact(m+1));
    1410                 :         14 :   return gerepileupto(av, gdiv(g,s));
    1411                 :            : }
    1412                 :            : 
    1413                 :            : #ifdef _MSC_VER
    1414                 :            : #pragma optimize("g",on)
    1415                 :            : #endif
    1416                 :            : GEN
    1417                 :         21 : sumalt(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1418                 :            : {
    1419                 :            :   ulong k, N;
    1420                 :         21 :   pari_sp av = avma, av2, lim;
    1421                 :            :   GEN s, az, c, e1, d;
    1422                 :            : 
    1423         [ -  + ]:         21 :   if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
    1424                 :         21 :   e1 = addsr(3, sqrtr(stor(8,prec)));
    1425                 :         21 :   N = (ulong)(0.4*(prec2nbits(prec)+ 7));
    1426                 :         21 :   d = powru(e1,N);
    1427                 :         21 :   d = shiftr(addrr(d, invr(d)),-1);
    1428                 :         21 :   a = setloop(a);
    1429                 :         21 :   az = gen_m1; c = d;
    1430                 :         21 :   s = gen_0;
    1431                 :         21 :   av2 = avma; lim = stack_lim(av,4);
    1432                 :         21 :   for (k=0; ; k++) /* k < N */
    1433                 :            :   {
    1434                 :       3640 :     c = addir(az,c); s = gadd(s, gmul(c, eval(E, a)));
    1435         [ +  + ]:       3640 :     if (k==N-1) break;
    1436                 :       3619 :     az = diviuuexact(muluui((N-k)<<1,N+k,az), k+1, (k<<1)+1);
    1437                 :       3619 :     a = incloop(a); /* in place! */
    1438         [ -  + ]:       3619 :     if (low_stack(lim, stack_lim(av,4)))
    1439                 :            :     {
    1440         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumalt, k = %ld/%ld", k,N-1);
    1441                 :          0 :       gerepileall(av2, 3, &az,&c,&s);
    1442                 :            :     }
    1443                 :       3619 :   }
    1444                 :         21 :   return gerepileupto(av, gdiv(s,d));
    1445                 :            : }
    1446                 :            : 
    1447                 :            : GEN
    1448                 :         14 : sumalt2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1449                 :            : {
    1450                 :            :   long k, N;
    1451                 :         14 :   pari_sp av = avma, av2, lim;
    1452                 :            :   GEN s, dn, pol;
    1453                 :            : 
    1454         [ -  + ]:         14 :   if (typ(a) != t_INT) pari_err_TYPE("sumalt",a);
    1455                 :         14 :   N = (long)(0.31*(prec2nbits(prec) + 5));
    1456                 :         14 :   pol = polzagreel(N,N>>1,prec+EXTRAPRECWORD);
    1457                 :         14 :   pol = RgX_div_by_X_x(pol, gen_1, &dn);
    1458                 :         14 :   a = setloop(a);
    1459                 :         14 :   N = degpol(pol);
    1460                 :         14 :   s = gen_0;
    1461                 :         14 :   av2 = avma; lim = stack_lim(av,4);
    1462         [ +  - ]:        574 :   for (k=0; k<=N; k++)
    1463                 :            :   {
    1464                 :        574 :     s = gadd(s, gmul(gel(pol,k+2), eval(E, a)));
    1465         [ +  + ]:        574 :     if (k == N) break;
    1466                 :        560 :     a = incloop(a); /* in place! */
    1467         [ -  + ]:        560 :     if (low_stack(lim, stack_lim(av,4)))
    1468                 :            :     {
    1469         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"sumalt2, k = %ld/%ld", k,N-1);
    1470                 :          0 :       s = gerepileupto(av2, s);
    1471                 :            :     }
    1472                 :            :   }
    1473                 :         14 :   return gerepileupto(av, gdiv(s,dn));
    1474                 :            : }
    1475                 :            : 
    1476                 :            : GEN
    1477                 :         35 : sumalt0(GEN a, GEN code, long flag, long prec)
    1478                 :            : {
    1479      [ +  +  - ]:         35 :   switch(flag)
    1480                 :            :   {
    1481                 :         21 :     case 0: EXPR_WRAP(code, sumalt (EXPR_ARG,a,prec));
    1482                 :         14 :     case 1: EXPR_WRAP(code, sumalt2(EXPR_ARG,a,prec));
    1483                 :          0 :     default: pari_err_FLAG("sumalt");
    1484                 :            :   }
    1485                 :         35 :   return NULL; /* not reached */
    1486                 :            : }
    1487                 :            : 
    1488                 :            : GEN
    1489                 :         42 : sumpos(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1490                 :            : {
    1491                 :            :   long k, kk, N, G;
    1492                 :         42 :   pari_sp av = avma;
    1493                 :            :   GEN r, reel, s, az, c, x, e1, d, *stock;
    1494                 :            : 
    1495         [ -  + ]:         42 :   if (typ(a) != t_INT) pari_err_TYPE("sumpos",a);
    1496                 :         42 :   a = subis(a,1); reel = cgetr(prec);
    1497                 :         42 :   e1 = addsr(3, sqrtr(stor(8,prec)));
    1498                 :         42 :   N = (long)(0.4*(prec2nbits(prec) + 7));
    1499                 :         42 :   d = powru(e1,N);
    1500                 :         42 :   d = shiftr(addrr(d, invr(d)),-1);
    1501                 :         42 :   az = gen_m1; c = d;
    1502                 :         42 :   s = gen_0;
    1503                 :            : 
    1504                 :         42 :   G = -prec2nbits(prec) - 5;
    1505         [ +  + ]:       9828 :   stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;
    1506         [ +  + ]:       9828 :   for (k=0; k<N; k++)
    1507                 :            :   {
    1508 [ +  + ][ +  + ]:       9786 :     if (odd(k) && stock[k]) x = stock[k];
    1509                 :            :     else
    1510                 :            :     {
    1511                 :       6531 :       pari_sp av2 = avma;
    1512                 :       6531 :       x = gen_0; r = utoipos(2*k+2);
    1513                 :       6531 :       for(kk=0;;kk++)
    1514                 :            :       {
    1515                 :            : 
    1516                 :            :         long ex;
    1517                 :    3021004 :         affgr(eval(E, addii(r,a)), reel);
    1518         [ -  + ]:    3021004 :         if (!signe(reel)) break;
    1519                 :    3021004 :         ex = expo(reel) + kk; shiftr_inplace(reel, kk);
    1520 [ +  + ][ +  + ]:    3021004 :         x = mpadd(x,reel); if (kk && ex < G) break;
    1521                 :    3014473 :         r = shifti(r,1);
    1522                 :    3014473 :       }
    1523                 :       6531 :       x = gerepileupto(av2, x);
    1524         [ +  + ]:       6531 :       if (2*k < N) stock[2*k+1] = x;
    1525                 :       6531 :       affgr(eval(E, addsi(k+1,a)), reel);
    1526                 :       6531 :       x = addrr(reel, gmul2n(x,1));
    1527                 :            :     }
    1528                 :       9786 :     c = addir(az,c);
    1529         [ +  + ]:       9786 :     s = mpadd(s, mulrr(x, k&1? negr(c): c));
    1530                 :       9786 :     az = diviiexact(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1));
    1531                 :            :   }
    1532                 :         42 :   return gerepileupto(av, gdiv(s,d));
    1533                 :            : }
    1534                 :            : 
    1535                 :            : GEN
    1536                 :          0 : sumpos2(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)
    1537                 :            : {
    1538                 :            :   long k, kk, N, G;
    1539                 :          0 :   pari_sp av = avma;
    1540                 :            :   GEN r, reel, s, pol, dn, x, *stock;
    1541                 :            : 
    1542         [ #  # ]:          0 :   if (typ(a) != t_INT) pari_err_TYPE("sumpos2",a);
    1543                 :          0 :   a = subis(a,1); reel = cgetr(prec);
    1544                 :          0 :   N = (long)(0.31*(prec2nbits(prec) + 5));
    1545                 :            : 
    1546                 :          0 :   G = -prec2nbits(prec) - 5;
    1547         [ #  # ]:          0 :   stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;
    1548         [ #  # ]:          0 :   for (k=1; k<=N; k++)
    1549 [ #  # ][ #  # ]:          0 :     if (odd(k) || !stock[k])
    1550                 :            :     {
    1551                 :          0 :       pari_sp av2 = avma;
    1552                 :          0 :       x = gen_0; r = utoipos(2*k);
    1553                 :          0 :       for(kk=0;;kk++)
    1554                 :            :       {
    1555                 :            :         long ex;
    1556                 :          0 :         affgr(eval(E, addii(r,a)), reel);
    1557                 :          0 :         ex = expo(reel) + kk; shiftr_inplace(reel, kk);
    1558 [ #  # ][ #  # ]:          0 :         x = mpadd(x,reel); if (kk && ex < G) break;
    1559                 :          0 :         r = shifti(r,1);
    1560                 :          0 :       }
    1561                 :          0 :       x = gerepileupto(av2, x);
    1562         [ #  # ]:          0 :       if (2*k-1 < N) stock[2*k] = x;
    1563                 :          0 :       affgr(eval(E, addsi(k,a)), reel);
    1564                 :          0 :       stock[k] = addrr(reel, gmul2n(x,1));
    1565                 :            :     }
    1566                 :          0 :   s = gen_0;
    1567                 :          0 :   pol = polzagreel(N,N>>1,prec+EXTRAPRECWORD);
    1568                 :          0 :   pol = RgX_div_by_X_x(pol, gen_1, &dn);
    1569         [ #  # ]:          0 :   for (k=1; k<=lg(pol)-2; k++)
    1570                 :            :   {
    1571                 :          0 :     GEN p1 = gmul(gel(pol,k+1),stock[k]);
    1572         [ #  # ]:          0 :     if (!odd(k)) p1 = gneg_i(p1);
    1573                 :          0 :     s = gadd(s,p1);
    1574                 :            :   }
    1575                 :          0 :   return gerepileupto(av, gdiv(s,dn));
    1576                 :            : }
    1577                 :            : 
    1578                 :            : GEN
    1579                 :         42 : sumpos0(GEN a, GEN code, long flag, long prec)
    1580                 :            : {
    1581      [ +  -  - ]:         42 :   switch(flag)
    1582                 :            :   {
    1583                 :         42 :     case 0: EXPR_WRAP(code, sumpos (EXPR_ARG,a,prec));
    1584                 :          0 :     case 1: EXPR_WRAP(code, sumpos2(EXPR_ARG,a,prec));
    1585                 :          0 :     default: pari_err_FLAG("sumpos");
    1586                 :            :   }
    1587                 :         42 :   return NULL; /* not reached */
    1588                 :            : }
    1589                 :            : 
    1590                 :            : /********************************************************************/
    1591                 :            : /**                                                                **/
    1592                 :            : /**            SEARCH FOR REAL ZEROS of an expression              **/
    1593                 :            : /**                                                                **/
    1594                 :            : /********************************************************************/
    1595                 :            : /* Brent's method, [a,b] bracketing interval */
    1596                 :            : GEN
    1597                 :         35 : zbrent(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
    1598                 :            : {
    1599                 :            :   long sig, iter, itmax;
    1600                 :         35 :   pari_sp av = avma;
    1601                 :            :   GEN c, d, e, tol, fa, fb, fc;
    1602                 :            : 
    1603                 :         35 :   a = gtofp(a,prec);
    1604                 :         35 :   b = gtofp(b,prec); sig = cmprr(b,a);
    1605         [ -  + ]:         35 :   if (!sig) return gerepileupto(av, a);
    1606         [ -  + ]:         35 :   if (sig < 0) { c=a; a=b; b=c; } else c = b;
    1607                 :            : 
    1608                 :         35 :   fa = eval(E, a);
    1609                 :         35 :   fb = eval(E, b);
    1610         [ +  + ]:         35 :   if (gsigne(fa)*gsigne(fb) > 0)
    1611                 :          7 :     pari_err_DOMAIN("solve", "f(a)f(b)", ">", gen_0, mkvec2(fa,fb));
    1612                 :         28 :   itmax = prec2nbits(prec) * 2 + 1;
    1613                 :         28 :   tol = real2n(5-prec2nbits(prec), LOWDEFAULTPREC);
    1614                 :         28 :   fc = fb;
    1615                 :         28 :   e = d = NULL; /* gcc -Wall */
    1616         [ +  - ]:        609 :   for (iter=1; iter<=itmax; iter++)
    1617                 :            :   {
    1618                 :            :     GEN xm, tol1;
    1619         [ +  + ]:        609 :     if (gsigne(fb)*gsigne(fc) > 0)
    1620                 :            :     {
    1621                 :        185 :       c = a; fc = fa; e = d = subrr(b,a);
    1622                 :            :     }
    1623         [ +  + ]:        609 :     if (gcmp(gabs(fc,0),gabs(fb,0)) < 0)
    1624                 :            :     {
    1625                 :         15 :       a = b; b = c; c = a; fa = fb; fb = fc; fc = fa;
    1626                 :            :     }
    1627         [ -  + ]:        609 :     tol1 = absr_cmp(tol, b) > 0? sqrr(tol): mulrr(tol, absr(b));
    1628                 :        609 :     xm = shiftr(subrr(c,b),-1);
    1629 [ +  + ][ +  + ]:        609 :     if (absr_cmp(xm,tol1) <= 0 || gequal0(fb)) break; /* SUCCESS */
    1630                 :            : 
    1631 [ +  - ][ +  - ]:        581 :     if (absr_cmp(e,tol1) >= 0 && gcmp(gabs(fa,0),gabs(fb,0)) > 0)
    1632                 :        581 :     { /* attempt interpolation */
    1633                 :        581 :       GEN min1, min2, p, q, s = gdiv(fb,fa);
    1634         [ +  + ]:        581 :       if (cmprr(a,c)==0)
    1635                 :            :       {
    1636                 :        170 :         p = gmul2n(gmul(xm,s),1);
    1637                 :        170 :         q = gsubsg(1,s);
    1638                 :            :       }
    1639                 :            :       else
    1640                 :            :       {
    1641                 :        411 :         GEN r = gdiv(fb,fc);
    1642                 :        411 :         q = gdiv(fa,fc);
    1643                 :        411 :         p = gmul2n(gmul(gsub(q,r),gmul(xm,q)),1);
    1644                 :        411 :         p = gmul(s, gsub(p, gmul(gsub(b,a),gsubgs(r,1))));
    1645                 :        411 :         q = gmul(gmul(gsubgs(q,1),gsubgs(r,1)),gsubgs(s,1));
    1646                 :            :       }
    1647         [ +  + ]:        581 :       if (gsigne(p) > 0) q = gneg_i(q); else p = gneg_i(p);
    1648                 :        581 :       min1 = gsub(gmulsg(3,gmul(xm,q)), gabs(gmul(q,tol1),0));
    1649                 :        581 :       min2 = gabs(gmul(e,q),0);
    1650         [ +  + ]:        581 :       if (gcmp(gmul2n(p,1), gmin(min1,min2)) < 0)
    1651                 :        210 :         { e = d; d = gdiv(p,q); } /* interpolation OK */
    1652                 :            :       else
    1653                 :        371 :         { d = xm; e = d; } /* failed, use bisection */
    1654                 :            :     }
    1655                 :          0 :     else { d = xm; e = d; } /* bound decreasing too slowly, use bisection */
    1656                 :        581 :     a = b; fa = fb;
    1657         [ +  + ]:        581 :     if (gcmp(gabs(d,0),tol1) > 0) b = gadd(b,d);
    1658         [ +  + ]:         29 :     else if (gsigne(xm) > 0)      b = addrr(b,tol1);
    1659                 :         27 :     else                          b = subrr(b,tol1);
    1660                 :        581 :     fb = eval(E, b);
    1661                 :            :   }
    1662         [ -  + ]:         28 :   if (iter > itmax) pari_err_IMPL("solve recovery [too many iterations]");
    1663                 :         28 :   return gerepileuptoleaf(av, rcopy(b));
    1664                 :            : }
    1665                 :            : 
    1666                 :            : GEN
    1667                 :         21 : zbrent0(GEN a, GEN b, GEN code, long prec)
    1668                 :         21 : { EXPR_WRAP(code, zbrent(EXPR_ARG, a,b, prec)); }
    1669                 :            : 
    1670                 :            : /* x = solve_start(&D, a, b, prec)
    1671                 :            :  * while (x) {
    1672                 :            :  *   y = ...(x);
    1673                 :            :  *   x = solve_next(&D, y);
    1674                 :            :  * }
    1675                 :            :  * return D.res; */
    1676                 :            : 
    1677                 :            : /********************************************************************/
    1678                 :            : /**                                                                **/
    1679                 :            : /**            Numerical derivation                                **/
    1680                 :            : /**                                                                **/
    1681                 :            : /********************************************************************/
    1682                 :            : 
    1683                 :            : /* Rationale: (f(2^-e) - f(-2^-e) + O(2^-pr)) / (2 * 2^-e) = f'(0) + O(2^-2e)
    1684                 :            :  * since 2nd derivatives cancel.
    1685                 :            :  *   prec(LHS) = pr - e
    1686                 :            :  *   prec(RHS) = 2e, equal when  pr = 3e = 3/2 fpr (fpr = required final prec)
    1687                 :            :  *
    1688                 :            :  * For f'(x), x far from 0: prec(LHS) = pr - e - expo(x)
    1689                 :            :  * --> pr = 3/2 fpr + expo(x) */
    1690                 :            : GEN
    1691                 :         21 : derivnum(void *E, GEN (*eval)(void *, GEN), GEN x, long prec)
    1692                 :            : {
    1693                 :            :   GEN eps,a,b, y;
    1694                 :            :   long pr, l, e, ex;
    1695                 :         21 :   pari_sp av = avma;
    1696                 :         21 :   long p = precision(x);
    1697         [ -  + ]:         21 :   long fpr = p ? prec2nbits(p): prec2nbits(prec);
    1698                 :         21 :   ex = gexpo(x);
    1699         [ -  + ]:         21 :   if (ex < 0) ex = 0; /* near 0 */
    1700                 :         21 :   pr = (long)ceil(fpr * 1.5 + ex);
    1701                 :         21 :   l = nbits2prec(pr);
    1702         [ -  + ]:         21 :   switch(typ(x))
    1703                 :            :   {
    1704                 :            :     case t_REAL:
    1705                 :            :     case t_COMPLEX:
    1706                 :          0 :       x = gprec_w(x, l + nbits2extraprec(ex + BITS_IN_LONG));
    1707                 :            :   }
    1708                 :            : 
    1709                 :         21 :   e = fpr/2; /* 1/2 required prec (in sig. bits) */
    1710                 :         21 :   eps = real2n(-e, l);
    1711                 :         21 :   a = eval(E, gsub(x, eps));
    1712                 :         21 :   b = eval(E, gadd(x, eps));
    1713                 :         21 :   y = gmul2n(gsub(b,a), e-1);
    1714                 :         21 :   return gerepileupto(av, gprec_w(y, nbits2prec(fpr)));
    1715                 :            : }
    1716                 :            : 
    1717                 :            : GEN
    1718                 :         28 : derivfun(void *E, GEN (*eval)(void *, GEN), GEN x, long prec)
    1719                 :            : {
    1720                 :         28 :   pari_sp av = avma;
    1721                 :            :   long vx;
    1722   [ +  -  +  - ]:         28 :   switch(typ(x))
    1723                 :            :   {
    1724                 :            :   case t_REAL: case t_INT: case t_FRAC: case t_COMPLEX:
    1725                 :         21 :     return derivnum(E,eval, x, prec);
    1726                 :            :   case t_POL:
    1727                 :          0 :     x = RgX_to_ser(x, precdl+2+1); /* +1 because deriv reduce the precision by 1 */
    1728                 :            :   case t_SER: /* FALL THROUGH */
    1729                 :          7 :     vx = varn(x);
    1730                 :          7 :     return gerepileupto(av, gdiv(deriv(eval(E, x),vx), deriv(x,vx)));
    1731                 :          0 :   default: pari_err_TYPE("formal derivation",x);
    1732                 :         28 :     return NULL; /*NOT REACHED*/
    1733                 :            :   }
    1734                 :            : }
    1735                 :            : 
    1736                 :            : GEN
    1737                 :         14 : derivnum0(GEN a, GEN code, long prec)
    1738                 :            : {
    1739                 :         14 :   EXPR_WRAP(code, derivfun (EXPR_ARG,a,prec));
    1740                 :            : }
    1741                 :            : 
    1742                 :            : struct deriv_data
    1743                 :            : {
    1744                 :            :   GEN code;
    1745                 :            :   GEN args;
    1746                 :            : };
    1747                 :            : 
    1748                 :         28 : static GEN deriv_eval(void *E, GEN x)
    1749                 :            : {
    1750                 :         28 :  struct deriv_data *data=(struct deriv_data *)E;
    1751                 :         28 :  gel(data->args,1)=x;
    1752                 :         28 :  return closure_callgenvec(data->code, data->args);
    1753                 :            : }
    1754                 :            : 
    1755                 :            : GEN
    1756                 :         14 : derivfun0(GEN code, GEN args, long prec)
    1757                 :            : {
    1758                 :            :   struct deriv_data E;
    1759                 :         14 :   E.code=code; E.args=args;
    1760                 :         14 :   return derivfun((void*)&E, deriv_eval, gel(args,1), prec);
    1761                 :            : }

Generated by: LCOV version 1.9