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, _prodprimeslim;
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 1887 : initprimes1(ulong size, long *lenp, pari_prime *p1)
34 : {
35 1887 : pari_sp av = avma;
36 : long k;
37 1887 : byteptr q, r, s, p = (byteptr)stack_calloc(size+2), fin = p + size;
38 : pari_prime *re;
39 :
40 22644 : for (r=q=p,k=1; r<=fin; )
41 : {
42 33966 : do { r+=k; k+=2; r+=k; } while (*++q);
43 911421 : for (s=r; s<=fin; s+=k) *s = 1;
44 : }
45 1887 : re = p1; *re++ = 2; *re++ = 3; /* 2 and 3 */
46 1887 : for (s=q=p+1; ; s=q)
47 : {
48 964257 : do q++; while (*q);
49 322677 : if (q > fin) break;
50 320790 : *re++ = (pari_prime) 2*(q-p)+1;
51 : }
52 1887 : *re++ = 0;
53 1887 : *lenp = re - p1;
54 1887 : set_avma(av);
55 1887 : }
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 1887 : good_arena_size(ulong slow2_size, ulong total, ulong fixed_to_cache,
201 : cache_model_t *cache_model)
202 : {
203 1887 : ulong asize, cache_arena = cache_model->arena;
204 : double Xmin, Xmax, A, B, C1, C2, D, V;
205 1887 : 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 1887 : asize = cache_arena - fixed_to_cache;
236 1887 : if (total <= asize) return total;
237 : /* The simple case: fitting into cache doesn't slow us down more than 10% */
238 1887 : 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 7251 : sieve_chunk(pari_prime *known_primes, ulong s, byteptr data, ulong n)
343 : {
344 7251 : ulong p, cnt = n-1, start = s;
345 : pari_prime *q;
346 :
347 7251 : memset(data, 0, n);
348 7251 : start >>= 1; /* (start - 1)/2 */
349 7251 : start += n; /* Corresponds to the end */
350 : /* data corresponds to start, q runs over primediffs */
351 1041348 : 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 1034097 : long off = cnt - ((start+(p>>1)) % p);
358 1647201847 : while (off >= 0) { data[off] = 1; off -= p; }
359 : }
360 7251 : }
361 :
362 : static long
363 1887 : ZV_size(GEN x)
364 : {
365 1887 : long i, l = lg(x), s = l;
366 26418 : for (i = 1; i < l; i++) s += lgefint(gel(x,i));
367 1887 : return s;
368 : }
369 : /* avoid gcopy_avma so as to deallocate using free() in pari_close_primes */
370 : static GEN
371 1887 : ZV_copy_alloc(GEN x)
372 : {
373 1887 : long i, l = lg(x), s = ZV_size(x);
374 1887 : GEN z = (GEN)pari_malloc(s * sizeof(long)), av = z + s;
375 26418 : for (i = 1; i < l; i++) gel(z,i) = av = icopy_avma(gel(x,i), (pari_sp)av);
376 1887 : z[0] = x[0] & (TYPBITS|LGBITS); return z;
377 : }
378 :
379 : static void
380 1887 : set_prodprimes(ulong M)
381 : {
382 1887 : pari_sp av = avma;
383 1887 : ulong b = 1UL << 8;
384 1887 : GEN W, v = primes_interval_zv(3, M);
385 1887 : long u, j, jold, l = lg(v);
386 :
387 1887 : W = cgetg(64+1, t_VEC);
388 1887 : M = v[l-1]; /* = precprime(M) */
389 1887 : _prodprimeslim = cgetalloc(64+1, t_VECSMALL);
390 1887 : if (b > M) b = M;
391 154781175 : for (jold = j = u = 1; j < l; j++)
392 154779288 : if (uel(v,j) >= b) /* if j = l-1, then b = M = v[j] */
393 : {
394 24531 : long lw = (j == l-1? l: j) - jold + 1;
395 24531 : GEN w = v+jold-1;
396 24531 : w[0] = evaltyp(t_VECSMALL) | _evallg(lw);
397 24531 : _prodprimeslim[u] = w[lw - 1];
398 : /* product of primes from
399 : * p_jold = 3 if first entry, else nextprime(_prodprime_lim[u - 1] + 1)
400 : * to
401 : * p_{j-1} = _prodprimeslim[u] = precprime(M or 2^{u+7}) */
402 24531 : gel(W,u++) = zv_prod_Z(w);
403 24531 : jold = j; b *= 2;
404 24531 : if (b > M) b = M; /* truncate last run */
405 : }
406 24531 : for (j = 2; j < u; j++) gel(W,j) = mulii(gel(W,j-1), gel(W,j));
407 1887 : setlg(W, u); _prodprimes = ZV_copy_alloc(W);
408 1887 : setlg(_prodprimeslim, u); set_avma(av);
409 1887 : }
410 :
411 : static void
412 1887 : initprimes0(ulong maxnum, long *lenp, pari_prime *p1)
413 : {
414 1887 : pari_sp av = avma, bot = pari_mainstack->bot;
415 : long alloced, psize;
416 : byteptr q, end, p;
417 : ulong remains, curlow, rootnum, asize, prime_above, last;
418 : pari_prime *end1, *curdiff, *p_prime_above;
419 :
420 1887 : if (!odd(maxnum)) maxnum--; /* make it odd. */
421 : /* base case */
422 1887 : if (maxnum < 1ul<<17) { initprimes1(maxnum>>1, lenp, p1); return; }
423 :
424 : /* Checked to be enough up to 40e6, attained at 155893 */
425 1887 : rootnum = usqrt(maxnum) | 1;
426 1887 : initprimes1(rootnum>>1, &psize, p1);
427 1887 : last = rootnum;
428 1887 : end1 = p1 + psize - 1;
429 1887 : remains = (maxnum - last) >> 1; /* number of odd numbers to check */
430 : /* we access primes array of psize too; but we access it consecutively,
431 : * thus we do not include it in fixed_to_cache */
432 1887 : asize = good_arena_size((ulong)(rootnum * slow2_in_roots), remains+1, 0,
433 : &cache_model) - 1;
434 : /* enough room on the stack ? */
435 1887 : alloced = (((byteptr)avma) <= ((byteptr)bot) + asize);
436 1887 : p = (byteptr)(alloced? pari_malloc(asize+1): stack_malloc(asize+1));
437 1887 : end = p + asize; /* the 0 sentinel goes at end. */
438 1887 : curlow = last + 2; /* First candidate: know primes up to last (odd). */
439 1887 : curdiff = end1;
440 :
441 : /* During each iteration p..end-1 represents a range of odd
442 : numbers. */
443 1887 : p_prime_above = p1 + 2;
444 1887 : prime_above = 3;
445 9138 : while (remains)
446 : { /* cycle over arenas; performance not crucial */
447 : pari_prime was_delta;
448 7251 : if (asize > remains) { asize = remains; end = p + asize; }
449 : /* Fake the upper limit appropriate for the given arena */
450 329928 : while (prime_above*prime_above <= curlow + (asize << 1) && *p_prime_above)
451 322677 : prime_above = *p_prime_above++;
452 7251 : was_delta = *p_prime_above;
453 7251 : *p_prime_above = 0; /* sentinel for sieve_chunk */
454 7251 : sieve_chunk(p1, curlow, p, asize);
455 7251 : *p_prime_above = was_delta; /* restore */
456 :
457 7251 : p[asize] = 0; /* sentinel */
458 7251 : for (q = p; ; q++)
459 : { /* q runs over addresses corresponding to primes */
460 988372563 : while (*q) q++; /* use sentinel at end */
461 154463862 : if (q >= end) break;
462 154456611 : *curdiff++ = (pari_prime) 2*(q-p) + curlow;
463 : }
464 7251 : remains -= asize;
465 7251 : curlow += (asize<<1);
466 : }
467 1887 : *curdiff++ = 0; /* sentinel */
468 1887 : *lenp = curdiff - p1;
469 1887 : if (alloced) pari_free(p); else set_avma(av);
470 : }
471 :
472 : ulong
473 36499206 : maxprime(void) { return pari_PRIMES? pari_PRIMES[pari_PRIMES[0]]: 0; }
474 : ulong
475 64097379 : maxprimelim(void) { return pari_PRIMES? _maxprimelim: 0; }
476 : ulong
477 196 : maxprimeN(void) { return pari_PRIMES? pari_PRIMES[0]: 0; }
478 : GEN
479 1106514 : prodprimes(void) { return pari_PRIMES? _prodprimes: NULL; }
480 : GEN
481 1106514 : prodprimeslim(void) { return pari_PRIMES? _prodprimeslim: NULL; }
482 : void
483 0 : maxprime_check(ulong c) { if (maxprime() < c) pari_err_MAXPRIME(c); }
484 :
485 : static pari_prime*
486 1887 : initprimes(ulong maxnum)
487 : {
488 : pari_prime *t;
489 : long len;
490 : ulong N;
491 1887 : if (maxnum < 65537)
492 : {
493 0 : maxnum = 65537;
494 0 : N = 6543;
495 : }
496 : else
497 1887 : N = (long) ceil(primepi_upper_bound((double)maxnum));
498 1887 : t = (pari_prime*) pari_malloc(sizeof(*t) * (N+2));
499 1887 : initprimes0(maxnum, &len, t+1); t[0] = (pari_prime)(len-1);
500 1887 : _maxprimelim = maxnum;
501 1887 : return (pari_prime*) pari_realloc(t, sizeof(*t) * (len+1));
502 : }
503 :
504 : void
505 1887 : initprimetable(ulong maxnum, ulong factorlimit)
506 : {
507 1887 : pari_prime *old = pari_PRIMES;
508 : #ifdef LONG_IS_64BIT
509 1622 : maxnum = minuu(maxnum, 4294967295);
510 : #endif
511 1887 : pari_PRIMES = initprimes(maxnum);
512 1887 : if (old) free(old);
513 1887 : set_prodprimes(minuu(maxprime(), factorlimit));
514 1887 : }
515 :
516 : /**********************************************************************/
517 : /*** ***/
518 : /*** forprime ***/
519 : /*** ***/
520 : /**********************************************************************/
521 :
522 : /* return good chunk size for sieve, 16 | chunk + 2 */
523 : static ulong
524 8843446 : optimize_chunk(ulong a, ulong b)
525 : {
526 : /* TODO: Optimize size (surely < 512k to stay in L2 cache, but not so large
527 : * as to force recalculating too often). */
528 : /* bigarena is in bytes, we use bits, and only odds */
529 8843446 : ulong defchunk = (a>>31) > 1 ? 0x80000UL : 0x8000;
530 8843446 : ulong chunk = (cache_model.bigarena ? cache_model.bigarena : defchunk)<<4;
531 8843446 : ulong tmp = (b - a) / chunk + 1;
532 :
533 8843446 : if (tmp == 1)
534 195 : chunk = b - a + 16;
535 : else
536 8843251 : chunk = (b - a) / tmp + 15;
537 : /* ensure 16 | chunk + 2 */
538 8843446 : return (((chunk + 2)>>4)<<4) - 2;
539 : }
540 : static void
541 8843446 : sieve_init(forprime_t *T, ulong a, ulong b)
542 : {
543 8843446 : T->sieveb = b;
544 8843446 : T->chunk = optimize_chunk(a, b);
545 : /* >> 1 [only odds] + 3 [convert from bits to bytes] */
546 8843446 : T->isieve = (unsigned char*)stack_malloc(((T->chunk+2) >> 4) + 1);
547 8843446 : T->cache[0] = 0;
548 8843446 : T->a = a;
549 8843446 : T->end = minuu(a + T->chunk, b);
550 8843446 : T->pos = T->maxpos = 0;
551 8843446 : }
552 :
553 : enum {PRST_none, PRST_diffptr, PRST_sieve, PRST_unextprime, PRST_nextprime};
554 :
555 : static void
556 8843575 : u_forprime_set_prime_table(forprime_t *T, ulong a)
557 : {
558 8843575 : T->strategy = PRST_diffptr;
559 8843575 : if (a < 3)
560 : {
561 2785147 : T->p = 0;
562 2785147 : T->n = 0;
563 : }
564 : else
565 : {
566 6058428 : long n = PRIMES_search(a - 1);
567 6058428 : if (n < 0) n = - n - 1;
568 6058428 : T->n = n;
569 6058428 : T->p = pari_PRIMES[n];
570 : }
571 8843575 : }
572 :
573 : /* Set p so that p + q the smallest integer = c (mod q) and > original p.
574 : * Assume 0 < c < q. */
575 : static void
576 102080 : arith_set(forprime_t *T)
577 : {
578 102080 : ulong r = T->p % T->q; /* 0 <= r <= min(p, q-1) */
579 102080 : pari_sp av = avma;
580 102080 : GEN d = adduu(T->p - r, T->c); /* = c mod q */
581 102080 : if (T->c > r) d = subiu(d, T->q);
582 : /* d = c mod q, d = c > r? p-r+c-q: p-r+c, so that
583 : * d <= p and d+q = c>r? p-r+c : p-r+c+q > p */
584 102080 : if (signe(d) <= 0)
585 : {
586 20 : T->p = 0;
587 20 : T->strategy = PRST_nextprime;
588 20 : affii(d, T->pp);
589 : }
590 : else
591 102060 : T->p = itou_or_0(d);
592 102080 : set_avma(av);
593 102080 : }
594 :
595 : /* Run through primes in arithmetic progression = c (mod q).
596 : * Warning: b = ULONG_MAX may signal that we are called by higher level
597 : * function handling a continuation for larger b; this sentinel value
598 : * must not be modified */
599 : static int
600 25530517 : u_forprime_sieve_arith_init(forprime_t *T, struct pari_sieve *psieve,
601 : ulong a, ulong b, ulong c, ulong q)
602 : {
603 : #ifdef LONG_IS_64BIT
604 21902746 : const ulong UPRIME_MAX = 18446744073709551557UL;
605 : #else
606 3627771 : const ulong UPRIME_MAX = 4294967291UL;
607 : #endif
608 : ulong Plim, P, P2, Y, sieveb;
609 :
610 25530517 : if (!odd(b) && b > 2) b--;
611 25530517 : if (a > b || b < 2)
612 : {
613 883670 : T->strategy = PRST_diffptr; /* paranoia */
614 883670 : T->p = 0; /* empty */
615 883670 : T->b = 0; /* empty */
616 883670 : T->n = 0;
617 883670 : return 0;
618 : }
619 24646847 : P = maxprime();
620 24646847 : if (b != ULONG_MAX && b > UPRIME_MAX) b = UPRIME_MAX;
621 24646847 : if (q != 1)
622 : {
623 : ulong D;
624 544245 : c %= q; D = ugcd(c, q);
625 544245 : if (D != 1) { a = maxuu(a,D); if (b != ULONG_MAX) b = minuu(b,D); }
626 544245 : if (odd(q) && (a > 2 || c != 2))
627 : { /* only *odd* primes. If a <= c = 2, then p = 2 must be included :-( */
628 459587 : if (!odd(c)) c += q;
629 459587 : q <<= 1;
630 : }
631 : }
632 24646847 : T->q = q;
633 24646847 : T->c = c;
634 24646847 : T->strategy = PRST_none; /* unknown */
635 24646847 : T->psieve = psieve; /* unused for now */
636 24646847 : T->isieve = NULL; /* unused for now */
637 24646847 : T->b = b;
638 24646847 : if (P >= b) { /* [a,b] \subset prime table */
639 4944713 : u_forprime_set_prime_table(T, a);
640 4944713 : return 1;
641 : }
642 : /* b > P */
643 19702134 : if (a >= P)
644 : {
645 15803272 : T->p = a - 1;
646 15803272 : if (T->q != 1) arith_set(T);
647 : }
648 : else
649 3898862 : u_forprime_set_prime_table(T, a);
650 19702134 : if (T->strategy == PRST_none) T->strategy = PRST_unextprime;
651 : /* now strategy is either PRST_diffptr or PRST_unextprime */
652 :
653 19702134 : P2 = (P & HIGHMASK)? 0 : P*P;
654 19702134 : sieveb = b; if (P2 && P2 < b) sieveb = P2;
655 : /* maxprime^2 >= sieveb */
656 19702134 : Plim = maxprimelim();
657 19702134 : if (a <= Plim) a = Plim + 1;
658 19702134 : if (sieveb < a + 16) return 1;
659 9302196 : Y = sieveb - a + 1; /* number of integers in sievable interval > 16 */
660 9302196 : P = usqrt(sieveb); /* largest sieving prime */
661 : /* FIXME: should sieve as well if q != 1, adapt sieve code */
662 9302196 : if (q == 1 && (!P2 || P2 > a) && 3/M_LN2 * Y >= uprimepi(P))
663 : /* Sieve implemented & possible & not too costly. Cost model is
664 : * - nextprime: about Y / log(b) primes to test [neglect cost for composites]
665 : * individual cost average = 3 log2(b) mulmod, total = 3 Y / log(2) mulmod
666 : * - sieve: pi(P) mod + Y loglog(b) add
667 : * Since loglog(b) < 4, and add < 10*mulmod, we neglect the Y loglog(b) term.
668 : * We have mod < mulmod < 2*mod; for now, assume mulmod ~ mod. */
669 : {
670 8843446 : if (T->strategy == PRST_unextprime) T->strategy = PRST_sieve;
671 8843446 : sieve_init(T, a, sieveb);
672 : }
673 9302196 : return 1;
674 : }
675 :
676 : int
677 20136090 : u_forprime_arith_init(forprime_t *T, ulong a, ulong b, ulong c, ulong q)
678 20136090 : { return u_forprime_sieve_arith_init(T, NULL, a, b, c, q); }
679 :
680 : /* will run through primes in [a,b] */
681 : int
682 19589237 : u_forprime_init(forprime_t *T, ulong a, ulong b)
683 19589237 : { return u_forprime_arith_init(T, a,b, 0,1); }
684 :
685 : /* will run through primes in [a,b] */
686 : static int
687 5394427 : u_forprime_sieve_init(forprime_t *T, struct pari_sieve *s, ulong b)
688 5394427 : { return u_forprime_sieve_arith_init(T, s, s->start, b, s->c, s->q); }
689 :
690 : /* now only run through primes <= c; assume c <= b above */
691 : void
692 63 : u_forprime_restrict(forprime_t *T, ulong c) { T->b = c; }
693 :
694 : /* b = NULL: loop forever */
695 : int
696 2791 : forprimestep_init(forprime_t *T, GEN a, GEN b, GEN q)
697 : {
698 2791 : GEN c = NULL;
699 : long lb;
700 :
701 2791 : a = gceil(a); if (typ(a) != t_INT) pari_err_TYPE("forprime_init",a);
702 2791 : T->qq = NULL; T->q = 1; T->c = 0;
703 2791 : if (q)
704 : {
705 133 : switch(typ(q))
706 : {
707 56 : case t_INT:
708 56 : c = a; break;
709 77 : case t_INTMOD:
710 77 : c = gel(q,2); q = gel(q,1);
711 : /* first int >= initial a which is = c (mod q) */
712 77 : a = addii(a, modii(subii(c,a), q)); break;
713 0 : default: pari_err_TYPE("forprimestep_init",q);
714 : }
715 133 : if (signe(q) <= 0) pari_err_TYPE("forprimestep_init (q <= 0)",q);
716 133 : if (equali1(q)) c = q = NULL;
717 : else
718 : {
719 133 : GEN D = gcdii(c, q);
720 133 : if (!is_pm1(D))
721 : { /* at most one prime: c */
722 42 : if (cmpii(a, D) < 0) a = D;
723 42 : if (gcmp(b, D) > 0) b = D;
724 : }
725 133 : if ((T->q = itou_or_0(q)))
726 125 : T->c = umodiu(c, T->q);
727 : else
728 8 : T->qq = q;
729 : }
730 : }
731 2791 : if (signe(a) <= 0) a = q? modii(a, q): gen_1;
732 2791 : if (b && typ(b) != t_INFINITY)
733 : {
734 1384 : b = gfloor(b);
735 1384 : if (typ(b) != t_INT) pari_err_TYPE("forprime_init",b);
736 1384 : if (signe(b) < 0 || cmpii(a,b) > 0)
737 : {
738 21 : T->strategy = PRST_nextprime; /* paranoia */
739 21 : T->bb = T->pp = gen_0; return 0;
740 : }
741 1363 : lb = lgefint(b);
742 1363 : T->bb = b;
743 : }
744 1407 : else if (!b || inf_get_sign(b) > 0)
745 : {
746 1407 : lb = lgefint(a) + 4;
747 1407 : T->bb = NULL;
748 : }
749 : else /* b == -oo */
750 : {
751 0 : T->strategy = PRST_nextprime; /* paranoia */
752 0 : T->bb = T->pp = gen_0; return 0;
753 : }
754 2770 : T->pp = cgeti(T->qq? maxuu(lb, lgefint(T->qq)): lb);
755 : /* a, b are positive integers, a <= b */
756 2770 : if (!T->qq && lgefint(a) == 3) /* lb == 3 implies b != NULL */
757 2626 : return u_forprime_arith_init(T, uel(a,2), lb == 3? uel(b,2): ULONG_MAX,
758 : T->c, T->q);
759 144 : T->strategy = PRST_nextprime;
760 144 : affii(T->qq? subii(a,T->qq): subiu(a,T->q), T->pp); return 1;
761 : }
762 : int
763 1315 : forprime_init(forprime_t *T, GEN a, GEN b)
764 1315 : { return forprimestep_init(T,a,b,NULL); }
765 :
766 : /* assume a <= b <= maxprime()^2, a,b odd, sieve[n] corresponds to
767 : * a+16*n, a+16*n+2, ..., a+16*n+14 (bits 0 to 7)
768 : * maxpos = index of last sieve cell.
769 : * b-a+2 must be divisible by 16 for use by u_forprime_next */
770 : static void
771 9238 : sieve_block(ulong a, ulong b, ulong maxpos, unsigned char* sieve)
772 : {
773 9238 : ulong i, N = pari_PRIMES[0], lim = usqrt(b), sz = (b-a) >> 1;
774 9238 : (void)memset(sieve, 0, maxpos+1);
775 25143422 : for (i = 2; i <= N; i++)
776 : { /* p is odd */
777 25143416 : ulong k, r, p = pari_PRIMES[i]; /* starts at p = 3 */
778 25143416 : if (p > lim) break;
779 :
780 : /* solve a + 2k = 0 (mod p) */
781 25134184 : r = a % p;
782 25134184 : if (r == 0)
783 16223 : k = 0;
784 : else
785 : {
786 25117961 : k = p - r;
787 25117961 : if (odd(k)) k += p;
788 25117961 : k >>= 1;
789 : }
790 : /* m = a + 2k is the smallest odd m >= a, p | m */
791 : /* position n (corresponds to a+2n) is sieve[n>>3], bit n&7 */
792 5796635715 : while (k <= sz) { sieve[k>>3] |= 1 << (k&7); k += p; /* 2k += 2p */ }
793 : }
794 9238 : }
795 :
796 : static void
797 1887 : pari_sieve_init(struct pari_sieve *s, ulong a, ulong b)
798 : {
799 1887 : ulong maxpos= (b - a) >> 4;
800 1887 : s->start = a; s->end = b;
801 1887 : s->sieve = (unsigned char*) pari_malloc(maxpos+1);
802 1887 : s->c = 0; s->q = 1;
803 1887 : sieve_block(a, b, maxpos, s->sieve);
804 1887 : s->maxpos = maxpos; /* must be last in case of SIGINT */
805 1887 : }
806 :
807 : static struct pari_sieve pari_sieve_modular;
808 : /* p and f are allowed to be 0 (default) */
809 : static void
810 1887 : synchronize_primelimits(ulong p, ulong f)
811 : {
812 1887 : ulong DFT = 1UL << 20;
813 1887 : if (!p && !f) p = f = DFT;
814 0 : else if (p <= f) p = f; /* includes p = 0 */
815 0 : else if (!f) f = minuu(p, DFT);
816 1887 : GP_DATA->primelimit = p;
817 1887 : GP_DATA->factorlimit = f; /* <= primelimit */
818 1887 : }
819 :
820 : /* also initializes GP_DATA->(primelimit | factorlimit) */
821 : void
822 1887 : pari_init_primes(ulong maxprime, ulong factorlimit)
823 : {
824 1887 : ulong a = (1UL<<31) + 1, b = a + (1UL<<20)-2;
825 1887 : synchronize_primelimits(maxprime, factorlimit);
826 1887 : initprimetable(GP_DATA->primelimit, GP_DATA->factorlimit);
827 1887 : pari_sieve_init(&pari_sieve_modular, a, b);
828 1887 : }
829 :
830 : void
831 1887 : pari_close_primes(void)
832 : {
833 1887 : if (pari_PRIMES)
834 : {
835 1887 : pari_free(pari_PRIMES);
836 1887 : pari_free(_prodprimes);
837 1887 : pari_free(_prodprimeslim);
838 : }
839 1887 : pari_free(pari_sieve_modular.sieve);
840 1887 : }
841 :
842 : void
843 4278320 : init_modular_small(forprime_t *S)
844 : {
845 : #ifdef LONG_IS_64BIT
846 3677246 : u_forprime_sieve_init(S, &pari_sieve_modular, ULONG_MAX);
847 : #else
848 601074 : ulong a = (1UL<<((BITS_IN_LONG-2)>>1))+1;
849 601074 : u_forprime_init(S, a, ULONG_MAX);
850 : #endif
851 4278320 : }
852 :
853 : void
854 12020644 : init_modular_big(forprime_t *S)
855 : {
856 : #ifdef LONG_IS_64BIT
857 10303463 : u_forprime_init(S, HIGHBIT + 1, ULONG_MAX);
858 : #else
859 1717181 : u_forprime_sieve_init(S, &pari_sieve_modular, ULONG_MAX);
860 : #endif
861 12020644 : }
862 :
863 : /* T->cache is a 0-terminated list of primes, return the first one and
864 : * remove it from list. Most of the time the list contains a single prime */
865 : static ulong
866 132523313 : shift_cache(forprime_t *T)
867 : {
868 : long i;
869 132523313 : T->p = T->cache[0];
870 177296710 : for (i = 1;; i++) /* remove one prime from cache */
871 177296710 : if (! (T->cache[i-1] = T->cache[i]) ) break;
872 132523313 : return T->p;
873 : }
874 :
875 : ulong
876 219640779 : u_forprime_next(forprime_t *T)
877 : {
878 219640779 : if (T->strategy == PRST_diffptr)
879 : {
880 : for(;;)
881 : {
882 228419315 : if (++T->n <= pari_PRIMES[0])
883 : {
884 228419154 : T->p = pari_PRIMES[T->n];
885 228419154 : if (T->p > T->b) return 0;
886 228217006 : if (T->q == 1 || T->p % T->q == T->c) return T->p;
887 : }
888 : else
889 : { /* beyond the table */
890 161 : T->strategy = T->isieve? PRST_sieve: PRST_unextprime;
891 161 : if (T->q != 1) { arith_set(T); if (!T->p) return 0; }
892 : /* T->p possibly not a prime if q != 1 */
893 161 : break;
894 : }
895 : }
896 : }
897 148873793 : if (T->strategy == PRST_sieve)
898 : { /* require sieveb - a >= 16 */
899 : ulong n;
900 132523738 : if (T->cache[0]) return shift_cache(T);
901 94703926 : NEXT_CHUNK:
902 94711277 : if (T->psieve)
903 : {
904 5394423 : T->sieve = T->psieve->sieve;
905 5394423 : T->end = T->psieve->end;
906 5394423 : if (T->end > T->sieveb) T->end = T->sieveb;
907 5394423 : T->maxpos = T->psieve->maxpos;
908 5394423 : T->pos = 0;
909 5394423 : T->psieve = NULL;
910 : }
911 143243326 : for (n = T->pos; n < T->maxpos; n++)
912 143232728 : if (T->sieve[n] != 0xFF)
913 : {
914 94700679 : unsigned char mask = T->sieve[n];
915 94700679 : ulong p = T->a + (n<<4);
916 94700679 : long i = 0;
917 94700679 : T->pos = n;
918 94700679 : if (!(mask & 1)) T->cache[i++] = p;
919 94700679 : if (!(mask & 2)) T->cache[i++] = p+2;
920 94700679 : if (!(mask & 4)) T->cache[i++] = p+4;
921 94700679 : if (!(mask & 8)) T->cache[i++] = p+6;
922 94700679 : if (!(mask & 16)) T->cache[i++] = p+8;
923 94700679 : if (!(mask & 32)) T->cache[i++] = p+10;
924 94700679 : if (!(mask & 64)) T->cache[i++] = p+12;
925 94700679 : if (!(mask &128)) T->cache[i++] = p+14;
926 94700679 : T->cache[i] = 0;
927 94700679 : T->pos = n+1;
928 94700679 : return shift_cache(T);
929 : }
930 : /* n = T->maxpos, last cell: check p <= b */
931 10598 : if (T->maxpos && n == T->maxpos && T->sieve[n] != 0xFF)
932 : {
933 2930 : unsigned char mask = T->sieve[n];
934 2930 : ulong p = T->a + (n<<4);
935 2930 : long i = 0;
936 2930 : T->pos = n;
937 2930 : if (!(mask & 1) && p <= T->sieveb) T->cache[i++] = p;
938 2930 : if (!(mask & 2) && p <= T->sieveb-2) T->cache[i++] = p+2;
939 2930 : if (!(mask & 4) && p <= T->sieveb-4) T->cache[i++] = p+4;
940 2930 : if (!(mask & 8) && p <= T->sieveb-6) T->cache[i++] = p+6;
941 2930 : if (!(mask & 16) && p <= T->sieveb-8) T->cache[i++] = p+8;
942 2930 : if (!(mask & 32) && p <= T->sieveb-10) T->cache[i++] = p+10;
943 2930 : if (!(mask & 64) && p <= T->sieveb-12) T->cache[i++] = p+12;
944 2930 : if (!(mask &128) && p <= T->sieveb-14) T->cache[i++] = p+14;
945 2930 : if (i)
946 : {
947 2822 : T->cache[i] = 0;
948 2822 : T->pos = n+1;
949 2822 : return shift_cache(T);
950 : }
951 : }
952 :
953 7776 : if (T->maxpos && T->end >= T->sieveb) /* done with sieves ? */
954 : {
955 425 : if (T->sieveb == T->b && T->b != ULONG_MAX) return 0;
956 7 : T->strategy = PRST_unextprime;
957 : }
958 : else
959 : { /* initialize next chunk */
960 7351 : T->sieve = T->isieve;
961 7351 : if (T->maxpos == 0)
962 3379 : T->a |= 1; /* first time; ensure odd */
963 : else
964 3972 : T->a = (T->end + 2) | 1;
965 7351 : T->end = T->a + T->chunk; /* may overflow */
966 7351 : if (T->end < T->a || T->end > T->sieveb) T->end = T->sieveb;
967 : /* end and a are odd; sieve[k] contains the a + 8*2k + (0,2,...,14).
968 : * The largest k is (end-a) >> 4 */
969 7351 : T->pos = 0;
970 7351 : T->maxpos = (T->end - T->a) >> 4; /* >= 1 */
971 7351 : sieve_block(T->a, T->end, T->maxpos, T->sieve);
972 7351 : goto NEXT_CHUNK;
973 : }
974 : }
975 16350062 : if (T->strategy == PRST_unextprime)
976 : {
977 16350062 : if (T->q == 1)
978 : {
979 : #ifdef LONG_IS_64BIT
980 16196526 : switch(T->p)
981 : {
982 : #define retp(x) return T->p = (HIGHBIT+x <= T->b)? HIGHBIT+x: 0
983 10303439 : case HIGHBIT: retp(29);
984 3199949 : case HIGHBIT + 29: retp(99);
985 344278 : case HIGHBIT + 99: retp(123);
986 189037 : case HIGHBIT +123: retp(131);
987 153158 : case HIGHBIT +131: retp(155);
988 135726 : case HIGHBIT +155: retp(255);
989 115105 : case HIGHBIT +255: retp(269);
990 105106 : case HIGHBIT +269: retp(359);
991 85021 : case HIGHBIT +359: retp(435);
992 65461 : case HIGHBIT +435: retp(449);
993 58511 : case HIGHBIT +449: retp(453);
994 55529 : case HIGHBIT +453: retp(485);
995 49423 : case HIGHBIT +485: retp(491);
996 46246 : case HIGHBIT +491: retp(543);
997 43840 : case HIGHBIT +543: retp(585);
998 41213 : case HIGHBIT +585: retp(599);
999 34916 : case HIGHBIT +599: retp(753);
1000 33968 : case HIGHBIT +753: retp(849);
1001 32504 : case HIGHBIT +849: retp(879);
1002 30680 : case HIGHBIT +879: retp(885);
1003 29774 : case HIGHBIT +885: retp(903);
1004 29216 : case HIGHBIT +903: retp(995);
1005 : #undef retp
1006 : }
1007 : #endif
1008 1014498 : T->p = unextprime(T->p + 1);
1009 1014498 : if (T->p > T->b) return 0;
1010 : }
1011 : else do {
1012 2802903 : T->p += T->q;
1013 2802903 : if (T->p < T->q || T->p > T->b) { T->p = 0; break; } /* overflow */
1014 2802877 : } while (!uisprime(T->p));
1015 1167899 : if (T->p && T->p <= T->b) return T->p;
1016 : /* overflow ulong, switch to GEN */
1017 48 : T->strategy = PRST_nextprime;
1018 : }
1019 48 : return 0; /* overflow */
1020 : }
1021 :
1022 : GEN
1023 47789149 : forprime_next(forprime_t *T)
1024 : {
1025 : pari_sp av;
1026 : GEN p;
1027 47789149 : if (T->strategy != PRST_nextprime)
1028 : {
1029 47777350 : ulong u = u_forprime_next(T);
1030 47777350 : if (u) { affui(u, T->pp); return T->pp; }
1031 : /* failure */
1032 1184 : if (T->strategy != PRST_nextprime) return NULL; /* we're done */
1033 : /* overflow ulong, switch to GEN */
1034 48 : u = ULONG_MAX;
1035 48 : if (T->q > 1) u -= (ULONG_MAX-T->c) % T->q;
1036 48 : affui(u, T->pp);
1037 : }
1038 11847 : av = avma; p = T->pp;
1039 11847 : if (T->q == 1)
1040 : {
1041 11669 : p = nextprime(addiu(p, 1));
1042 11669 : if (T->bb && abscmpii(p, T->bb) > 0) return gc_NULL(av);
1043 : } else do {
1044 3341 : p = T->qq? addii(p, T->qq): addiu(p, T->q);
1045 3341 : if (T->bb && abscmpii(p, T->bb) > 0) return gc_NULL(av);
1046 3285 : } while (!BPSW_psp(p));
1047 11663 : affii(p, T->pp); return gc_const(av, T->pp);
1048 : }
1049 :
1050 : void
1051 1456 : forprimestep(GEN a, GEN b, GEN q, GEN code)
1052 : {
1053 1456 : pari_sp av = avma;
1054 : forprime_t T;
1055 :
1056 1456 : if (!forprimestep_init(&T, a,b,q)) { set_avma(av); return; }
1057 :
1058 1442 : push_lex(T.pp,code);
1059 353929 : while(forprime_next(&T))
1060 : {
1061 352914 : closure_evalvoid(code); if (loop_break()) break;
1062 : /* p changed in 'code', complain */
1063 352494 : if (get_lex(-1) != T.pp)
1064 7 : pari_err(e_MISC,"prime index read-only: was changed to %Ps", get_lex(-1));
1065 : }
1066 1435 : pop_lex(1); set_avma(av);
1067 : }
1068 : void
1069 1330 : forprime(GEN a, GEN b, GEN code) { return forprimestep(a,b,NULL,code); }
1070 :
1071 : int
1072 70 : forcomposite_init(forcomposite_t *C, GEN a, GEN b)
1073 : {
1074 70 : pari_sp av = avma;
1075 70 : a = gceil(a);
1076 70 : if (typ(a)!=t_INT) pari_err_TYPE("forcomposite",a);
1077 70 : if (b) {
1078 63 : if (typ(b) == t_INFINITY) b = NULL;
1079 : else
1080 : {
1081 56 : b = gfloor(b);
1082 56 : if (typ(b)!=t_INT) pari_err_TYPE("forcomposite",b);
1083 : }
1084 : }
1085 70 : if (signe(a) < 0) pari_err_DOMAIN("forcomposite", "a", "<", gen_0, a);
1086 70 : if (abscmpiu(a, 4) < 0) a = utoipos(4);
1087 70 : C->first = 1;
1088 70 : if (!forprime_init(&C->T, a,b) && cmpii(a,b) > 0)
1089 : {
1090 7 : C->n = gen_1; /* in case caller forgets to check the return value */
1091 7 : C->b = gen_0; return gc_bool(av,0);
1092 : }
1093 63 : C->n = setloop(a);
1094 63 : C->b = b;
1095 63 : C->p = NULL; return 1;
1096 : }
1097 :
1098 : GEN
1099 238 : forcomposite_next(forcomposite_t *C)
1100 : {
1101 238 : if (C->first) /* first call ever */
1102 : {
1103 63 : C->first = 0;
1104 63 : C->p = forprime_next(&C->T);
1105 : }
1106 : else
1107 175 : C->n = incloop(C->n);
1108 238 : if (C->p)
1109 : {
1110 161 : if (cmpii(C->n, C->p) < 0) return C->n;
1111 77 : C->n = incloop(C->n);
1112 : /* n = p+1 */
1113 77 : C->p = forprime_next(&C->T); /* nextprime(p) > n */
1114 77 : if (C->p) return C->n;
1115 : }
1116 105 : if (!C->b || cmpii(C->n, C->b) <= 0) return C->n;
1117 42 : return NULL;
1118 : }
1119 :
1120 : void
1121 70 : forcomposite(GEN a, GEN b, GEN code)
1122 : {
1123 70 : pari_sp av = avma;
1124 : forcomposite_t T;
1125 : GEN n;
1126 70 : if (!forcomposite_init(&T,a,b)) return;
1127 63 : push_lex(T.n,code);
1128 238 : while((n = forcomposite_next(&T)))
1129 : {
1130 196 : closure_evalvoid(code); if (loop_break()) break;
1131 : /* n changed in 'code', complain */
1132 182 : if (get_lex(-1) != n)
1133 7 : pari_err(e_MISC,"index read-only: was changed to %Ps", get_lex(-1));
1134 : }
1135 56 : pop_lex(1); set_avma(av);
1136 : }
|
