Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - language - forprime.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.0 lcov report (development 29685-209b4888d8) Lines: 444 510 87.1 %
Date: 2024-11-05 09:09:42 Functions: 36 38 94.7 %
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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : /**********************************************************************/
      19             : /***                                                                ***/
      20             : /***                     Public prime table                         ***/
      21             : /***                                                                ***/
      22             : /**********************************************************************/
      23             : 
      24             : static ulong _maxprimelim = 0;
      25             : static GEN _prodprimes,_prodprimes_addr;
      26             : typedef unsigned char *byteptr;
      27             : 
      28             : /* Build/Rebuild table of prime differences. The actual work is done by the
      29             :  * following two subroutines;  the user entry point is the function
      30             :  * initprimes() below; initprimes1() is the basecase, called when
      31             :  * maxnum (size) is moderate. Must be called after pari_init_stack() )*/
      32             : static void
      33        1848 : initprimes1(ulong size, long *lenp, pari_prime *p1)
      34             : {
      35        1848 :   pari_sp av = avma;
      36             :   long k;
      37        1848 :   byteptr q, r, s, p = (byteptr)stack_calloc(size+2), fin = p + size;
      38             :   pari_prime *re;
      39             : 
      40       22176 :   for (r=q=p,k=1; r<=fin; )
      41             :   {
      42       33264 :     do { r+=k; k+=2; r+=k; } while (*++q);
      43      892584 :     for (s=r; s<=fin; s+=k) *s = 1;
      44             :   }
      45        1848 :   re = p1; *re++ = 2; *re++ = 3; /* 2 and 3 */
      46        1848 :   for (s=q=p+1; ; s=q)
      47             :   {
      48      944328 :     do q++; while (*q);
      49      316008 :     if (q > fin) break;
      50      314160 :     *re++ = (pari_prime) 2*(q-p)+1;
      51             :   }
      52        1848 :   *re++ = 0;
      53        1848 :   *lenp = re - p1;
      54        1848 :   set_avma(av);
      55        1848 : }
      56             : 
      57             : /*  Timing in ms (Athlon/850; reports 512K of secondary cache; looks
      58             :     like there is 64K of quickier cache too).
      59             : 
      60             :       arena|    30m     100m    300m    1000m    2000m  <-- primelimit
      61             :       =================================================
      62             :       16K       1.1053  1.1407  1.2589  1.4368   1.6086
      63             :       24K       1.0000  1.0625  1.1320  1.2443   1.3095
      64             :       32K       1.0000  1.0469  1.0761  1.1336   1.1776
      65             :       48K       1.0000  1.0000  1.0254  1.0445   1.0546
      66             :       50K       1.0000  1.0000  1.0152  1.0345   1.0464
      67             :       52K       1.0000  1.0000  1.0203  1.0273   1.0362
      68             :       54K       1.0000  1.0000  1.0812  1.0216   1.0281
      69             :       56K       1.0526  1.0000  1.0051  1.0144   1.0205
      70             :       58K       1.0000  1.0000  1.0000  1.0086   1.0123
      71             :       60K       0.9473  0.9844  1.0051  1.0014   1.0055
      72             :       62K       1.0000  0.9844  0.9949  0.9971   0.9993
      73             :       64K       1.0000  1.0000  1.0000  1.0000   1.0000
      74             :       66K       1.2632  1.2187  1.2183  1.2055   1.1953
      75             :       68K       1.4211  1.4844  1.4721  1.4425   1.4188
      76             :       70K       1.7368  1.7188  1.7107  1.6767   1.6421
      77             :       72K       1.9474  1.9531  1.9594  1.9023   1.8573
      78             :       74K       2.2105  2.1875  2.1827  2.1207   2.0650
      79             :       76K       2.4211  2.4219  2.4010  2.3305   2.2644
      80             :       78K       2.5789  2.6250  2.6091  2.5330   2.4571
      81             :       80K       2.8421  2.8125  2.8223  2.7213   2.6380
      82             :       84K       3.1053  3.1875  3.1776  3.0819   2.9802
      83             :       88K       3.5263  3.5312  3.5228  3.4124   3.2992
      84             :       92K       3.7895  3.8438  3.8375  3.7213   3.5971
      85             :       96K       4.0000  4.1093  4.1218  3.9986   3.9659
      86             :       112K      4.3684  4.5781  4.5787  4.4583   4.6115
      87             :       128K      4.7368  4.8750  4.9188  4.8075   4.8997
      88             :       192K      5.5263  5.7188  5.8020  5.6911   5.7064
      89             :       256K      6.0000  6.2187  6.3045  6.1954   6.1033
      90             :       384K      6.7368  6.9531  7.0405  6.9181   6.7912
      91             :       512K      7.3158  7.5156  7.6294  7.5000   7.4654
      92             :       768K      9.1579  9.4531  9.6395  9.5014   9.1075
      93             :       1024K    10.368  10.7497 10.9999 10.878   10.8201
      94             :       1536K    12.579  13.3124 13.7660 13.747   13.4739
      95             :       2048K    13.737  14.4839 15.0509 15.151   15.1282
      96             :       3076K    14.789  15.5780 16.2993 16.513   16.3365
      97             : 
      98             :     Now the same number relative to the model
      99             : 
     100             :     (1 + 0.36*sqrt(primelimit)/arena) * (arena <= 64 ? 1.05 : (arena-64)**0.38)
     101             : 
     102             :      [SLOW2_IN_ROOTS = 0.36, ALPHA = 0.38]
     103             : 
     104             :       arena|    30m     100m    300m    1000m    2000m  <-- primelimit
     105             :       =================================================
     106             :         16K    1.014    0.9835  0.9942  0.9889  1.004
     107             :         24K    0.9526   0.9758  0.9861  0.9942  0.981
     108             :         32K    0.971    0.9939  0.9884  0.9849  0.9806
     109             :         48K    0.9902   0.9825  0.996   0.9945  0.9885
     110             :         50K    0.9917   0.9853  0.9906  0.9926  0.9907
     111             :         52K    0.9932   0.9878  0.9999  0.9928  0.9903
     112             :         54K    0.9945   0.9902  1.064   0.9939  0.9913
     113             :         56K    1.048    0.9924  0.9925  0.993   0.9921
     114             :         58K    0.9969   0.9945  0.9909  0.9932  0.9918
     115             :         60K    0.9455   0.9809  0.9992  0.9915  0.9923
     116             :         62K    0.9991   0.9827  0.9921  0.9924  0.9929
     117             :         64K    1        1       1       1       1
     118             :         66K    1.02     0.9849  0.9857  0.9772  0.9704
     119             :         68K    0.8827   0.9232  0.9176  0.9025  0.8903
     120             :         70K    0.9255   0.9177  0.9162  0.9029  0.8881
     121             :         72K    0.9309   0.936   0.9429  0.9219  0.9052
     122             :         74K    0.9715   0.9644  0.967   0.9477  0.9292
     123             :         76K    0.9935   0.9975  0.9946  0.9751  0.9552
     124             :         78K    0.9987   1.021   1.021   1.003   0.9819
     125             :         80K    1.047    1.041   1.052   1.027   1.006
     126             :         84K    1.052    1.086   1.092   1.075   1.053
     127             :         88K    1.116    1.125   1.133   1.117   1.096
     128             :         92K    1.132    1.156   1.167   1.155   1.134
     129             :         96K    1.137    1.177   1.195   1.185   1.196
     130             :        112K    1.067    1.13    1.148   1.15    1.217
     131             :        128K    1.04     1.083   1.113   1.124   1.178
     132             :        192K    0.9368   0.985   1.025   1.051   1.095
     133             :        256K    0.8741   0.9224  0.9619  0.995   1.024
     134             :        384K    0.8103   0.8533  0.8917  0.9282  0.9568
     135             :        512K    0.7753   0.8135  0.8537  0.892   0.935
     136             :        768K    0.8184   0.8638  0.9121  0.9586  0.9705
     137             :       1024K    0.8241   0.8741  0.927   0.979   1.03
     138             :       1536K    0.8505   0.9212  0.9882  1.056   1.096
     139             :       2048K    0.8294   0.8954  0.9655  1.041   1.102
     140             : 
     141             : */
     142             : 
     143             : #ifndef SLOW2_IN_ROOTS
     144             :   /* SLOW2_IN_ROOTS below 3: some slowdown starts to be noticable
     145             :    * when things fit into the cache on Sparc.
     146             :    * The choice of 2.6 gives a slowdown of 1-2% on UltraSparcII,
     147             :    * but makes calculations for "maximum" of 436273009
     148             :    * fit into 256K cache (still common for some architectures).
     149             :    *
     150             :    * One may change it when small caches become uncommon, but the gain
     151             :    * is not going to be very noticable... */
     152             : #  ifdef i386           /* gcc defines this? */
     153             : #    define SLOW2_IN_ROOTS      0.36
     154             : #  else
     155             : #    define SLOW2_IN_ROOTS      2.6
     156             : #  endif
     157             : #endif
     158             : #ifndef CACHE_ARENA
     159             : #  ifdef i386           /* gcc defines this? */
     160             :    /* Due to smaller SLOW2_IN_ROOTS, smaller arena is OK; fit L1 cache */
     161             : #    define CACHE_ARENA (63 * 1024UL) /* No slowdown even with 64K L1 cache */
     162             : #  else
     163             : #    define CACHE_ARENA (200 * 1024UL) /* No slowdown even with 256K L2 cache */
     164             : #  endif
     165             : #endif
     166             : 
     167             : #define CACHE_ALPHA     (0.38)          /* Cache performance model parameter */
     168             : #define CACHE_CUTOFF    (0.018)         /* Cache performance not smooth here */
     169             : 
     170             : static double slow2_in_roots = SLOW2_IN_ROOTS;
     171             : 
     172             : typedef struct {
     173             :     ulong arena;
     174             :     double power;
     175             :     double cutoff;
     176             :     ulong bigarena;
     177             : } cache_model_t;
     178             : 
     179             : static cache_model_t cache_model = { CACHE_ARENA, CACHE_ALPHA, CACHE_CUTOFF, 0 };
     180             : 
     181             : /* Assume that some calculation requires a chunk of memory to be
     182             :    accessed often in more or less random fashion (as in sieving).
     183             :    Assume that the calculation can be done in steps by subdividing the
     184             :    chunk into smaller subchunks (arenas) and treating them
     185             :    separately.  Assume that the overhead of subdivision is equivalent
     186             :    to the number of arenas.
     187             : 
     188             :    Find an optimal size of the arena taking into account the overhead
     189             :    of subdivision, and the overhead of arena not fitting into the
     190             :    cache.  Assume that arenas of size slow2_in_roots slows down the
     191             :    calculation 2x (comparing to very big arenas; when cache hits do
     192             :    not matter).  Since cache performance varies wildly with
     193             :    architecture, load, and wheather (especially with cache coloring
     194             :    enabled), use an idealized cache model based on benchmarks above.
     195             : 
     196             :    Assume that an independent region of FIXED_TO_CACHE bytes is accessed
     197             :    very often concurrently with the arena access.
     198             :  */
     199             : static ulong
     200        1848 : good_arena_size(ulong slow2_size, ulong total, ulong fixed_to_cache,
     201             :                 cache_model_t *cache_model)
     202             : {
     203        1848 :   ulong asize, cache_arena = cache_model->arena;
     204             :   double Xmin, Xmax, A, B, C1, C2, D, V;
     205        1848 :   double alpha = cache_model->power, cut_off = cache_model->cutoff;
     206             : 
     207             :   /* Estimated relative slowdown,
     208             :      with overhead = max((fixed_to_cache+arena)/cache_arena - 1, 0):
     209             : 
     210             :      1 + slow2_size/arena due to initialization overhead;
     211             : 
     212             :      max(1, 4.63 * overhead^0.38 ) due to footprint > cache size.
     213             : 
     214             :      [The latter is hard to substantiate theoretically, but this
     215             :      function describes benchmarks pretty close; it does not hurt that
     216             :      one can minimize it explicitly too ;-).  The switch between
     217             :      different choices of max() happens when overhead=0.018.]
     218             : 
     219             :      Thus the problem is minimizing (1 + slow2_size/arena)*overhead**0.29.
     220             :      This boils down to F=((X+A)/(X+B))X^alpha, X=overhead,
     221             :      B = (1 - fixed_to_cache/cache_arena), A = B + slow2_size/cache_arena,
     222             :      alpha = 0.38, and X>=0.018, X>-B.
     223             : 
     224             :      We need to find the rightmost root of (X+A)*(X+B) - alpha(A-B)X to the
     225             :      right of 0.018 (if such exists and is below Xmax).  Then we manually
     226             :      check the remaining region [0, 0.018].
     227             : 
     228             :      Since we cannot trust the purely-experimental cache-hit slowdown
     229             :      function, as a sanity check always prefer fitting into the
     230             :      cache (or "almost fitting") if F-law predicts that the larger
     231             :      value of the arena provides less than 10% speedup.
     232             :    */
     233             : 
     234             :   /* The simplest case: we fit into cache */
     235        1848 :   asize = cache_arena - fixed_to_cache;
     236        1848 :   if (total <= asize) return total;
     237             :   /* The simple case: fitting into cache doesn't slow us down more than 10% */
     238        1848 :   if (asize > 10 * slow2_size) return asize;
     239             :   /* Slowdown of not fitting into cache is significant.  Try to optimize.
     240             :      Do not be afraid to spend some time on optimization - in trivial
     241             :      cases we do not reach this point; any gain we get should
     242             :      compensate the time spent on optimization.  */
     243             : 
     244           0 :   B = (1 - ((double)fixed_to_cache)/cache_arena);
     245           0 :   A = B + ((double)slow2_size)/cache_arena;
     246           0 :   C2 = A*B;
     247           0 :   C1 = (A + B - 1/alpha*(A - B))/2;
     248           0 :   D = C1*C1 - C2;
     249           0 :   if (D > 0)
     250           0 :     V = cut_off*cut_off + 2*C1*cut_off + C2; /* Value at CUT_OFF */
     251             :   else
     252           0 :     V = 0; /* Peacify the warning */
     253           0 :   Xmin = cut_off;
     254           0 :   Xmax = ((double)total - fixed_to_cache)/cache_arena; /* Two candidates */
     255             : 
     256           0 :   if ( D <= 0 || (V >= 0 && C1 + cut_off >= 0) ) /* slowdown increasing */
     257           0 :     Xmax = cut_off; /* Only one candidate */
     258           0 :   else if (V >= 0 && /* slowdown concave down */
     259           0 :            ((Xmax + C1) <= 0 || (Xmax*Xmax + 2*C1*Xmax + C2) <= 0))
     260             :       /* DO NOTHING */;  /* Keep both candidates */
     261           0 :   else if (V <= 0 && (Xmax*Xmax + 2*C1*Xmax + C2) <= 0) /*slowdown decreasing*/
     262           0 :       Xmin = cut_off; /* Only one candidate */
     263             :   else /* Now we know: 2 roots, the largest is in CUT_OFF..Xmax */
     264           0 :       Xmax = sqrt(D) - C1;
     265           0 :   if (Xmax != Xmin) { /* Xmin == CUT_OFF; Check which one is better */
     266           0 :     double v1 = (cut_off + A)/(cut_off + B);
     267           0 :     double v2 = 2.33 * (Xmax + A)/(Xmax + B) * pow(Xmax, alpha);
     268             : 
     269           0 :     if (1.1 * v2 >= v1) /* Prefer fitting into the cache if slowdown < 10% */
     270           0 :       V = v1;
     271             :     else
     272           0 :     { Xmin = Xmax; V = v2; }
     273           0 :   } else if (B > 0) /* We need V */
     274           0 :     V = 2.33 * (Xmin + A)/(Xmin + B) * pow(Xmin, alpha);
     275           0 :   if (B > 0 && 1.1 * V > A/B)  /* Now Xmin is the minumum.  Compare with 0 */
     276           0 :     Xmin = 0;
     277             : 
     278           0 :   asize = (ulong)((1 + Xmin)*cache_arena - fixed_to_cache);
     279           0 :   if (asize > total) asize = total; /* May happen due to approximations */
     280           0 :   return asize;
     281             : }
     282             : 
     283             : /* Use as in
     284             :     install(set_optimize,lLDG)          \\ Through some M too?
     285             :     set_optimize(2,1) \\ disable dependence on limit
     286             :     \\ 1: how much cache usable, 2: slowdown of setup, 3: alpha, 4: cutoff,
     287             :     \\ 5: cache size (typically whole L2 or L3) in bytes to use in forprime()
     288             :     \\ 2,3,4 are in units of 0.001
     289             : 
     290             :     { time_primes_arena(ar,limit) =     \\ ar = arena size in K
     291             :         set_optimize(1,floor(ar*1024));
     292             :         default(primelimit, 200 000);   \\ 100000 results in *larger* malloc()!
     293             :         gettime;
     294             :         default(primelimit, floor(limit));
     295             :         if(ar >= 1, ar=floor(ar));
     296             :         print("arena "ar"K => "gettime"ms");
     297             :     }
     298             : */
     299             : long
     300           0 : set_optimize(long what, GEN g)
     301             : {
     302           0 :   long ret = 0;
     303             : 
     304           0 :   switch (what) {
     305           0 :   case 1:
     306           0 :     ret = (long)cache_model.arena;
     307           0 :     break;
     308           0 :   case 2:
     309           0 :     ret = (long)(slow2_in_roots * 1000);
     310           0 :     break;
     311           0 :   case 3:
     312           0 :     ret = (long)(cache_model.power * 1000);
     313           0 :     break;
     314           0 :   case 4:
     315           0 :     ret = (long)(cache_model.cutoff * 1000);
     316           0 :     break;
     317           0 :   case 5:
     318           0 :     ret = (long)(cache_model.bigarena);
     319           0 :     break;
     320           0 :   default:
     321           0 :     pari_err_BUG("set_optimize");
     322           0 :     break;
     323             :   }
     324           0 :   if (g != NULL) {
     325           0 :     ulong val = itou(g);
     326             : 
     327           0 :     switch (what) {
     328           0 :     case 1: cache_model.arena = val; break;
     329           0 :     case 2: slow2_in_roots     = (double)val / 1000.; break;
     330           0 :     case 3: cache_model.power  = (double)val / 1000.; break;
     331           0 :     case 4: cache_model.cutoff = (double)val / 1000.; break;
     332           0 :     case 5: cache_model.bigarena = val; break;
     333             :     }
     334             :   }
     335           0 :   return ret;
     336             : }
     337             : 
     338             : /* s is odd; prime (starting from 3 = known_primes[2]), terminated by a 0 byte.
     339             :  * Checks n odd numbers starting at 'start', setting bytes to 0 (composite)
     340             :  * or 1 (prime), starting at data */
     341             : static void
     342        7080 : sieve_chunk(pari_prime *known_primes, ulong s, byteptr data, ulong n)
     343             : {
     344        7080 :   ulong p, cnt = n-1, start = s;
     345             :   pari_prime *q;
     346             : 
     347        7080 :   memset(data, 0, n);
     348        7080 :   start >>= 1;  /* (start - 1)/2 */
     349        7080 :   start += n; /* Corresponds to the end */
     350             :   /* data corresponds to start, q runs over primediffs */
     351     1017120 :   for (q = known_primes + 1, p = 3; p; p = *++q)
     352             :   { /* first odd number >= start > p and divisible by p
     353             :        = last odd number <= start + 2p - 2 and 0 (mod p)
     354             :        = p + last number <= start + p - 2 and 0 (mod 2p)
     355             :        = p + start+p-2 - (start+p-2) % 2p
     356             :        = start + 2(p - 1 - ((start-1)/2 + (p-1)/2) % p). */
     357     1010040 :     long off = cnt - ((start+(p>>1)) % p);
     358  1613194432 :     while (off >= 0) { data[off] = 1; off -= p; }
     359             :   }
     360        7080 : }
     361             : 
     362             : static void
     363        1848 : set_prodprimes(void)
     364             : {
     365        1848 :   pari_sp ltop = avma, av;
     366        1848 :   ulong b = 1UL << 8, M = minuu(maxprime(), GP_DATA->factorlimit);
     367        1848 :   GEN W, w, v = primes_interval_zv(3, M);
     368        1848 :   long s, u, j, jold, l = lg(v);
     369             : 
     370        1848 :   W = cgetg(64+1, t_VEC);
     371   151582200 :   for (jold = j = u = 1; j < l; j++)
     372   151580352 :     if (j==l-1 || uel(v,j) >= b)
     373             :     {
     374       24024 :       long lw = (j == l-1? l: j) - jold + 1;
     375       24024 :       w = v+jold-1; w[0] = evaltyp(t_VECSMALL) | _evallg(lw);
     376       24024 :       gel(W,u++) = zv_prod_Z(w); /* p_jold ... p_{j-1} */
     377       24024 :       jold = j; b *= 2;
     378       24024 :       if (b > M) b = M; /* truncate last run */
     379             :     }
     380        1848 :   setlg(W, u);
     381       24024 :   for (j = 2; j < u; j++) gel(W,j) = mulii(gel(W,j-1), gel(W,j));
     382        1848 :   s = gsizeword(W);
     383        1848 :   w = (GEN)pari_malloc(s*sizeof(long));
     384        1848 :   av = (pari_sp)(w + s);
     385        1848 :   _prodprimes_addr = w;
     386        1848 :   _prodprimes = gcopy_avma(W, &av);
     387        1848 :   set_avma(ltop);
     388        1848 : }
     389             : 
     390             : static void
     391        1848 : initprimes0(ulong maxnum, long *lenp, pari_prime *p1)
     392             : {
     393        1848 :   pari_sp av = avma, bot = pari_mainstack->bot;
     394             :   long alloced, psize;
     395             :   byteptr q, end, p;
     396             :   ulong remains, curlow, rootnum, asize, prime_above, last;
     397             :   pari_prime *end1, *curdiff, *p_prime_above;
     398             : 
     399        1848 :   if (!odd(maxnum)) maxnum--; /* make it odd. */
     400             :   /* base case */
     401        1848 :   if (maxnum < 1ul<<17) { initprimes1(maxnum>>1, lenp, p1); return; }
     402             : 
     403             :   /* Checked to be enough up to 40e6, attained at 155893 */
     404        1848 :   rootnum = usqrt(maxnum) | 1;
     405        1848 :   initprimes1(rootnum>>1, &psize, p1);
     406        1848 :   last = rootnum;
     407        1848 :   end1 = p1 + psize - 1;
     408        1848 :   remains = (maxnum - last) >> 1; /* number of odd numbers to check */
     409             :   /* we access primes array of psize too; but we access it consecutively,
     410             :    * thus we do not include it in fixed_to_cache */
     411        1848 :   asize = good_arena_size((ulong)(rootnum * slow2_in_roots), remains+1, 0,
     412             :                           &cache_model) - 1;
     413             :   /* enough room on the stack ? */
     414        1848 :   alloced = (((byteptr)avma) <= ((byteptr)bot) + asize);
     415        1848 :   p = (byteptr)(alloced? pari_malloc(asize+1): stack_malloc(asize+1));
     416        1848 :   end = p + asize; /* the 0 sentinel goes at end. */
     417        1848 :   curlow = last + 2; /* First candidate: know primes up to last (odd). */
     418        1848 :   curdiff = end1;
     419             : 
     420             :   /* During each iteration p..end-1 represents a range of odd
     421             :      numbers.   */
     422        1848 :   p_prime_above = p1 + 2;
     423        1848 :   prime_above = 3;
     424        8928 :   while (remains)
     425             :   { /* cycle over arenas; performance not crucial */
     426             :     pari_prime was_delta;
     427        7080 :     if (asize > remains) { asize = remains; end = p + asize; }
     428             :     /* Fake the upper limit appropriate for the given arena */
     429      323088 :     while (prime_above*prime_above <= curlow + (asize << 1) && *p_prime_above)
     430      316008 :       prime_above = *p_prime_above++;
     431        7080 :     was_delta = *p_prime_above;
     432        7080 :     *p_prime_above = 0; /* sentinel for sieve_chunk */
     433        7080 :     sieve_chunk(p1, curlow, p, asize);
     434        7080 :     *p_prime_above = was_delta; /* restore */
     435             : 
     436        7080 :     p[asize] = 0; /* sentinel */
     437        7080 :     for (q = p; ; q++)
     438             :     { /* q runs over addresses corresponding to primes */
     439   967945128 :       while (*q) q++; /* use sentinel at end */
     440   151271424 :       if (q >= end) break;
     441   151264344 :       *curdiff++ = (pari_prime) 2*(q-p) + curlow;
     442             :     }
     443        7080 :     remains -= asize;
     444        7080 :     curlow += (asize<<1);
     445             :   }
     446        1848 :   *curdiff++ = 0; /* sentinel */
     447        1848 :   *lenp = curdiff - p1;
     448        1848 :   if (alloced) pari_free(p); else set_avma(av);
     449             : }
     450             : 
     451             : ulong
     452    46711173 : maxprime(void) { return pari_PRIMES? pari_PRIMES[pari_PRIMES[0]]: 0; }
     453             : ulong
     454    70312701 : maxprimelim(void) { return pari_PRIMES? _maxprimelim: 0; }
     455             : ulong
     456         196 : maxprimeN(void) { return pari_PRIMES? pari_PRIMES[0]: 0; }
     457             : GEN
     458     2691311 : prodprimes(void) { return pari_PRIMES? _prodprimes: NULL; }
     459             : void
     460           0 : maxprime_check(ulong c) { if (maxprime() < c) pari_err_MAXPRIME(c); }
     461             : 
     462             : static pari_prime*
     463        1848 : initprimes(ulong maxnum)
     464             : {
     465             :   pari_prime *t;
     466             :   long len;
     467             :   ulong N;
     468        1848 :   if (maxnum < 65537)
     469             :   {
     470           0 :     maxnum = 65537;
     471           0 :     N = 6543;
     472             :   }
     473             :   else
     474        1848 :     N = (long) ceil(primepi_upper_bound((double)maxnum));
     475        1848 :   t = (pari_prime*) pari_malloc(sizeof(*t) * (N+2));
     476        1848 :   initprimes0(maxnum, &len, t+1); t[0] = (pari_prime)(len-1);
     477        1848 :   _maxprimelim = maxnum;
     478        1848 :   return (pari_prime*) pari_realloc(t, sizeof(*t) * (len+1));
     479             : }
     480             : 
     481             : void
     482        1848 : initprimetable(ulong maxnum)
     483             : {
     484        1848 :   pari_prime *old = pari_PRIMES;
     485             : #ifdef LONG_IS_64BIT
     486        1592 :   maxnum = minuu(maxnum, 4294967295);
     487             : #endif
     488        1848 :   pari_PRIMES = initprimes(maxnum);
     489        1848 :   if (old) free(old);
     490        1848 :   set_prodprimes();
     491        1848 : }
     492             : 
     493             : /**********************************************************************/
     494             : /***                                                                ***/
     495             : /***                     forprime                                   ***/
     496             : /***                                                                ***/
     497             : /**********************************************************************/
     498             : 
     499             : /* return good chunk size for sieve, 16 | chunk + 2 */
     500             : static ulong
     501     8476266 : optimize_chunk(ulong a, ulong b)
     502             : {
     503             :   /* TODO: Optimize size (surely < 512k to stay in L2 cache, but not so large
     504             :    * as to force recalculating too often). */
     505             :   /* bigarena is in bytes, we use bits, and only odds */
     506     8476266 :   ulong defchunk = (a>>31) > 1 ? 0x80000UL : 0x8000;
     507     8476266 :   ulong chunk = (cache_model.bigarena ? cache_model.bigarena : defchunk)<<4;
     508     8476266 :   ulong tmp = (b - a) / chunk + 1;
     509             : 
     510     8476266 :   if (tmp == 1)
     511         196 :     chunk = b - a + 16;
     512             :   else
     513     8476070 :     chunk = (b - a) / tmp + 15;
     514             :   /* ensure 16 | chunk + 2 */
     515     8476266 :   return (((chunk + 2)>>4)<<4) - 2;
     516             : }
     517             : static void
     518     8476267 : sieve_init(forprime_t *T, ulong a, ulong b)
     519             : {
     520     8476267 :   T->sieveb = b;
     521     8476267 :   T->chunk = optimize_chunk(a, b);
     522             :   /* >> 1 [only odds] + 3 [convert from bits to bytes] */
     523     8476263 :   T->isieve = (unsigned char*)stack_malloc(((T->chunk+2) >> 4) + 1);
     524     8476233 :   T->cache[0] = 0;
     525     8476233 :   T->a = a;
     526     8476233 :   T->end = minuu(a + T->chunk, b);
     527     8476225 :   T->pos = T->maxpos = 0;
     528     8476225 : }
     529             : 
     530             : enum {PRST_none, PRST_diffptr, PRST_sieve, PRST_unextprime, PRST_nextprime};
     531             : 
     532             : static void
     533    13682066 : u_forprime_set_prime_table(forprime_t *T, ulong a)
     534             : {
     535    13682066 :   T->strategy = PRST_diffptr;
     536    13682066 :   if (a < 3)
     537             :   {
     538     2393967 :     T->p = 0;
     539     2393967 :     T->n = 0;
     540             :   }
     541             :   else
     542             :   {
     543    11288099 :     long n = PRIMES_search(a - 1);
     544    11290864 :     if (n < 0) n = - n - 1;
     545    11290864 :     T->n = n;
     546    11290864 :     T->p = pari_PRIMES[n];
     547             :   }
     548    13684831 : }
     549             : 
     550             : /* Set p so that p + q the smallest integer = c (mod q) and > original p.
     551             :  * Assume 0 < c < q. */
     552             : static void
     553      101982 : arith_set(forprime_t *T)
     554             : {
     555      101982 :   ulong r = T->p % T->q; /* 0 <= r <= min(p, q-1) */
     556      101982 :   pari_sp av = avma;
     557      101982 :   GEN d = adduu(T->p - r, T->c); /* = c mod q */
     558      101982 :   if (T->c > r) d = subiu(d, T->q);
     559             :   /* d = c mod q,  d = c > r? p-r+c-q: p-r+c, so that
     560             :    *  d <= p  and  d+q = c>r? p-r+c  : p-r+c+q > p */
     561      101982 :   if (signe(d) <= 0)
     562             :   {
     563          20 :     T->p = 0;
     564          20 :     T->strategy = PRST_nextprime;
     565          20 :     affii(d, T->pp);
     566             :   }
     567             :   else
     568      101962 :     T->p = itou_or_0(d);
     569      101982 :   set_avma(av);
     570      101982 : }
     571             : 
     572             : /* Run through primes in arithmetic progression = c (mod q).
     573             :  * Warning: b = ULONG_MAX may signal that we are called by higher level
     574             :  * function handling a continuation for larger b; this sentinel value
     575             :  * must not be modified */
     576             : static int
     577    29000365 : u_forprime_sieve_arith_init(forprime_t *T, struct pari_sieve *psieve,
     578             :                             ulong a, ulong b, ulong c, ulong q)
     579             : {
     580             : #ifdef LONG_IS_64BIT
     581    24867274 :   const ulong UPRIME_MAX = 18446744073709551557UL;
     582             : #else
     583     4133091 :   const ulong UPRIME_MAX = 4294967291UL;
     584             : #endif
     585             :   ulong Plim, P, P2, Y, sieveb;
     586             : 
     587    29000365 :   if (!odd(b) && b > 2) b--;
     588    28999970 :   if (a > b || b < 2)
     589             :   {
     590      883385 :     T->strategy = PRST_diffptr; /* paranoia */
     591      883385 :     T->p = 0; /* empty */
     592      883385 :     T->b = 0; /* empty */
     593      883385 :     T->n = 0;
     594      883385 :     return 0;
     595             :   }
     596    28116585 :   P = maxprime();
     597    28115886 :   if (b != ULONG_MAX && b > UPRIME_MAX) b = UPRIME_MAX;
     598    28115886 :   if (q != 1)
     599             :   {
     600             :     ulong D;
     601      587768 :     c %= q; D = ugcd(c, q);
     602      587765 :     if (D != 1) { a = maxuu(a,D); if (b != ULONG_MAX) b = minuu(b,D); }
     603      587765 :     if (odd(q) && (a > 2 || c != 2))
     604             :     { /* only *odd* primes. If a <= c = 2, then p = 2 must be included :-( */
     605      510372 :       if (!odd(c)) c += q;
     606      510468 :       q <<= 1;
     607             :     }
     608             :   }
     609    28115966 :   T->q = q;
     610    28115966 :   T->c = c;
     611    28115966 :   T->strategy = PRST_none; /* unknown */
     612    28115966 :   T->psieve = psieve; /* unused for now */
     613    28115966 :   T->isieve = NULL; /* unused for now */
     614    28115966 :   T->b = b;
     615    28115966 :   if (P >= b) { /* [a,b] \subset prime table */
     616    10079299 :     u_forprime_set_prime_table(T, a);
     617    10081273 :     return 1;
     618             :   }
     619             :   /* b > P */
     620    18036667 :   if (a >= P)
     621             :   {
     622    14433814 :     T->p = a - 1;
     623    14433814 :     if (T->q != 1) arith_set(T);
     624             :   }
     625             :   else
     626     3602853 :     u_forprime_set_prime_table(T, a);
     627    18036654 :   if (T->strategy == PRST_none) T->strategy = PRST_unextprime;
     628             :   /* now strategy is either PRST_diffptr or PRST_unextprime */
     629             : 
     630    18036654 :   P2 = (P & HIGHMASK)? 0 : P*P;
     631    18036654 :   sieveb = b; if (P2 && P2 < b) sieveb = P2;
     632             :   /* maxprime^2 >= sieveb */
     633    18036654 :   Plim = maxprimelim();
     634    18036632 :   if (a <= Plim) a = Plim + 1;
     635    18036632 :   if (sieveb < a + 16) return 1;
     636     8978244 :   Y = sieveb - a + 1; /* number of integers in sievable interval > 16 */
     637     8978244 :   P = usqrt(sieveb); /* largest sieving prime */
     638             :   /* FIXME: should sieve as well if q != 1, adapt sieve code */
     639     8978407 :   if (q == 1 && (!P2 || P2 > a) && 3/M_LN2 * Y >= uprimepi(P))
     640             :   /* Sieve implemented & possible & not too costly. Cost model is
     641             :    * - nextprime: about Y / log(b) primes to test [neglect cost for composites]
     642             :    *   individual cost average = 3 log2(b) mulmod, total = 3 Y / log(2) mulmod
     643             :    * - sieve: pi(P) mod + Y loglog(b) add
     644             :    * Since loglog(b) < 4, and add < 10*mulmod, we neglect the Y loglog(b) term.
     645             :    * We have mod < mulmod < 2*mod; for now, assume mulmod ~ mod. */
     646             :   {
     647     8476271 :     if (T->strategy == PRST_unextprime) T->strategy = PRST_sieve;
     648     8476271 :     sieve_init(T, a, sieveb);
     649             :   }
     650     8978484 :   return 1;
     651             : }
     652             : 
     653             : int
     654    23635735 : u_forprime_arith_init(forprime_t *T, ulong a, ulong b, ulong c, ulong q)
     655    23635735 : { return u_forprime_sieve_arith_init(T, NULL, a, b, c, q); }
     656             : 
     657             : /* will run through primes in [a,b] */
     658             : int
     659    23043437 : u_forprime_init(forprime_t *T, ulong a, ulong b)
     660    23043437 : { return u_forprime_arith_init(T, a,b, 0,1); }
     661             : 
     662             : /* will run through primes in [a,b] */
     663             : static int
     664     5356836 : u_forprime_sieve_init(forprime_t *T, struct pari_sieve *s, ulong b)
     665     5356836 : { return u_forprime_sieve_arith_init(T, s, s->start, b, s->c, s->q); }
     666             : 
     667             : /* now only run through primes <= c; assume c <= b above */
     668             : void
     669          63 : u_forprime_restrict(forprime_t *T, ulong c) { T->b = c; }
     670             : 
     671             : /* b = NULL: loop forever */
     672             : int
     673        2420 : forprimestep_init(forprime_t *T, GEN a, GEN b, GEN q)
     674             : {
     675        2420 :   GEN c = NULL;
     676             :   long lb;
     677             : 
     678        2420 :   a = gceil(a); if (typ(a) != t_INT) pari_err_TYPE("forprime_init",a);
     679        2420 :   T->qq = NULL; T->q = 1; T->c = 0;
     680        2420 :   if (q)
     681             :   {
     682         133 :     switch(typ(q))
     683             :     {
     684          56 :       case t_INT:
     685          56 :         c = a; break;
     686          77 :       case t_INTMOD:
     687          77 :         c = gel(q,2); q = gel(q,1);
     688             :         /* first int >= initial a which is = c (mod q) */
     689          77 :         a = addii(a, modii(subii(c,a), q)); break;
     690           0 :       default: pari_err_TYPE("forprimestep_init",q);
     691             :     }
     692         133 :     if (signe(q) <= 0) pari_err_TYPE("forprimestep_init (q <= 0)",q);
     693         133 :     if (equali1(q)) c = q = NULL;
     694             :     else
     695             :     {
     696         133 :       GEN D = gcdii(c, q);
     697         133 :       if (!is_pm1(D))
     698             :       { /* at most one prime: c */
     699          42 :         if (cmpii(a, D) < 0) a = D;
     700          42 :         if (gcmp(b, D) > 0) b = D;
     701             :       }
     702         133 :       if ((T->q = itou_or_0(q)))
     703         125 :         T->c = umodiu(c, T->q);
     704             :       else
     705           8 :         T->qq = q;
     706             :     }
     707             :   }
     708        2420 :   if (signe(a) <= 0) a = q? modii(a, q): gen_1;
     709        2420 :   if (b && typ(b) != t_INFINITY)
     710             :   {
     711        1013 :     b = gfloor(b);
     712        1013 :     if (typ(b) != t_INT) pari_err_TYPE("forprime_init",b);
     713        1013 :     if (signe(b) < 0 || cmpii(a,b) > 0)
     714             :     {
     715          21 :       T->strategy = PRST_nextprime; /* paranoia */
     716          21 :       T->bb = T->pp = gen_0; return 0;
     717             :     }
     718         992 :     lb = lgefint(b);
     719         992 :     T->bb = b;
     720             :   }
     721        1407 :   else if (!b || inf_get_sign(b) > 0)
     722             :   {
     723        1407 :     lb = lgefint(a) + 4;
     724        1407 :     T->bb = NULL;
     725             :   }
     726             :   else /* b == -oo */
     727             :   {
     728           0 :     T->strategy = PRST_nextprime; /* paranoia */
     729           0 :     T->bb = T->pp = gen_0; return 0;
     730             :   }
     731        2399 :   T->pp = cgeti(T->qq? maxuu(lb, lgefint(T->qq)): lb);
     732             :   /* a, b are positive integers, a <= b */
     733        2399 :   if (!T->qq && lgefint(a) == 3) /* lb == 3 implies b != NULL */
     734        2256 :     return u_forprime_arith_init(T, uel(a,2), lb == 3? uel(b,2): ULONG_MAX,
     735             :                                     T->c, T->q);
     736         143 :   T->strategy = PRST_nextprime;
     737         143 :   affii(T->qq? subii(a,T->qq): subiu(a,T->q), T->pp); return 1;
     738             : }
     739             : int
     740        1315 : forprime_init(forprime_t *T, GEN a, GEN b)
     741        1315 : { return forprimestep_init(T,a,b,NULL); }
     742             : 
     743             : /* assume a <= b <= maxprime()^2, a,b odd, sieve[n] corresponds to
     744             :  *   a+16*n, a+16*n+2, ..., a+16*n+14 (bits 0 to 7)
     745             :  * maxpos = index of last sieve cell.
     746             :  * b-a+2 must be divisible by 16 for use by u_forprime_next */
     747             : static void
     748        9115 : sieve_block(ulong a, ulong b, ulong maxpos, unsigned char* sieve)
     749             : {
     750        9115 :   ulong i, lim = usqrt(b), sz = (b-a) >> 1;
     751        9115 :   (void)memset(sieve, 0, maxpos+1);
     752        9115 :   for (i = 2;; i++)
     753    24415778 :   { /* p is odd */
     754    24424893 :     ulong k, r, p = pari_PRIMES[i]; /* starts at p = 3 */
     755    24424893 :     if (p > lim) break;
     756             : 
     757             :     /* solve a + 2k = 0 (mod p) */
     758    24415778 :     r = a % p;
     759    24415778 :     if (r == 0)
     760       16028 :       k = 0;
     761             :     else
     762             :     {
     763    24399750 :       k = p - r;
     764    24399750 :       if (odd(k)) k += p;
     765    24399750 :       k >>= 1;
     766             :     }
     767             :     /* m = a + 2k is the smallest odd m >= a, p | m */
     768             :     /* position n (corresponds to a+2n) is sieve[n>>3], bit n&7 */
     769  5716476299 :     while (k <= sz) { sieve[k>>3] |= 1 << (k&7); k += p; /* 2k += 2p */ }
     770             :   }
     771        9115 : }
     772             : 
     773             : static void
     774        1848 : pari_sieve_init(struct pari_sieve *s, ulong a, ulong b)
     775             : {
     776        1848 :   ulong maxpos= (b - a) >> 4;
     777        1848 :   s->start = a; s->end = b;
     778        1848 :   s->sieve = (unsigned char*) pari_malloc(maxpos+1);
     779        1848 :   s->c = 0; s->q = 1;
     780        1848 :   sieve_block(a, b, maxpos, s->sieve);
     781        1848 :   s->maxpos = maxpos; /* must be last in case of SIGINT */
     782        1848 : }
     783             : 
     784             : static struct pari_sieve pari_sieve_modular;
     785             : 
     786             : #ifdef LONG_IS_64BIT
     787             : #define PARI_MODULAR_BASE ((1UL<<((BITS_IN_LONG-2)>>1))+1)
     788             : #else
     789             : #define PARI_MODULAR_BASE ((1UL<<(BITS_IN_LONG-1))+1)
     790             : #endif
     791             : 
     792             : void
     793        1848 : pari_init_primes(ulong maxprime)
     794             : {
     795        1848 :   ulong a = PARI_MODULAR_BASE, b = a + (1UL<<20)-2;
     796        1848 :   initprimetable(maxprime);
     797        1848 :   pari_sieve_init(&pari_sieve_modular, a, b);
     798        1848 : }
     799             : 
     800             : void
     801        1848 : pari_close_primes(void)
     802             : {
     803        1848 :   if (pari_PRIMES)
     804             :   {
     805        1848 :     pari_free(pari_PRIMES);
     806        1848 :     pari_free(_prodprimes_addr);
     807             :   }
     808        1848 :   pari_free(pari_sieve_modular.sieve);
     809        1848 : }
     810             : 
     811             : void
     812     4503084 : init_modular_small(forprime_t *S)
     813             : {
     814             : #ifdef LONG_IS_64BIT
     815     3859836 :   u_forprime_sieve_init(S, &pari_sieve_modular, ULONG_MAX);
     816             : #else
     817      643248 :   ulong a = (1UL<<((BITS_IN_LONG-2)>>1))+1;
     818      643248 :   u_forprime_init(S, a, ULONG_MAX);
     819             : #endif
     820     4503107 : }
     821             : 
     822             : void
     823    10468724 : init_modular_big(forprime_t *S)
     824             : {
     825             : #ifdef LONG_IS_64BIT
     826     8971723 :   u_forprime_init(S, HIGHBIT + 1, ULONG_MAX);
     827             : #else
     828     1497001 :   u_forprime_sieve_init(S, &pari_sieve_modular, ULONG_MAX);
     829             : #endif
     830    10468820 : }
     831             : 
     832             : /* T->cache is a 0-terminated list of primes, return the first one and
     833             :  * remove it from list. Most of the time the list contains a single prime */
     834             : static ulong
     835   129994938 : shift_cache(forprime_t *T)
     836             : {
     837             :   long i;
     838   129994938 :   T->p = T->cache[0];
     839   173799704 :   for (i = 1;; i++)  /* remove one prime from cache */
     840   173799704 :     if (! (T->cache[i-1] = T->cache[i]) ) break;
     841   129994938 :   return T->p;
     842             : }
     843             : 
     844             : ulong
     845   207632412 : u_forprime_next(forprime_t *T)
     846             : {
     847   207632412 :   if (T->strategy == PRST_diffptr)
     848             :   {
     849             :     for(;;)
     850             :     {
     851   221015442 :       if (++T->n <= pari_PRIMES[0])
     852             :       {
     853   221015281 :         T->p = pari_PRIMES[T->n];
     854   221015281 :         if (T->p > T->b) return 0;
     855   220828318 :         if (T->q == 1 || T->p % T->q == T->c) return T->p;
     856             :       }
     857             :       else
     858             :       { /* beyond the table */
     859         161 :         T->strategy = T->isieve? PRST_sieve: PRST_unextprime;
     860         161 :         if (T->q != 1) { arith_set(T); if (!T->p) return 0; }
     861             :         /* T->p possibly not a prime if q != 1 */
     862         161 :         break;
     863             :       }
     864             :     }
     865             :   }
     866   144750369 :   if (T->strategy == PRST_sieve)
     867             :   { /* require sieveb - a >= 16 */
     868             :     ulong n;
     869   129995499 :     if (T->cache[0]) return shift_cache(T);
     870    92980384 : NEXT_CHUNK:
     871    92988519 :     if (T->psieve)
     872             :     {
     873     5356855 :       T->sieve = T->psieve->sieve;
     874     5356855 :       T->end = T->psieve->end;
     875     5356855 :       if (T->end > T->sieveb) T->end = T->sieveb;
     876     5356855 :       T->maxpos = T->psieve->maxpos;
     877     5356855 :       T->pos = 0;
     878     5356855 :       T->psieve = NULL;
     879             :     }
     880   140784477 :     for (n = T->pos; n < T->maxpos; n++)
     881   140774003 :       if (T->sieve[n] != 0xFF)
     882             :       {
     883    92978045 :         unsigned char mask = T->sieve[n];
     884    92978045 :         ulong p = T->a + (n<<4);
     885    92978045 :         long i = 0;
     886    92978045 :         T->pos = n;
     887    92978045 :         if (!(mask &  1)) T->cache[i++] = p;
     888    92978045 :         if (!(mask &  2)) T->cache[i++] = p+2;
     889    92978045 :         if (!(mask &  4)) T->cache[i++] = p+4;
     890    92978045 :         if (!(mask &  8)) T->cache[i++] = p+6;
     891    92978045 :         if (!(mask & 16)) T->cache[i++] = p+8;
     892    92978045 :         if (!(mask & 32)) T->cache[i++] = p+10;
     893    92978045 :         if (!(mask & 64)) T->cache[i++] = p+12;
     894    92978045 :         if (!(mask &128)) T->cache[i++] = p+14;
     895    92978045 :         T->cache[i] = 0;
     896    92978045 :         T->pos = n+1;
     897    92978045 :         return shift_cache(T);
     898             :       }
     899             :     /* n = T->maxpos, last cell: check p <= b */
     900       10474 :     if (T->maxpos && n == T->maxpos && T->sieve[n] != 0xFF)
     901             :     {
     902        2873 :       unsigned char mask = T->sieve[n];
     903        2873 :       ulong p = T->a + (n<<4);
     904        2873 :       long i = 0;
     905        2873 :       T->pos = n;
     906        2873 :       if (!(mask &  1) && p <= T->sieveb) T->cache[i++] = p;
     907        2873 :       if (!(mask &  2) && p <= T->sieveb-2) T->cache[i++] = p+2;
     908        2873 :       if (!(mask &  4) && p <= T->sieveb-4) T->cache[i++] = p+4;
     909        2873 :       if (!(mask &  8) && p <= T->sieveb-6) T->cache[i++] = p+6;
     910        2873 :       if (!(mask & 16) && p <= T->sieveb-8) T->cache[i++] = p+8;
     911        2873 :       if (!(mask & 32) && p <= T->sieveb-10) T->cache[i++] = p+10;
     912        2873 :       if (!(mask & 64) && p <= T->sieveb-12) T->cache[i++] = p+12;
     913        2873 :       if (!(mask &128) && p <= T->sieveb-14) T->cache[i++] = p+14;
     914        2873 :       if (i)
     915             :       {
     916        2772 :         T->cache[i] = 0;
     917        2772 :         T->pos = n+1;
     918        2772 :         return shift_cache(T);
     919             :       }
     920             :     }
     921             : 
     922        7702 :     if (T->maxpos && T->end >= T->sieveb) /* done with sieves ? */
     923             :     {
     924         419 :       if (T->sieveb == T->b && T->b != ULONG_MAX) return 0;
     925           1 :       T->strategy = PRST_unextprime;
     926             :     }
     927             :     else
     928             :     { /* initialize next chunk */
     929        7283 :       T->sieve = T->isieve;
     930        7283 :       if (T->maxpos == 0)
     931        3366 :         T->a |= 1; /* first time; ensure odd */
     932             :       else
     933        3917 :         T->a = (T->end + 2) | 1;
     934        7283 :       T->end = T->a + T->chunk; /* may overflow */
     935        7283 :       if (T->end < T->a || T->end > T->sieveb) T->end = T->sieveb;
     936             :       /* end and a are odd; sieve[k] contains the a + 8*2k + (0,2,...,14).
     937             :        * The largest k is (end-a) >> 4 */
     938        7283 :       T->pos = 0;
     939        7283 :       T->maxpos = (T->end - T->a) >> 4; /* >= 1 */
     940        7283 :       sieve_block(T->a, T->end, T->maxpos, T->sieve);
     941        8135 :       goto NEXT_CHUNK;
     942             :     }
     943             :   }
     944    14754871 :   if (T->strategy == PRST_unextprime)
     945             :   {
     946    14753827 :     if (T->q == 1)
     947             :     {
     948             : #ifdef LONG_IS_64BIT
     949    14600706 :       switch(T->p)
     950             :       {
     951             : #define retp(x) return T->p = (HIGHBIT+x <= T->b)? HIGHBIT+x: 0
     952     8971788 :         case HIGHBIT: retp(29);
     953     3180270 :         case HIGHBIT + 29: retp(99);
     954      354950 :         case HIGHBIT + 99: retp(123);
     955      204687 :         case HIGHBIT +123: retp(131);
     956      145524 :         case HIGHBIT +131: retp(155);
     957      124662 :         case HIGHBIT +155: retp(255);
     958      103368 :         case HIGHBIT +255: retp(269);
     959       93603 :         case HIGHBIT +269: retp(359);
     960       75202 :         case HIGHBIT +359: retp(435);
     961       55933 :         case HIGHBIT +435: retp(449);
     962       49121 :         case HIGHBIT +449: retp(453);
     963       46243 :         case HIGHBIT +453: retp(485);
     964       40256 :         case HIGHBIT +485: retp(491);
     965       37171 :         case HIGHBIT +491: retp(543);
     966       34772 :         case HIGHBIT +543: retp(585);
     967       32167 :         case HIGHBIT +585: retp(599);
     968       28170 :         case HIGHBIT +599: retp(753);
     969       27390 :         case HIGHBIT +753: retp(849);
     970       26334 :         case HIGHBIT +849: retp(879);
     971       24716 :         case HIGHBIT +879: retp(885);
     972       24015 :         case HIGHBIT +885: retp(903);
     973       23535 :         case HIGHBIT +903: retp(995);
     974             : #undef retp
     975             :       }
     976             : #endif
     977      896901 :       T->p = unextprime(T->p + 1);
     978      896940 :       if (T->p > T->b) return 0;
     979             :     }
     980             :     else do {
     981     2798471 :       T->p += T->q;
     982     2798471 :       if (T->p < T->q || T->p > T->b) { T->p = 0; break; } /* overflow */
     983     2798445 :     } while (!uisprime(T->p));
     984     1049932 :     if (T->p && T->p <= T->b) return T->p;
     985             :     /* overflow ulong, switch to GEN */
     986          52 :     T->strategy = PRST_nextprime;
     987             :   }
     988        1096 :   return 0; /* overflow */
     989             : }
     990             : 
     991             : GEN
     992    45385959 : forprime_next(forprime_t *T)
     993             : {
     994             :   pari_sp av;
     995             :   GEN p;
     996    45385959 :   if (T->strategy != PRST_nextprime)
     997             :   {
     998    45378080 :     ulong u = u_forprime_next(T);
     999    45378080 :     if (u) { affui(u, T->pp); return T->pp; }
    1000             :     /* failure */
    1001         814 :     if (T->strategy != PRST_nextprime) return NULL; /* we're done */
    1002             :     /* overflow ulong, switch to GEN */
    1003          48 :     u = ULONG_MAX;
    1004          48 :     if (T->q > 1) u -= (ULONG_MAX-T->c) % T->q;
    1005          48 :     affui(u, T->pp);
    1006             :   }
    1007        7927 :   av = avma; p = T->pp;
    1008        7927 :   if (T->q == 1)
    1009             :   {
    1010        7749 :     p = nextprime(addiu(p, 1));
    1011        7749 :     if (T->bb && abscmpii(p, T->bb) > 0) return gc_NULL(av);
    1012             :   } else do {
    1013        3341 :     p = T->qq? addii(p, T->qq): addiu(p, T->q);
    1014        3341 :     if (T->bb && abscmpii(p, T->bb) > 0) return gc_NULL(av);
    1015        3285 :   } while (!BPSW_psp(p));
    1016        7744 :   affii(p, T->pp); return gc_const(av, T->pp);
    1017             : }
    1018             : 
    1019             : void
    1020        1085 : forprimestep(GEN a, GEN b, GEN q, GEN code)
    1021             : {
    1022        1085 :   pari_sp av = avma;
    1023             :   forprime_t T;
    1024             : 
    1025        1085 :   if (!forprimestep_init(&T, a,b,q)) { set_avma(av); return; }
    1026             : 
    1027        1071 :   push_lex(T.pp,code);
    1028      309822 :   while(forprime_next(&T))
    1029             :   {
    1030      309178 :     closure_evalvoid(code); if (loop_break()) break;
    1031             :     /* p changed in 'code', complain */
    1032      308758 :     if (get_lex(-1) != T.pp)
    1033           7 :       pari_err(e_MISC,"prime index read-only: was changed to %Ps", get_lex(-1));
    1034             :   }
    1035        1064 :   pop_lex(1); set_avma(av);
    1036             : }
    1037             : void
    1038         959 : forprime(GEN a, GEN b, GEN code) { return forprimestep(a,b,NULL,code); }
    1039             : 
    1040             : int
    1041          70 : forcomposite_init(forcomposite_t *C, GEN a, GEN b)
    1042             : {
    1043          70 :   pari_sp av = avma;
    1044          70 :   a = gceil(a);
    1045          70 :   if (typ(a)!=t_INT) pari_err_TYPE("forcomposite",a);
    1046          70 :   if (b) {
    1047          63 :     if (typ(b) == t_INFINITY) b = NULL;
    1048             :     else
    1049             :     {
    1050          56 :       b = gfloor(b);
    1051          56 :       if (typ(b)!=t_INT) pari_err_TYPE("forcomposite",b);
    1052             :     }
    1053             :   }
    1054          70 :   if (signe(a) < 0) pari_err_DOMAIN("forcomposite", "a", "<", gen_0, a);
    1055          70 :   if (abscmpiu(a, 4) < 0) a = utoipos(4);
    1056          70 :   C->first = 1;
    1057          70 :   if (!forprime_init(&C->T, a,b) && cmpii(a,b) > 0)
    1058             :   {
    1059           7 :     C->n = gen_1; /* in case caller forgets to check the return value */
    1060           7 :     C->b = gen_0; return gc_bool(av,0);
    1061             :   }
    1062          63 :   C->n = setloop(a);
    1063          63 :   C->b = b;
    1064          63 :   C->p = NULL; return 1;
    1065             : }
    1066             : 
    1067             : GEN
    1068         238 : forcomposite_next(forcomposite_t *C)
    1069             : {
    1070         238 :   if (C->first) /* first call ever */
    1071             :   {
    1072          63 :     C->first = 0;
    1073          63 :     C->p = forprime_next(&C->T);
    1074             :   }
    1075             :   else
    1076         175 :     C->n = incloop(C->n);
    1077         238 :   if (C->p)
    1078             :   {
    1079         161 :     if (cmpii(C->n, C->p) < 0) return C->n;
    1080          77 :     C->n = incloop(C->n);
    1081             :     /* n = p+1 */
    1082          77 :     C->p = forprime_next(&C->T); /* nextprime(p) > n */
    1083          77 :     if (C->p) return C->n;
    1084             :   }
    1085         105 :   if (!C->b || cmpii(C->n, C->b) <= 0) return C->n;
    1086          42 :   return NULL;
    1087             : }
    1088             : 
    1089             : void
    1090          70 : forcomposite(GEN a, GEN b, GEN code)
    1091             : {
    1092          70 :   pari_sp av = avma;
    1093             :   forcomposite_t T;
    1094             :   GEN n;
    1095          70 :   if (!forcomposite_init(&T,a,b)) return;
    1096          63 :   push_lex(T.n,code);
    1097         238 :   while((n = forcomposite_next(&T)))
    1098             :   {
    1099         196 :     closure_evalvoid(code); if (loop_break()) break;
    1100             :     /* n changed in 'code', complain */
    1101         182 :     if (get_lex(-1) != n)
    1102           7 :       pari_err(e_MISC,"index read-only: was changed to %Ps", get_lex(-1));
    1103             :   }
    1104          56 :   pop_lex(1); set_avma(av);
    1105             : }

Generated by: LCOV version 1.16