Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - language - forprime.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21353-12523aa) Lines: 370 440 84.1 %
Date: 2017-11-24 06:20:58 Functions: 32 37 86.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2016  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /**********************************************************************/
      18             : /***                                                                ***/
      19             : /***                     Public prime table                         ***/
      20             : /***                                                                ***/
      21             : /**********************************************************************/
      22             : 
      23             : static ulong _maxprime = 0;
      24             : static ulong diffptrlen;
      25             : 
      26             : /* Building/Rebuilding the diffptr table. The actual work is done by the
      27             :  * following two subroutines;  the user entry point is the function
      28             :  * initprimes() below.  initprimes1() is the old algorithm, called when
      29             :  * maxnum (size) is moderate. Must be called after pari_init_stack() )*/
      30             : static void
      31        1536 : initprimes1(ulong size, long *lenp, ulong *lastp, byteptr p1)
      32             : {
      33        1536 :   pari_sp av = avma;
      34             :   long k;
      35        1536 :   byteptr q, r, s, p = (byteptr)stack_calloc(size+2), fin = p + size;
      36             : 
      37       16896 :   for (r=q=p,k=1; r<=fin; )
      38             :   {
      39       21504 :     do { r+=k; k+=2; r+=k; } while (*++q);
      40       13824 :     for (s=r; s<=fin; s+=k) *s = 1;
      41             :   }
      42        1536 :   r = p1; *r++ = 2; *r++ = 1; /* 2 and 3 */
      43      192000 :   for (s=q=p+1; ; s=q)
      44             :   {
      45      542208 :     do q++; while (*q);
      46      192000 :     if (q > fin) break;
      47      190464 :     *r++ = (unsigned char) ((q-s) << 1);
      48      190464 :   }
      49        1536 :   *r++ = 0;
      50        1536 :   *lenp = r - p1;
      51        1536 :   *lastp = ((s - p) << 1) + 1;
      52        1536 :   avma = av;
      53        1536 : }
      54             : 
      55             : /*  Timing in ms (Athlon/850; reports 512K of secondary cache; looks
      56             :     like there is 64K of quickier cache too).
      57             : 
      58             :       arena|    30m     100m    300m    1000m    2000m  <-- primelimit
      59             :       =================================================
      60             :       16K       1.1053  1.1407  1.2589  1.4368   1.6086
      61             :       24K       1.0000  1.0625  1.1320  1.2443   1.3095
      62             :       32K       1.0000  1.0469  1.0761  1.1336   1.1776
      63             :       48K       1.0000  1.0000  1.0254  1.0445   1.0546
      64             :       50K       1.0000  1.0000  1.0152  1.0345   1.0464
      65             :       52K       1.0000  1.0000  1.0203  1.0273   1.0362
      66             :       54K       1.0000  1.0000  1.0812  1.0216   1.0281
      67             :       56K       1.0526  1.0000  1.0051  1.0144   1.0205
      68             :       58K       1.0000  1.0000  1.0000  1.0086   1.0123
      69             :       60K       0.9473  0.9844  1.0051  1.0014   1.0055
      70             :       62K       1.0000  0.9844  0.9949  0.9971   0.9993
      71             :       64K       1.0000  1.0000  1.0000  1.0000   1.0000
      72             :       66K       1.2632  1.2187  1.2183  1.2055   1.1953
      73             :       68K       1.4211  1.4844  1.4721  1.4425   1.4188
      74             :       70K       1.7368  1.7188  1.7107  1.6767   1.6421
      75             :       72K       1.9474  1.9531  1.9594  1.9023   1.8573
      76             :       74K       2.2105  2.1875  2.1827  2.1207   2.0650
      77             :       76K       2.4211  2.4219  2.4010  2.3305   2.2644
      78             :       78K       2.5789  2.6250  2.6091  2.5330   2.4571
      79             :       80K       2.8421  2.8125  2.8223  2.7213   2.6380
      80             :       84K       3.1053  3.1875  3.1776  3.0819   2.9802
      81             :       88K       3.5263  3.5312  3.5228  3.4124   3.2992
      82             :       92K       3.7895  3.8438  3.8375  3.7213   3.5971
      83             :       96K       4.0000  4.1093  4.1218  3.9986   3.9659
      84             :       112K      4.3684  4.5781  4.5787  4.4583   4.6115
      85             :       128K      4.7368  4.8750  4.9188  4.8075   4.8997
      86             :       192K      5.5263  5.7188  5.8020  5.6911   5.7064
      87             :       256K      6.0000  6.2187  6.3045  6.1954   6.1033
      88             :       384K      6.7368  6.9531  7.0405  6.9181   6.7912
      89             :       512K      7.3158  7.5156  7.6294  7.5000   7.4654
      90             :       768K      9.1579  9.4531  9.6395  9.5014   9.1075
      91             :       1024K    10.368  10.7497 10.9999 10.878   10.8201
      92             :       1536K    12.579  13.3124 13.7660 13.747   13.4739
      93             :       2048K    13.737  14.4839 15.0509 15.151   15.1282
      94             :       3076K    14.789  15.5780 16.2993 16.513   16.3365
      95             : 
      96             :     Now the same number relative to the model
      97             : 
      98             :     (1 + 0.36*sqrt(primelimit)/arena) * (arena <= 64 ? 1.05 : (arena-64)**0.38)
      99             : 
     100             :      [SLOW2_IN_ROOTS = 0.36, ALPHA = 0.38]
     101             : 
     102             :       arena|    30m     100m    300m    1000m    2000m  <-- primelimit
     103             :       =================================================
     104             :         16K    1.014    0.9835  0.9942  0.9889  1.004
     105             :         24K    0.9526   0.9758  0.9861  0.9942  0.981
     106             :         32K    0.971    0.9939  0.9884  0.9849  0.9806
     107             :         48K    0.9902   0.9825  0.996   0.9945  0.9885
     108             :         50K    0.9917   0.9853  0.9906  0.9926  0.9907
     109             :         52K    0.9932   0.9878  0.9999  0.9928  0.9903
     110             :         54K    0.9945   0.9902  1.064   0.9939  0.9913
     111             :         56K    1.048    0.9924  0.9925  0.993   0.9921
     112             :         58K    0.9969   0.9945  0.9909  0.9932  0.9918
     113             :         60K    0.9455   0.9809  0.9992  0.9915  0.9923
     114             :         62K    0.9991   0.9827  0.9921  0.9924  0.9929
     115             :         64K    1        1       1       1       1
     116             :         66K    1.02     0.9849  0.9857  0.9772  0.9704
     117             :         68K    0.8827   0.9232  0.9176  0.9025  0.8903
     118             :         70K    0.9255   0.9177  0.9162  0.9029  0.8881
     119             :         72K    0.9309   0.936   0.9429  0.9219  0.9052
     120             :         74K    0.9715   0.9644  0.967   0.9477  0.9292
     121             :         76K    0.9935   0.9975  0.9946  0.9751  0.9552
     122             :         78K    0.9987   1.021   1.021   1.003   0.9819
     123             :         80K    1.047    1.041   1.052   1.027   1.006
     124             :         84K    1.052    1.086   1.092   1.075   1.053
     125             :         88K    1.116    1.125   1.133   1.117   1.096
     126             :         92K    1.132    1.156   1.167   1.155   1.134
     127             :         96K    1.137    1.177   1.195   1.185   1.196
     128             :        112K    1.067    1.13    1.148   1.15    1.217
     129             :        128K    1.04     1.083   1.113   1.124   1.178
     130             :        192K    0.9368   0.985   1.025   1.051   1.095
     131             :        256K    0.8741   0.9224  0.9619  0.995   1.024
     132             :        384K    0.8103   0.8533  0.8917  0.9282  0.9568
     133             :        512K    0.7753   0.8135  0.8537  0.892   0.935
     134             :        768K    0.8184   0.8638  0.9121  0.9586  0.9705
     135             :       1024K    0.8241   0.8741  0.927   0.979   1.03
     136             :       1536K    0.8505   0.9212  0.9882  1.056   1.096
     137             :       2048K    0.8294   0.8954  0.9655  1.041   1.102
     138             : 
     139             : */
     140             : 
     141             : #ifndef SLOW2_IN_ROOTS
     142             :   /* SLOW2_IN_ROOTS below 3: some slowdown starts to be noticable
     143             :    * when things fit into the cache on Sparc.
     144             :    * The choice of 2.6 gives a slowdown of 1-2% on UltraSparcII,
     145             :    * but makes calculations for "maximum" of 436273009
     146             :    * fit into 256K cache (still common for some architectures).
     147             :    *
     148             :    * One may change it when small caches become uncommon, but the gain
     149             :    * is not going to be very noticable... */
     150             : #  ifdef i386           /* gcc defines this? */
     151             : #    define SLOW2_IN_ROOTS      0.36
     152             : #  else
     153             : #    define SLOW2_IN_ROOTS      2.6
     154             : #  endif
     155             : #endif
     156             : #ifndef CACHE_ARENA
     157             : #  ifdef i386           /* gcc defines this? */
     158             :    /* Due to smaller SLOW2_IN_ROOTS, smaller arena is OK; fit L1 cache */
     159             : #    define CACHE_ARENA (63 * 1024UL) /* No slowdown even with 64K L1 cache */
     160             : #  else
     161             : #    define CACHE_ARENA (200 * 1024UL) /* No slowdown even with 256K L2 cache */
     162             : #  endif
     163             : #endif
     164             : 
     165             : #define CACHE_ALPHA     (0.38)          /* Cache performance model parameter */
     166             : #define CACHE_CUTOFF    (0.018)         /* Cache performance not smooth here */
     167             : 
     168             : static double slow2_in_roots = SLOW2_IN_ROOTS;
     169             : 
     170             : typedef struct {
     171             :     ulong arena;
     172             :     double power;
     173             :     double cutoff;
     174             : } cache_model_t;
     175             : 
     176             : static cache_model_t cache_model = { CACHE_ARENA, CACHE_ALPHA, CACHE_CUTOFF };
     177             : 
     178             : /* Assume that some calculation requires a chunk of memory to be
     179             :    accessed often in more or less random fashion (as in sieving).
     180             :    Assume that the calculation can be done in steps by subdividing the
     181             :    chunk into smaller subchunks (arenas) and treating them
     182             :    separately.  Assume that the overhead of subdivision is equivalent
     183             :    to the number of arenas.
     184             : 
     185             :    Find an optimal size of the arena taking into account the overhead
     186             :    of subdivision, and the overhead of arena not fitting into the
     187             :    cache.  Assume that arenas of size slow2_in_roots slows down the
     188             :    calculation 2x (comparing to very big arenas; when cache hits do
     189             :    not matter).  Since cache performance varies wildly with
     190             :    architecture, load, and wheather (especially with cache coloring
     191             :    enabled), use an idealized cache model based on benchmarks above.
     192             : 
     193             :    Assume that an independent region of FIXED_TO_CACHE bytes is accessed
     194             :    very often concurrently with the arena access.
     195             :  */
     196             : static ulong
     197        1536 : good_arena_size(ulong slow2_size, ulong total, ulong fixed_to_cache,
     198             :                 cache_model_t *cache_model)
     199             : {
     200        1536 :   ulong asize, cache_arena = cache_model->arena;
     201             :   double Xmin, Xmax, A, B, C1, C2, D, V;
     202        1536 :   double alpha = cache_model->power, cut_off = cache_model->cutoff;
     203             : 
     204             :   /* Estimated relative slowdown,
     205             :      with overhead = max((fixed_to_cache+arena)/cache_arena - 1, 0):
     206             : 
     207             :      1 + slow2_size/arena due to initialization overhead;
     208             : 
     209             :      max(1, 4.63 * overhead^0.38 ) due to footprint > cache size.
     210             : 
     211             :      [The latter is hard to substantiate theoretically, but this
     212             :      function describes benchmarks pretty close; it does not hurt that
     213             :      one can minimize it explicitly too ;-).  The switch between
     214             :      different choices of max() happens when overhead=0.018.]
     215             : 
     216             :      Thus the problem is minimizing (1 + slow2_size/arena)*overhead**0.29.
     217             :      This boils down to F=((X+A)/(X+B))X^alpha, X=overhead,
     218             :      B = (1 - fixed_to_cache/cache_arena), A = B + slow2_size/cache_arena,
     219             :      alpha = 0.38, and X>=0.018, X>-B.
     220             : 
     221             :      We need to find the rightmost root of (X+A)*(X+B) - alpha(A-B)X to the
     222             :      right of 0.018 (if such exists and is below Xmax).  Then we manually
     223             :      check the remaining region [0, 0.018].
     224             : 
     225             :      Since we cannot trust the purely-experimental cache-hit slowdown
     226             :      function, as a sanity check always prefer fitting into the
     227             :      cache (or "almost fitting") if F-law predicts that the larger
     228             :      value of the arena provides less than 10% speedup.
     229             :    */
     230             : 
     231             :   /* The simplest case: we fit into cache */
     232        1536 :   asize = cache_arena - fixed_to_cache;
     233        1536 :   if (total <= asize) return total;
     234             :   /* The simple case: fitting into cache doesn't slow us down more than 10% */
     235        1536 :   if (asize > 10 * slow2_size) return asize;
     236             :   /* Slowdown of not fitting into cache is significant.  Try to optimize.
     237             :      Do not be afraid to spend some time on optimization - in trivial
     238             :      cases we do not reach this point; any gain we get should
     239             :      compensate the time spent on optimization.  */
     240             : 
     241           0 :   B = (1 - ((double)fixed_to_cache)/cache_arena);
     242           0 :   A = B + ((double)slow2_size)/cache_arena;
     243           0 :   C2 = A*B;
     244           0 :   C1 = (A + B - 1/alpha*(A - B))/2;
     245           0 :   D = C1*C1 - C2;
     246           0 :   if (D > 0)
     247           0 :     V = cut_off*cut_off + 2*C1*cut_off + C2; /* Value at CUT_OFF */
     248             :   else
     249           0 :     V = 0; /* Peacify the warning */
     250           0 :   Xmin = cut_off;
     251           0 :   Xmax = ((double)total - fixed_to_cache)/cache_arena; /* Two candidates */
     252             : 
     253           0 :   if ( D <= 0 || (V >= 0 && C1 + cut_off >= 0) ) /* slowdown increasing */
     254           0 :     Xmax = cut_off; /* Only one candidate */
     255           0 :   else if (V >= 0 && /* slowdown concave down */
     256           0 :            ((Xmax + C1) <= 0 || (Xmax*Xmax + 2*C1*Xmax + C2) <= 0))
     257             :       /* DO NOTHING */;  /* Keep both candidates */
     258           0 :   else if (V <= 0 && (Xmax*Xmax + 2*C1*Xmax + C2) <= 0) /*slowdown decreasing*/
     259           0 :       Xmin = cut_off; /* Only one candidate */
     260             :   else /* Now we know: 2 roots, the largest is in CUT_OFF..Xmax */
     261           0 :       Xmax = sqrt(D) - C1;
     262           0 :   if (Xmax != Xmin) { /* Xmin == CUT_OFF; Check which one is better */
     263           0 :     double v1 = (cut_off + A)/(cut_off + B);
     264           0 :     double v2 = 2.33 * (Xmax + A)/(Xmax + B) * pow(Xmax, alpha);
     265             : 
     266           0 :     if (1.1 * v2 >= v1) /* Prefer fitting into the cache if slowdown < 10% */
     267           0 :       V = v1;
     268             :     else
     269           0 :     { Xmin = Xmax; V = v2; }
     270           0 :   } else if (B > 0) /* We need V */
     271           0 :     V = 2.33 * (Xmin + A)/(Xmin + B) * pow(Xmin, alpha);
     272           0 :   if (B > 0 && 1.1 * V > A/B)  /* Now Xmin is the minumum.  Compare with 0 */
     273           0 :     Xmin = 0;
     274             : 
     275           0 :   asize = (ulong)((1 + Xmin)*cache_arena - fixed_to_cache);
     276           0 :   if (asize > total) asize = total; /* May happen due to approximations */
     277           0 :   return asize;
     278             : }
     279             : 
     280             : /* Use as in
     281             :     install(set_optimize,lLDG)          \\ Through some M too?
     282             :     set_optimize(2,1) \\ disable dependence on limit
     283             :     \\ 1: how much cache usable, 2: slowdown of setup, 3: alpha, 4: cutoff
     284             :     \\ 2,3,4 are in units of 0.001
     285             : 
     286             :     { time_primes_arena(ar,limit) =     \\ ar = arena size in K
     287             :         set_optimize(1,floor(ar*1024));
     288             :         default(primelimit, 200 000);   \\ 100000 results in *larger* malloc()!
     289             :         gettime;
     290             :         default(primelimit, floor(limit));
     291             :         if(ar >= 1, ar=floor(ar));
     292             :         print("arena "ar"K => "gettime"ms");
     293             :     }
     294             : */
     295             : long
     296           0 : set_optimize(long what, GEN g)
     297             : {
     298           0 :   long ret = 0;
     299             : 
     300           0 :   switch (what) {
     301             :   case 1:
     302           0 :     ret = (long)cache_model.arena;
     303           0 :     break;
     304             :   case 2:
     305           0 :     ret = (long)(slow2_in_roots * 1000);
     306           0 :     break;
     307             :   case 3:
     308           0 :     ret = (long)(cache_model.power * 1000);
     309           0 :     break;
     310             :   case 4:
     311           0 :     ret = (long)(cache_model.cutoff * 1000);
     312           0 :     break;
     313             :   default:
     314           0 :     pari_err_BUG("set_optimize");
     315           0 :     break;
     316             :   }
     317           0 :   if (g != NULL) {
     318           0 :     ulong val = itou(g);
     319             : 
     320           0 :     switch (what) {
     321           0 :     case 1: cache_model.arena = val; break;
     322           0 :     case 2: slow2_in_roots     = (double)val / 1000.; break;
     323           0 :     case 3: cache_model.power  = (double)val / 1000.; break;
     324           0 :     case 4: cache_model.cutoff = (double)val / 1000.; break;
     325             :     }
     326             :   }
     327           0 :   return ret;
     328             : }
     329             : 
     330             : /* s is odd; prime differences (starting from 5-3=2) start at known_primes[2],
     331             :   terminated by a 0 byte. Checks n odd numbers starting at 'start', setting
     332             :   bytes starting at data to 0 (composite) or 1 (prime) */
     333             : static void
     334        3496 : sieve_chunk(byteptr known_primes, ulong s, byteptr data, ulong n)
     335             : {
     336        3496 :   ulong p, cnt = n-1, start = s, delta = 1;
     337             :   byteptr q;
     338             : 
     339        3496 :   memset(data, 0, n);
     340        3496 :   start >>= 1;  /* (start - 1)/2 */
     341        3496 :   start += n; /* Corresponds to the end */
     342             :   /* data corresponds to start, q runs over primediffs */
     343      409076 :   for (q = known_primes + 1, p = 3; delta; delta = *++q, p += delta)
     344             :   { /* first odd number >= start > p and divisible by p
     345             :        = last odd number <= start + 2p - 2 and 0 (mod p)
     346             :        = p + last number <= start + p - 2 and 0 (mod 2p)
     347             :        = p + start+p-2 - (start+p-2) % 2p
     348             :        = start + 2(p - 1 - ((start-1)/2 + (p-1)/2) % p). */
     349      405580 :     long off = cnt - ((start+(p>>1)) % p);
     350      405580 :     while (off >= 0) { data[off] = 1; off -= p; }
     351             :   }
     352        3496 : }
     353             : 
     354             : /* assume maxnum <= 436273289 < 2^29 */
     355             : static void
     356        1536 : initprimes0(ulong maxnum, long *lenp, ulong *lastp, byteptr p1)
     357             : {
     358        1536 :   pari_sp av = avma, bot = pari_mainstack->bot;
     359             :   long alloced, psize;
     360             :   byteptr q, end, p, end1, plast, curdiff;
     361             :   ulong last, remains, curlow, rootnum, asize;
     362             :   ulong prime_above;
     363             :   byteptr p_prime_above;
     364             : 
     365        1536 :   maxnum |= 1; /* make it odd. */
     366             :   /* base case */
     367        3072 :   if (maxnum < 1ul<<17) { initprimes1(maxnum>>1, lenp, lastp, p1); return; }
     368             : 
     369             :   /* Checked to be enough up to 40e6, attained at 155893 */
     370        1536 :   rootnum = usqrt(maxnum) | 1;
     371        1536 :   initprimes1(rootnum>>1, &psize, &last, p1);
     372        1536 :   end1 = p1 + psize - 1;
     373        1536 :   remains = (maxnum - last) >> 1; /* number of odd numbers to check */
     374             : 
     375             :   /* we access primes array of psize too; but we access it consecutively,
     376             :    * thus we do not include it in fixed_to_cache */
     377        1536 :   asize = good_arena_size((ulong)(rootnum * slow2_in_roots), remains+1, 0,
     378             :                           &cache_model) - 1;
     379             :   /* enough room on the stack ? */
     380        1536 :   alloced = (((byteptr)avma) <= ((byteptr)bot) + asize);
     381        1536 :   if (alloced)
     382           0 :     p = (byteptr)pari_malloc(asize+1);
     383             :   else
     384        1536 :     p = (byteptr)stack_malloc(asize+1);
     385        1536 :   end = p + asize; /* the 0 sentinel goes at end. */
     386        1536 :   curlow = last + 2; /* First candidate: know primes up to last (odd). */
     387        1536 :   curdiff = end1;
     388             : 
     389             :   /* During each iteration p..end-1 represents a range of odd
     390             :      numbers.  plast is a pointer which represents the last prime seen,
     391             :      it may point before p..end-1. */
     392        1536 :   plast = p - 1;
     393        1536 :   p_prime_above = p1 + 2;
     394        1536 :   prime_above = 3;
     395        6568 :   while (remains)
     396             :   { /* cycle over arenas; performance not crucial */
     397             :     unsigned char was_delta;
     398        3496 :     if (asize > remains) { asize = remains; end = p + asize; }
     399             :     /* Fake the upper limit appropriate for the given arena */
     400      198992 :     while (prime_above*prime_above <= curlow + (asize << 1) && *p_prime_above)
     401      192000 :       prime_above += *p_prime_above++;
     402        3496 :     was_delta = *p_prime_above;
     403        3496 :     *p_prime_above = 0; /* sentinel for sieve_chunk */
     404        3496 :     sieve_chunk(p1, curlow, p, asize);
     405        3496 :     *p_prime_above = was_delta; /* restore */
     406             : 
     407        3496 :     p[asize] = 0; /* sentinel */
     408    63612328 :     for (q = p; ; plast = q++)
     409             :     { /* q runs over addresses corresponding to primes */
     410    63612328 :       while (*q) q++; /* use sentinel at end */
     411    63612328 :       if (q >= end) break;
     412    63608832 :       *curdiff++ = (unsigned char)(q-plast) << 1; /* < 255 for q < 436273291 */
     413    63608832 :     }
     414        3496 :     plast -= asize;
     415        3496 :     remains -= asize;
     416        3496 :     curlow += (asize<<1);
     417             :   }
     418        1536 :   last = curlow - ((p - plast) << 1);
     419        1536 :   *curdiff++ = 0; /* sentinel */
     420        1536 :   *lenp = curdiff - p1;
     421        1536 :   *lastp = last;
     422        1536 :   if (alloced) pari_free(p); else avma = av;
     423             : }
     424             : 
     425             : ulong
     426    36588128 : maxprime(void) { return diffptr ? _maxprime : 0; }
     427             : 
     428             : void
     429           0 : maxprime_check(ulong c) { if (_maxprime < c) pari_err_MAXPRIME(c); }
     430             : 
     431             : /* We ensure 65302 <= maxnum <= 436273289: the LHS ensures modular function
     432             :  * have enough fast primes to work, the RHS ensures that p_{n+1} - p_n < 255
     433             :  * (N.B. RHS would be incorrect since initprimes0 would make it odd, thereby
     434             :  * increasing it by 1) */
     435             : byteptr
     436        1536 : initprimes(ulong maxnum, long *lenp, ulong *lastp)
     437             : {
     438             :   byteptr t;
     439        1536 :   if (maxnum < 65537)
     440           0 :     maxnum = 65537;
     441        1536 :   else if (maxnum > 436273289)
     442           0 :     maxnum = 436273289;
     443        1536 :   t = (byteptr)pari_malloc((size_t) (1.09 * maxnum/log((double)maxnum)) + 146);
     444        1536 :   initprimes0(maxnum, lenp, lastp, t);
     445        1536 :   return (byteptr)pari_realloc(t, *lenp);
     446             : }
     447             : 
     448             : void
     449        1536 : initprimetable(ulong maxnum)
     450             : {
     451             :   long len;
     452             :   ulong last;
     453        1536 :   byteptr p = initprimes(maxnum, &len, &last), old = diffptr;
     454        1536 :   diffptrlen = minss(diffptrlen, len);
     455        1536 :   _maxprime  = minss(_maxprime,last); /*Protect against ^C*/
     456        1536 :   diffptr = p; diffptrlen = len; _maxprime = last;
     457        1536 :   if (old) free(old);
     458        1536 : }
     459             : 
     460             : /* all init_primepointer_xx routines set *ptr to the corresponding place
     461             :  * in prime table */
     462             : /* smallest p >= a */
     463             : ulong
     464           0 : init_primepointer_geq(ulong a, byteptr *pd)
     465             : {
     466             :   ulong n, p;
     467           0 :   prime_table_next_p(a, pd, &p, &n);
     468           0 :   return p;
     469             : }
     470             : /* largest p < a */
     471             : ulong
     472    14490108 : init_primepointer_lt(ulong a, byteptr *pd)
     473             : {
     474             :   ulong n, p;
     475    14490108 :   prime_table_next_p(a, pd, &p, &n);
     476    14490120 :   PREC_PRIME_VIADIFF(p, *pd);
     477    14490120 :   return p;
     478             : }
     479             : /* largest p <= a */
     480             : ulong
     481           0 : init_primepointer_leq(ulong a, byteptr *pd)
     482             : {
     483             :   ulong n, p;
     484           0 :   prime_table_next_p(a, pd, &p, &n);
     485           0 :   if (p != a) PREC_PRIME_VIADIFF(p, *pd);
     486           0 :   return p;
     487             : }
     488             : /* smallest p > a */
     489             : ulong
     490           0 : init_primepointer_gt(ulong a, byteptr *pd)
     491             : {
     492             :   ulong n, p;
     493           0 :   prime_table_next_p(a, pd, &p, &n);
     494           0 :   if (p == a) NEXT_PRIME_VIADIFF(p, *pd);
     495           0 :   return p;
     496             : }
     497             : 
     498             : /**********************************************************************/
     499             : /***                                                                ***/
     500             : /***                     forprime                                   ***/
     501             : /***                                                                ***/
     502             : /**********************************************************************/
     503             : 
     504             : /* return good chunk size for sieve, 16 | chunk + 2 */
     505             : static ulong
     506     1525136 : optimize_chunk(ulong a, ulong b)
     507             : {
     508             :   /* TODO: Optimize size (surely < 512k to stay in L2 cache, but not so large
     509             :    * as to force recalculating too often). */
     510     1525136 :   ulong chunk = 0x80000UL;
     511     1525136 :   ulong tmp = (b - a) / chunk + 1;
     512             : 
     513     1525136 :   if (tmp == 1)
     514          28 :     chunk = b - a + 16;
     515             :   else
     516     1525108 :     chunk = (b - a) / tmp + 15;
     517             :   /* ensure 16 | chunk + 2 */
     518     1525136 :   return (((chunk + 2)>>4)<<4) - 2;
     519             : }
     520             : static void
     521     1525136 : sieve_init(forprime_t *T, ulong a, ulong b)
     522             : {
     523     1525136 :   T->sieveb = b;
     524     1525136 :   T->chunk = optimize_chunk(a, b);
     525     1525136 :   T->isieve = (unsigned char*)stack_malloc(((T->chunk+2) >> 4) + 1);
     526     1525136 :   T->cache[0] = 0;
     527             :   /* >> 1 [only odds] + 3 [convert from bits to bytes] */
     528     1525136 :   T->a = a;
     529     1525136 :   T->end = minuu(a + T->chunk, b);
     530     1525136 :   T->pos = T->maxpos = 0;
     531     1525136 : }
     532             : 
     533             : enum {PRST_none, PRST_diffptr, PRST_sieve, PRST_unextprime, PRST_nextprime};
     534             : 
     535             : static void
     536    17371278 : u_forprime_set_prime_table(forprime_t *T, ulong a)
     537             : {
     538    17371278 :   T->strategy = PRST_diffptr;
     539    17371278 :   if (a < 3)
     540             :   {
     541     2881109 :     T->p = 0;
     542     2881109 :     T->d = diffptr;
     543             :   }
     544             :   else
     545    14490169 :     T->p = init_primepointer_lt(a, &T->d);
     546    17371226 : }
     547             : 
     548             : /* Set p so that p + q the smallest integer = c (mod q) and > original p.
     549             :  * Assume 0 < c < q. Set p = 0 on overflow */
     550             : static void
     551         721 : arith_set(forprime_t *T)
     552             : {
     553         721 :   ulong r = T->p % T->q; /* 0 <= r <= min(p, q-1) */
     554         721 :   pari_sp av = avma;
     555         721 :   GEN d = adduu(T->p - r, T->c);
     556         721 :   if (T->c > r) d = subiu(d, T->q);
     557             :   /* d = c mod q,  d = c > r? p-r+c-q: p-r+c, so that
     558             :    *  d <= p  and  d+q = c>r? p-r+c  : p-r+c+q > p */
     559         721 :   T->p = itou_or_0(d); avma = av; /* d = 0 is impossible */
     560         721 : }
     561             : 
     562             : /* run through primes in arithmetic progression = c (mod q).
     563             :  * Assume (c,q)=1, 0 <= c < q */
     564             : static int
     565    22891836 : u_forprime_sieve_arith_init(forprime_t *T, struct pari_sieve *psieve,
     566             :                             ulong a, ulong b, ulong c, ulong q)
     567             : {
     568             :   ulong maxp, maxp2;
     569    22891836 :   if (!odd(b) && b > 2) b--;
     570    22891892 :   if (a > b || b < 2)
     571             :   {
     572      849505 :     T->strategy = PRST_diffptr; /* paranoia */
     573      849505 :     T->p = 0; /* empty */
     574      849505 :     T->b = 0; /* empty */
     575      849505 :     T->d = diffptr;
     576      849505 :     return 0;
     577             :   }
     578    22042387 :   maxp = maxprime();
     579    22042286 :   if (q != 1 && c != 2 && odd(q)) {
     580             :     /* only allow *odd* primes. If c = 2, then p = 2 must be included :-( */
     581      261975 :     if (!odd(c)) c += q;
     582      261975 :     q <<= 1;
     583             :   }
     584    22042286 :   T->q = q;
     585    22042286 :   T->c = c;
     586    22042286 :   T->strategy = PRST_none; /* unknown */
     587    22042286 :   T->psieve = psieve; /* unused for now */
     588    22042286 :   T->isieve = NULL; /* unused for now */
     589    22042286 :   T->b = b;
     590    22042286 :   if (maxp >= b) { /* [a,b] \subset prime table */
     591    16265162 :     u_forprime_set_prime_table(T, a);
     592    16265114 :     return 1;
     593             :   }
     594             :   /* b > maxp */
     595     5777124 :   if (a >= maxp)
     596             :   {
     597     4670998 :     T->p = a - 1;
     598     4670998 :     if (T->q > 1) arith_set(T);
     599             :   }
     600             :   else
     601     1106126 :     u_forprime_set_prime_table(T, a);
     602             : 
     603     5777124 :   maxp2 = (maxp & HIGHMASK)? 0 : maxp*maxp;
     604             :   /* FIXME: should sieve as well if q != 1, adapt sieve code */
     605     5777124 :   if (q != 1 || (maxp2 && maxp2 <= a)
     606     1540566 :              || T->b - maxuu(a,maxp) < maxp / expu(b))
     607     4251988 :   { if (T->strategy==PRST_none) T->strategy = PRST_unextprime; }
     608             :   else
     609             :   { /* worth sieving */
     610             : #ifdef LONG_IS_64BIT
     611      735332 :     const ulong UPRIME_MAX = 18446744073709551557UL;
     612             : #else
     613      789804 :     const ulong UPRIME_MAX = 4294967291UL;
     614             : #endif
     615             :     ulong sieveb;
     616     1525136 :     if (b > UPRIME_MAX) b = UPRIME_MAX;
     617     1525136 :     sieveb = b;
     618     1525136 :     if (maxp2 && maxp2 < b) sieveb = maxp2;
     619     1525136 :     if (T->strategy==PRST_none) T->strategy = PRST_sieve;
     620     1525136 :     sieve_init(T, maxuu(maxp+2, a), sieveb);
     621             :   }
     622     5777124 :   return 1;
     623             : }
     624             : 
     625             : int
     626    22196352 : u_forprime_arith_init(forprime_t *T, ulong a, ulong b, ulong c, ulong q)
     627    22196352 : { return u_forprime_sieve_arith_init(T, NULL, a, b, c, q); }
     628             : 
     629             : /* will run through primes in [a,b] */
     630             : int
     631    21934103 : u_forprime_init(forprime_t *T, ulong a, ulong b)
     632    21934103 : { return u_forprime_arith_init(T, a,b, 0,1); }
     633             : 
     634             : /* will run through primes in [a,b] */
     635             : static int
     636      695375 : u_forprime_sieve_init(forprime_t *T, struct pari_sieve *s, ulong b)
     637      695375 : { return u_forprime_sieve_arith_init(T, s, s->start, b, s->c, s->q); }
     638             : 
     639             : /* now only run through primes <= c; assume c <= b above */
     640             : void
     641          63 : u_forprime_restrict(forprime_t *T, ulong c) { T->b = c; }
     642             : 
     643             : /* b = NULL: loop forever */
     644             : int
     645         562 : forprime_init(forprime_t *T, GEN a, GEN b)
     646             : {
     647             :   long lb;
     648         562 :   a = gceil(a); if (typ(a) != t_INT) pari_err_TYPE("forprime_init",a);
     649         562 :   if (signe(a) <= 0) a = gen_1;
     650         562 :   if (b && typ(b) != t_INFINITY)
     651             :   {
     652         457 :     b = gfloor(b);
     653         457 :     if (typ(b) != t_INT) pari_err_TYPE("forprime_init",b);
     654         457 :     if (signe(b) < 0 || cmpii(a,b) > 0)
     655             :     {
     656           7 :       T->strategy = PRST_nextprime; /* paranoia */
     657           7 :       T->bb = gen_0;
     658           7 :       T->pp = gen_0;
     659           7 :       return 0;
     660             :     }
     661         450 :     lb = lgefint(b);
     662         450 :     T->bb = b;
     663             :   }
     664         105 :   else if (!b || inf_get_sign(b) > 0)
     665             :   {
     666         105 :     lb = lgefint(a) + 4;
     667         105 :     T->bb = NULL;
     668             :   }
     669             :   else /* b == -oo */
     670             :   {
     671           0 :     T->strategy = PRST_nextprime; /* paranoia */
     672           0 :     T->bb = gen_0;
     673           0 :     T->pp = gen_0;
     674           0 :     return 0;
     675             :   }
     676         555 :   T->pp = cgeti(lb);
     677             :   /* a, b are positive integers, a <= b */
     678         555 :   if (lgefint(a) == 3) /* lb == 3 implies b != NULL */
     679         440 :     return u_forprime_init(T, uel(a,2), lb == 3? uel(b,2): ULONG_MAX);
     680         115 :   T->strategy = PRST_nextprime;
     681         115 :   affii(subiu(a,1), T->pp);
     682         115 :   return 1;
     683             : }
     684             : 
     685             : /* assume a <= b <= maxprime()^2, a,b odd, sieve[n] corresponds to
     686             :  *   a+16*n, a+16*n+2, ..., a+16*n+14 (bits 0 to 7)
     687             :  * maxpos = index of last sieve cell.
     688             :  * b-a+2 must be divisible by 16 for use by u_forprime_next */
     689             : static void
     690        4602 : sieve_block(ulong a, ulong b, ulong maxpos, unsigned char* sieve)
     691             : {
     692        4602 :   ulong p = 2, lim = usqrt(b), sz = (b-a) >> 1;
     693        4602 :   byteptr d = diffptr+1;
     694        4602 :   (void)memset(sieve, 0, maxpos+1);
     695             :   for (;;)
     696             :   { /* p is odd */
     697             :     ulong k, r;
     698    15010244 :     NEXT_PRIME_VIADIFF(p, d); /* starts at p = 3 */
     699    15010244 :     if (p > lim) break;
     700             : 
     701             :     /* solve a + 2k = 0 (mod p) */
     702    15005642 :     r = a % p;
     703    15005642 :     if (r == 0)
     704        7273 :       k = 0;
     705             :     else
     706             :     {
     707    14998369 :       k = p - r;
     708    14998369 :       if (odd(k)) k += p;
     709    14998369 :       k >>= 1;
     710             :     }
     711             :     /* m = a + 2k is the smallest odd m >= a, p | m */
     712             :     /* position n (corresponds to a+2n) is sieve[n>>3], bit n&7 */
     713    15005642 :     while (k <= sz) { sieve[k>>3] |= 1 << (k&7); k += p; /* 2k += 2p */ }
     714    15005642 :   }
     715        4602 : }
     716             : 
     717             : static void
     718        1536 : pari_sieve_init(struct pari_sieve *s, ulong a, ulong b)
     719             : {
     720        1536 :   ulong maxpos= (b - a) >> 4;
     721        1536 :   s->start = a; s->end = b;
     722        1536 :   s->sieve = (unsigned char*) pari_malloc(maxpos+1);
     723        1536 :   s->c = 0; s->q = 1;
     724        1536 :   sieve_block(a, b, maxpos, s->sieve);
     725        1536 :   s->maxpos = maxpos; /* must be last in case of SIGINT */
     726        1536 : }
     727             : 
     728             : static struct pari_sieve pari_sieve_modular;
     729             : 
     730             : #ifdef LONG_IS_64BIT
     731             : #define PARI_MODULAR_BASE ((1UL<<((BITS_IN_LONG-2)>>1))+1)
     732             : #else
     733             : #define PARI_MODULAR_BASE ((1UL<<(BITS_IN_LONG-1))+1)
     734             : #endif
     735             : 
     736             : void
     737        1536 : pari_init_primes(ulong maxprime)
     738             : {
     739        1536 :   ulong a = PARI_MODULAR_BASE, b = a + (1UL<<20)-2;
     740        1536 :   initprimetable(maxprime);
     741        1536 :   pari_sieve_init(&pari_sieve_modular, a, b);
     742        1536 : }
     743             : 
     744             : void
     745        1536 : pari_close_primes(void)
     746             : {
     747        1536 :   pari_free(diffptr);
     748        1536 :   pari_free(pari_sieve_modular.sieve);
     749        1536 : }
     750             : 
     751             : void
     752       33668 : init_modular_small(forprime_t *S)
     753             : {
     754             : #ifdef LONG_IS_64BIT
     755       28860 :   u_forprime_sieve_init(S, &pari_sieve_modular, ULONG_MAX);
     756             : #else
     757        4808 :   ulong a = (1UL<<((BITS_IN_LONG-2)>>1))+1;
     758        4808 :   u_forprime_init(S, a, ULONG_MAX);
     759             : #endif
     760       33668 : }
     761             : 
     762             : void
     763     4640807 : init_modular_big(forprime_t *S)
     764             : {
     765             : #ifdef LONG_IS_64BIT
     766     3974292 :   ulong a = (1UL<<(BITS_IN_LONG-1))+1;
     767     3974292 :   u_forprime_init(S, a, ULONG_MAX);
     768             : #else
     769      666515 :   u_forprime_sieve_init(S, &pari_sieve_modular, ULONG_MAX);
     770             : #endif
     771     4640807 : }
     772             : 
     773             : void
     774       99915 : init_modular(forprime_t *S)
     775             : {
     776       99915 :   u_forprime_init(S, 27449, ULONG_MAX);
     777       99915 : }
     778             : 
     779             : /* T->cache is a 0-terminated list of primes, return the first one and
     780             :  * remove it from list. Most of the time the list contains a single prime */
     781             : static ulong
     782    86002691 : shift_cache(forprime_t *T)
     783             : {
     784             :   long i;
     785    86002691 :   T->p = T->cache[0];
     786   114953651 :   for (i = 1;; i++)  /* remove one prime from cache */
     787   114953651 :     if (! (T->cache[i-1] = T->cache[i]) ) break;
     788    28950960 :   return T->p;
     789             : }
     790             : 
     791             : ulong
     792   420138619 : u_forprime_next(forprime_t *T)
     793             : {
     794   420138619 :   if (T->strategy == PRST_diffptr)
     795             :   {
     796             :     for(;;)
     797             :     {
     798   479991735 :       if (!*(T->d))
     799             :       {
     800        5650 :         T->strategy = T->isieve? PRST_sieve: PRST_unextprime;
     801        5650 :         if (T->q != 1) { arith_set(T); if (!T->p) return 0; }
     802             :         /* T->p possibly not a prime if q != 1 */
     803        5650 :         break;
     804             :       }
     805             :       else
     806             :       {
     807   479986085 :         NEXT_PRIME_VIADIFF(T->p, T->d);
     808   479986085 :         if (T->p > T->b) return 0;
     809   477927678 :         if (T->q == 1 || T->p % T->q == T->c) return T->p;
     810             :       }
     811   149992173 :     }
     812             :   }
     813    90144707 :   if (T->strategy == PRST_sieve)
     814             :   {
     815             :     ulong n;
     816    86002866 :     if (T->cache[0]) return shift_cache(T);
     817             : NEXT_CHUNK:
     818    61352776 :     if (T->psieve)
     819             :     {
     820      695375 :       T->sieve = T->psieve->sieve;
     821      695375 :       T->end = T->psieve->end;
     822      695375 :       if (T->end > T->sieveb) T->end = T->sieveb;
     823      695375 :       T->maxpos = T->psieve->maxpos;
     824      695375 :       T->pos = 0;
     825      695375 :       T->psieve = NULL;
     826             :     }
     827    96776164 :     for (n = T->pos; n < T->maxpos; n++)
     828    96771065 :       if (T->sieve[n] != 0xFF)
     829             :       {
     830    61347677 :         unsigned char mask = T->sieve[n];
     831    61347677 :         ulong p = T->a + (n<<4);
     832    61347677 :         long i = 0;
     833    61347677 :         T->pos = n;
     834    61347677 :         if (!(mask &  1)) T->cache[i++] = p;
     835    61347677 :         if (!(mask &  2)) T->cache[i++] = p+2;
     836    61347677 :         if (!(mask &  4)) T->cache[i++] = p+4;
     837    61347677 :         if (!(mask &  8)) T->cache[i++] = p+6;
     838    61347677 :         if (!(mask & 16)) T->cache[i++] = p+8;
     839    61347677 :         if (!(mask & 32)) T->cache[i++] = p+10;
     840    61347677 :         if (!(mask & 64)) T->cache[i++] = p+12;
     841    61347677 :         if (!(mask &128)) T->cache[i++] = p+14;
     842    61347677 :         T->cache[i] = 0;
     843    61347677 :         T->pos = n+1;
     844    61347677 :         return shift_cache(T);
     845             :       }
     846             :     /* n = T->maxpos, last cell: check p <= b */
     847        5099 :     if (T->maxpos && n == T->maxpos && T->sieve[n] != 0xFF)
     848             :     {
     849        1935 :       unsigned char mask = T->sieve[n];
     850        1935 :       ulong p = T->a + (n<<4);
     851        1935 :       long i = 0;
     852        1935 :       T->pos = n;
     853        1935 :       if (!(mask &  1) && p <= T->sieveb) T->cache[i++] = p;
     854        1935 :       if (!(mask &  2) && p <= T->sieveb-2) T->cache[i++] = p+2;
     855        1935 :       if (!(mask &  4) && p <= T->sieveb-4) T->cache[i++] = p+4;
     856        1935 :       if (!(mask &  8) && p <= T->sieveb-6) T->cache[i++] = p+6;
     857        1935 :       if (!(mask & 16) && p <= T->sieveb-8) T->cache[i++] = p+8;
     858        1935 :       if (!(mask & 32) && p <= T->sieveb-10) T->cache[i++] = p+10;
     859        1935 :       if (!(mask & 64) && p <= T->sieveb-12) T->cache[i++] = p+12;
     860        1935 :       if (!(mask &128) && p <= T->sieveb-14) T->cache[i++] = p+14;
     861        1935 :       if (i)
     862             :       {
     863        1858 :         T->cache[i] = 0;
     864        1858 :         T->pos = n+1;
     865        1858 :         return shift_cache(T);
     866             :       }
     867             :     }
     868             : 
     869        3241 :     if (T->maxpos && T->end >= T->sieveb) /* done with sieves ? */
     870             :     {
     871         175 :       if (T->sieveb == T->b && T->b != ULONG_MAX) return 0;
     872           1 :       T->strategy = PRST_unextprime;
     873             :     }
     874             :     else
     875             :     { /* initialize next chunk */
     876        3066 :       T->sieve = T->isieve;
     877        3066 :       if (T->maxpos == 0)
     878         245 :         T->a |= 1; /* first time; ensure odd */
     879             :       else
     880        2821 :         T->a = (T->end + 2) | 1;
     881        3066 :       T->end = T->a + T->chunk; /* may overflow */
     882        3066 :       if (T->end < T->a || T->end > T->sieveb) T->end = T->sieveb;
     883             :       /* end and a are odd; sieve[k] contains the a + 8*2k + (0,2,...,14).
     884             :        * The largest k is (end-a) >> 4 */
     885        3066 :       T->pos = 0;
     886        3066 :       T->maxpos = (T->end - T->a) >> 4;
     887        3066 :       sieve_block(T->a, T->end, T->maxpos, T->sieve);
     888        3066 :       goto NEXT_CHUNK;
     889             :     }
     890             :   }
     891     4141842 :   if (T->strategy == PRST_unextprime)
     892             :   {
     893     4141842 :     if (T->q == 1)
     894             :     {
     895             : #ifdef LONG_IS_64BIT
     896     4130946 :       if (T->p == (1UL<<63)) return T->p = 9223372036854775837UL;
     897      156654 :       if (T->p == 9223372036854775837UL) return T->p = 9223372036854775907UL;
     898             : #endif
     899      134093 :       T->p = unextprime(T->p + 1);
     900             :     }
     901             :     else
     902             :     {
     903             :       do {
     904        2863 :         T->p += T->q;
     905        2863 :         if (T->p < T->q) return 0; /* overflow */
     906        2863 :       } while (!uisprime(T->p));
     907             :     }
     908      134814 :     if (!T->p) /* overflow ulong, switch to GEN */
     909          15 :       T->strategy = PRST_nextprime;
     910             :     else
     911             :     {
     912      134799 :       if (T->p > T->b) return 0;
     913      128805 :       return T->p;
     914             :     }
     915             :   }
     916          15 :   return 0; /* overflow */
     917             : }
     918             : 
     919             : GEN
     920     8401452 : forprime_next(forprime_t *T)
     921             : {
     922             :   pari_sp av;
     923             :   GEN p;
     924     8401452 :   if (T->strategy != PRST_nextprime)
     925             :   {
     926     8393839 :     ulong q = u_forprime_next(T);
     927     8393839 :     if (q) { affui(q, T->pp); return T->pp; }
     928             :     /* failure */
     929         351 :     if (T->strategy != PRST_nextprime) return NULL; /* we're done */
     930             :     /* overflow ulong, switch to GEN */
     931          15 :     affui(ULONG_MAX, T->pp); /* initialize */
     932             :   }
     933        7628 :   av = avma;
     934        7628 :   p = nextprime(addiu(T->pp, 1));
     935        7628 :   if (T->bb && abscmpii(p, T->bb) > 0) { avma = av; return NULL; }
     936        7512 :   affii(p, T->pp); avma = av; return T->pp;
     937             : }
     938             : 
     939             : void
     940         301 : forprime(GEN a, GEN b, GEN code)
     941             : {
     942         301 :   pari_sp av = avma;
     943             :   forprime_t T;
     944             : 
     945         595 :   if (!forprime_init(&T, a,b)) { avma = av; return; }
     946             : 
     947         301 :   push_lex(T.pp,code);
     948        4452 :   while(forprime_next(&T))
     949             :   {
     950        3885 :     closure_evalvoid(code); if (loop_break()) break;
     951             :     /* p changed in 'code', complain */
     952        3857 :     if (get_lex(-1) != T.pp)
     953           7 :       pari_err(e_MISC,"prime index read-only: was changed to %Ps", get_lex(-1));
     954             :   }
     955         294 :   pop_lex(1); avma = av;
     956             : }
     957             : 
     958             : int
     959          42 : forcomposite_init(forcomposite_t *C, GEN a, GEN b)
     960             : {
     961          42 :   pari_sp av = avma;
     962          42 :   a = gceil(a); if (typ(a)!=t_INT) pari_err_TYPE("forcomposite",a);
     963          42 :   if (b) {
     964          35 :     b = gfloor(b);if (typ(b)!=t_INT) pari_err_TYPE("forcomposite",b);
     965             :   }
     966          42 :   if (signe(a) < 0) pari_err_DOMAIN("forcomposite", "a", "<", gen_0, a);
     967          42 :   if (abscmpiu(a, 4) < 0) a = utoipos(4);
     968          42 :   C->first = 1;
     969          42 :   if (!forprime_init(&C->T, a,b))
     970             :   {
     971           7 :     C->n = gen_1; /* in case caller forgets to check the return value */
     972           7 :     C->b = gen_0;
     973           7 :     avma = av; return 0;
     974             :   }
     975          35 :   C->n = setloop(a);
     976          35 :   C->b = b;
     977          35 :   C->p = NULL; return 1;
     978             : }
     979             : 
     980             : GEN
     981         147 : forcomposite_next(forcomposite_t *C)
     982             : {
     983         147 :   if (C->first) /* first call ever */
     984             :   {
     985          35 :     C->first = 0;
     986          35 :     C->p = forprime_next(&C->T);
     987             :   }
     988             :   else
     989         112 :     C->n = incloop(C->n);
     990         147 :   if (C->p)
     991             :   {
     992         119 :     if (cmpii(C->n, C->p) < 0) return C->n;
     993          56 :     C->n = incloop(C->n);
     994             :     /* n = p+1 */
     995          56 :     C->p = forprime_next(&C->T); /* nextprime(p) > n */
     996          56 :     if (C->p) return C->n;
     997             :   }
     998          49 :   if (!C->b || cmpii(C->n, C->b) <= 0) return C->n;
     999          21 :   return NULL;
    1000             : }
    1001             : 
    1002             : void
    1003          42 : forcomposite(GEN a, GEN b, GEN code)
    1004             : {
    1005          42 :   pari_sp av = avma;
    1006             :   forcomposite_t T;
    1007             :   GEN n;
    1008          77 :   if (!forcomposite_init(&T,a,b)) return;
    1009          35 :   push_lex(T.n,code);
    1010         182 :   while((n = forcomposite_next(&T)))
    1011             :   {
    1012         126 :     closure_evalvoid(code); if (loop_break()) break;
    1013             :     /* n changed in 'code', complain */
    1014         119 :     if (get_lex(-1) != n)
    1015           7 :       pari_err(e_MISC,"index read-only: was changed to %Ps", get_lex(-1));
    1016             :   }
    1017          28 :   pop_lex(1); avma = av;
    1018             : }

Generated by: LCOV version 1.11