File Coverage

prime_nth_count.c
Criterion Covered Total %
statement 373 424 87.9
branch 388 548 70.8
condition n/a
subroutine n/a
pod n/a
total 761 972 78.2


line stmt bran cond sub pod time code
1             #include
2             #include
3             #include
4              
5             #define FUNC_popcnt 1
6             #define FUNC_isqrt 1
7             #include "ptypes.h"
8             #include "sieve.h"
9             #include "cache.h"
10             #include "lmo.h"
11             #include "constants.h"
12             #include "prime_nth_count.h"
13             #include "util.h"
14              
15             #include
16             #if _MSC_VER || defined(__IBMC__) || defined(__IBMCPP__) || (defined(__STDC_VERSION__) && __STDC_VERSION >= 199901L)
17             /* math.h should give us these as functions or macros.
18             *
19             * extern long double floorl(long double);
20             * extern long double ceill(long double);
21             * extern long double sqrtl(long double);
22             * extern long double logl(long double);
23             */
24             #else
25             #define floorl(x) (long double) floor( (double) (x) )
26             #define ceill(x) (long double) ceil( (double) (x) )
27             #define sqrtl(x) (long double) sqrt( (double) (x) )
28             #define logl(x) (long double) log( (double) (x) )
29             #endif
30              
31             #if defined(__GNUC__)
32             #define word_unaligned(m,wordsize) ((uintptr_t)m & (wordsize-1))
33             #else /* uintptr_t is part of C99 */
34             #define word_unaligned(m,wordsize) ((unsigned int)m & (wordsize-1))
35             #endif
36              
37             /* TODO: This data is duplicated in util.c. */
38             static const unsigned char prime_sieve30[] =
39             {0x01,0x20,0x10,0x81,0x49,0x24,0xc2,0x06,0x2a,0xb0,0xe1,0x0c,0x15,0x59,0x12,
40             0x61,0x19,0xf3,0x2c,0x2c,0xc4,0x22,0xa6,0x5a,0x95,0x98,0x6d,0x42,0x87,0xe1,
41             0x59,0xa9,0xa9,0x1c,0x52,0xd2,0x21,0xd5,0xb3,0xaa,0x26,0x5c,0x0f,0x60,0xfc,
42             0xab,0x5e,0x07,0xd1,0x02,0xbb,0x16,0x99,0x09,0xec,0xc5,0x47,0xb3,0xd4,0xc5,
43             0xba,0xee,0x40,0xab,0x73,0x3e,0x85,0x4c,0x37,0x43,0x73,0xb0,0xde,0xa7,0x8e,
44             0x8e,0x64,0x3e,0xe8,0x10,0xab,0x69,0xe5,0xf7,0x1a,0x7c,0x73,0xb9,0x8d,0x04,
45             0x51,0x9a,0x6d,0x70,0xa7,0x78,0x2d,0x6d,0x27,0x7e,0x9a,0xd9,0x1c,0x5f,0xee,
46             0xc7,0x38,0xd9,0xc3,0x7e,0x14,0x66,0x72,0xae,0x77,0xc1,0xdb,0x0c,0xcc,0xb2,
47             0xa5,0x74,0xe3,0x58,0xd5,0x4b,0xa7,0xb3,0xb1,0xd9,0x09,0xe6,0x7d,0x23,0x7c,
48             0x3c,0xd3,0x0e,0xc7,0xfd,0x4a,0x32,0x32,0xfd,0x4d,0xb5,0x6b,0xf3,0xa8,0xb3,
49             0x85,0xcf,0xbc,0xf4,0x0e,0x34,0xbb,0x93,0xdb,0x07,0xe6,0xfe,0x6a,0x57,0xa3,
50             0x8c,0x15,0x72,0xdb,0x69,0xd4,0xaf,0x59,0xdd,0xe1,0x3b,0x2e,0xb7,0xf9,0x2b,
51             0xc5,0xd0,0x8b,0x63,0xf8,0x95,0xfa,0x77,0x40,0x97,0xea,0xd1,0x9f,0xaa,0x1c,
52             0x48,0xae,0x67,0xf7,0xeb,0x79,0xa5,0x55,0xba,0xb2,0xb6,0x8f,0xd8,0x2d,0x6c,
53             0x2a,0x35,0x54,0xfd,0x7c,0x9e,0xfa,0xdb,0x31,0x78,0xdd,0x3d,0x56,0x52,0xe7,
54             0x73,0xb2,0x87,0x2e,0x76,0xe9,0x4f,0xa8,0x38,0x9d,0x5d,0x3f,0xcb,0xdb,0xad,
55             0x51,0xa5,0xbf,0xcd,0x72,0xde,0xf7,0xbc,0xcb,0x49,0x2d,0x49,0x26,0xe6,0x1e,
56             0x9f,0x98,0xe5,0xc6,0x9f,0x2f,0xbb,0x85,0x6b,0x65,0xf6,0x77,0x7c,0x57,0x8b,
57             0xaa,0xef,0xd8,0x5e,0xa2,0x97,0xe1,0xdc,0x37,0xcd,0x1f,0xe6,0xfc,0xbb,0x8c,
58             0xb7,0x4e,0xc7,0x3c,0x19,0xd5,0xa8,0x9e,0x67,0x4a,0xe3,0xf5,0x97,0x3a,0x7e,
59             0x70,0x53,0xfd,0xd6,0xe5,0xb8,0x1c,0x6b,0xee,0xb1,0x9b,0xd1,0xeb,0x34,0xc2,
60             0x23,0xeb,0x3a,0xf9,0xef,0x16,0xd6,0x4e,0x7d,0x16,0xcf,0xb8,0x1c,0xcb,0xe6,
61             0x3c,0xda,0xf5,0xcf};
62             #define NPRIME_SIEVE30 (sizeof(prime_sieve30)/sizeof(prime_sieve30[0]))
63              
64             static const unsigned short primes_small[] =
65             {0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
66             101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
67             193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,
68             293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,
69             409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499};
70             #define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))
71              
72              
73             static const unsigned char byte_zeros[256] =
74             {8,7,7,6,7,6,6,5,7,6,6,5,6,5,5,4,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,
75             7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
76             7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
77             6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
78             7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
79             6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
80             6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
81             5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,4,3,3,2,3,2,2,1,3,2,2,1,2,1,1,0};
82 23248           static UV count_zero_bits(const unsigned char* m, UV nbytes)
83             {
84 23248           UV count = 0;
85             #if BITS_PER_WORD == 64
86 23248 100         if (nbytes >= 16) {
87 156791 100         while ( word_unaligned(m,sizeof(UV)) && nbytes--)
    50          
88 135659           count += byte_zeros[*m++];
89 21132 50         if (nbytes >= 8) {
90 21132           UV* wordptr = (UV*)m;
91 21132           UV nwords = nbytes / 8;
92 21132           UV nzeros = nwords * 64;
93 21132           m += nwords * 8;
94 21132           nbytes %= 8;
95 299353 100         while (nwords--)
96 278221           nzeros -= popcnt(*wordptr++);
97 21132           count += nzeros;
98             }
99             }
100             #endif
101 117186 100         while (nbytes--)
102 93938           count += byte_zeros[*m++];
103 23248           return count;
104             }
105              
106             /* Given a sieve of size nbytes, walk it counting zeros (primes) until:
107             *
108             * (1) we counted them all: return the count, which will be less than maxcount.
109             *
110             * (2) we hit maxcount: set position to the index of the maxcount'th prime
111             * and return count (which will be equal to maxcount).
112             */
113 1086           static UV count_segment_maxcount(const unsigned char* sieve, UV base, UV nbytes, UV maxcount, UV* pos)
114             {
115 1086           UV count = 0;
116 1086           UV byte = 0;
117 1086           const unsigned char* sieveptr = sieve;
118 1086           const unsigned char* maxsieve = sieve + nbytes;
119              
120 1086 50         MPUassert(sieve != 0, "count_segment_maxcount incorrect args");
121 1086 50         MPUassert(pos != 0, "count_segment_maxcount incorrect args");
122 1086           *pos = 0;
123 1086 50         if ( (nbytes == 0) || (maxcount == 0) )
    50          
124 0           return 0;
125              
126             /* Do fixed-length word counts to start, with possible overcounting */
127 4499 100         while ((count+64) < maxcount && sieveptr < maxsieve) {
    50          
128 3413           UV top = base + 3*maxcount;
129 3413 100         UV div = (top < 8000) ? 8 : /* 8 cannot overcount */
    100          
    100          
130             (top < 1000000) ? 4 :
131             (top < 10000000) ? 3 : 2;
132 3413           UV minbytes = (maxcount-count)/div;
133 3413 50         if (minbytes > (UV)(maxsieve-sieveptr)) minbytes = maxsieve-sieveptr;
134 3413           count += count_zero_bits(sieveptr, minbytes);
135 3413           sieveptr += minbytes;
136             }
137             /* Count until we reach the end or >= maxcount */
138 15194 50         while ( (sieveptr < maxsieve) && (count < maxcount) )
    100          
139 14108           count += byte_zeros[*sieveptr++];
140             /* If we went too far, back up. */
141 2172 100         while (count >= maxcount)
142 1086           count -= byte_zeros[*--sieveptr];
143             /* We counted this many bytes */
144 1086           byte = sieveptr - sieve;
145              
146 1086 50         MPUassert(count < maxcount, "count_segment_maxcount wrong count");
147              
148 1086 50         if (byte == nbytes)
149 0           return count;
150              
151             /* The result is somewhere in the next byte */
152 16700 50         START_DO_FOR_EACH_SIEVE_PRIME(sieve, 0, byte*30+1, nbytes*30-1)
    50          
    100          
    50          
    50          
153 2774 100         if (++count == maxcount) { *pos = p; return count; }
154 0           END_DO_FOR_EACH_SIEVE_PRIME;
155              
156 0           MPUassert(0, "count_segment_maxcount failure");
157             return 0;
158             }
159              
160             /* Given a sieve of size nbytes, counting zeros (primes) but excluding the
161             * areas outside lowp and highp.
162             */
163 21715           static UV count_segment_ranged(const unsigned char* sieve, UV nbytes, UV lowp, UV highp)
164             {
165             UV count, hi_d, lo_d, lo_m;
166              
167 21715 50         MPUassert( sieve != 0, "count_segment_ranged incorrect args");
168 21715 50         if (nbytes == 0) return 0;
169              
170 21715           count = 0;
171 21715           hi_d = highp/30;
172              
173 21715 50         if (hi_d >= nbytes) {
174 0           hi_d = nbytes-1;
175 0           highp = hi_d*30+29;
176             }
177              
178 21715 50         if (highp < lowp)
179 0           return 0;
180              
181             #if 0
182             /* Dead simple way */
183             START_DO_FOR_EACH_SIEVE_PRIME(sieve, 0, lowp, highp)
184             count++;
185             END_DO_FOR_EACH_SIEVE_PRIME;
186             return count;
187             #endif
188              
189 21715           lo_d = lowp/30;
190 21715           lo_m = lowp - lo_d*30;
191             /* Count first fragment */
192 21715 100         if (lo_m > 1) {
193 21599           UV upper = (highp <= (lo_d*30+29)) ? highp : (lo_d*30+29);
194 191812 50         START_DO_FOR_EACH_SIEVE_PRIME(sieve, 0, lowp, upper)
    100          
    100          
    100          
    100          
195 148565           count++;
196 21599           END_DO_FOR_EACH_SIEVE_PRIME;
197 21599           lowp = upper+2;
198 21599           lo_d = lowp/30;
199             }
200 21715 100         if (highp < lowp)
201 1340           return count;
202              
203             /* Count bytes in the middle */
204             {
205 20375           UV hi_m = highp - hi_d*30;
206 20375           UV count_bytes = hi_d - lo_d + (hi_m == 29);
207 20375 100         if (count_bytes > 0) {
208 19835           count += count_zero_bits(sieve+lo_d, count_bytes);
209 19835           lowp += 30*count_bytes;
210             }
211             }
212 20375 100         if (highp < lowp)
213 1533           return count;
214              
215             /* Count last fragment */
216 343272 50         START_DO_FOR_EACH_SIEVE_PRIME(sieve, 0, lowp, highp)
    100          
    100          
    100          
    100          
217 41206           count++;
218 18842           END_DO_FOR_EACH_SIEVE_PRIME;
219              
220 18842           return count;
221             }
222              
223              
224             /*
225             * The pi(x) prime count functions. prime_count(x) gives an exact number,
226             * but requires determining all the primes up to x, so will be much slower.
227             *
228             * prime_count_lower(x) and prime_count_upper(x) give lower and upper limits,
229             * which will bound the exact value. These bounds should be fairly tight.
230             *
231             * pi_upper(x) - pi(x) pi_lower(x) - pi(x)
232             * < 10 for x < 5_371 < 10 for x < 9_437
233             * < 50 for x < 295_816 < 50 for x < 136_993
234             * < 100 for x < 1_761_655 < 100 for x < 909_911
235             * < 200 for x < 9_987_821 < 200 for x < 8_787_901
236             * < 400 for x < 34_762_891 < 400 for x < 30_332_723
237             * < 1000 for x < 372_748_528 < 1000 for x < 233_000_533
238             * < 5000 for x < 1_882_595_905 < 5000 for x < over 4300M
239             *
240             * The average of the upper and lower bounds is within 9 for all x < 15809, and
241             * within 50 for all x < 1_763_367.
242             *
243             * It is common to use the following Chebyshev inequality for x >= 17:
244             * 1*x/logx <-> 1.25506*x/logx
245             * but this gives terribly loose bounds.
246             *
247             * Rosser and Schoenfeld's bound for x >= 67 of
248             * x/(logx-1/2) <-> x/(logx-3/2)
249             * is much tighter. These bounds can be tightened even more.
250             *
251             * The formulas of Dusart for higher x are better yet. I recommend the paper
252             * by Burde for further information. Dusart's thesis is also a good resource.
253             *
254             * I have tweaked the bounds formulas for small (under 70_000M) numbers so they
255             * are tighter. These bounds are verified via trial. The Dusart bounds
256             * (1.8 and 2.51) are used for larger numbers since those are proven.
257             *
258             */
259              
260             #include "prime_count_tables.h"
261              
262 21768           UV segment_prime_count(UV low, UV high)
263             {
264             const unsigned char* cache_sieve;
265             unsigned char* segment;
266             UV segment_size, low_d, high_d;
267 21768           UV count = 0;
268              
269 21768 100         if ((low <= 2) && (high >= 2)) count++;
    100          
270 21768 100         if ((low <= 3) && (high >= 3)) count++;
    100          
271 21768 100         if ((low <= 5) && (high >= 5)) count++;
    100          
272 21768 100         if (low < 7) low = 7;
273              
274 21768 100         if (low > high) return count;
275              
276 21715 100         if (low == 7 && high <= 30*NPRIME_SIEVE30) {
    100          
277 21434           count += count_segment_ranged(prime_sieve30, NPRIME_SIEVE30, low, high);
278 21434           return count;
279             }
280              
281             /* If we have sparse prime count tables, use them here. These will adjust
282             * 'low' and 'count' appropriately for a value slightly less than ours.
283             * This should leave just a small amount of sieving left. They stop at
284             * some point, e.g. 3000M, so we'll get the answer to that point then have
285             * to sieve all the rest. We should be using LMO or Lehmer much earlier. */
286             #ifdef APPLY_TABLES
287 11129 100         APPLY_TABLES
    100          
    100          
    100          
    50          
    100          
    100          
    50          
    100          
    100          
    50          
    100          
    100          
    50          
    50          
    100          
    50          
    0          
    0          
    0          
    0          
    0          
288             #endif
289              
290 281           low_d = low/30;
291 281           high_d = high/30;
292              
293             /* Count full bytes only -- no fragments from primary cache */
294 281           segment_size = get_prime_cache(0, &cache_sieve) / 30;
295 281 100         if (segment_size < high_d) {
296             /* Expand sieve to sqrt(n) */
297 46 50         UV endp = (high_d >= (UV_MAX/30)) ? UV_MAX-2 : 30*high_d+29;
298             release_prime_cache(cache_sieve);
299 46           segment_size = get_prime_cache( isqrt(endp) + 1 , &cache_sieve) / 30;
300             }
301              
302 281 50         if ( (segment_size > 0) && (low_d <= segment_size) ) {
    100          
303             /* Count all the primes in the primary cache in our range */
304 235           count += count_segment_ranged(cache_sieve, segment_size, low, high);
305              
306 235 50         if (high_d < segment_size) {
307             release_prime_cache(cache_sieve);
308 235           return count;
309             }
310              
311 0           low_d = segment_size;
312 0 0         if (30*low_d > low) low = 30*low_d;
313             }
314             release_prime_cache(cache_sieve);
315              
316             /* More primes needed. Repeatedly segment sieve. */
317             {
318 46           void* ctx = start_segment_primes(low, high, &segment);
319             UV seg_base, seg_low, seg_high;
320 92 100         while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
321 46           segment_size = seg_high/30 - seg_low/30 + 1;
322 46           count += count_segment_ranged(segment, segment_size, seg_low-seg_base, seg_high-seg_base);
323             }
324 46           end_segment_primes(ctx);
325             }
326              
327 21768           return count;
328             }
329              
330 683           UV prime_count(UV lo, UV hi)
331             {
332             UV count, range_threshold;
333              
334 683 50         if (lo > hi || hi < 2)
    100          
335 1           return 0;
336              
337             /* We use table acceleration so this is preferable for small inputs */
338 682 100         if (hi < 66000000) return segment_prime_count(lo, hi);
339              
340             /* Rough empirical threshold for when segment faster than LMO */
341 29           range_threshold = hi / (isqrt(hi)/200);
342 29 100         if ( (hi-lo+1) < range_threshold )
343 10           return segment_prime_count(lo, hi);
344              
345 19 50         return LMO_prime_count(hi) - ((lo < 2) ? 0 : LMO_prime_count(lo-1));
346              
347             return count;
348             }
349              
350 36           UV prime_count_approx(UV n)
351             {
352 36 100         if (n < 3000000) return segment_prime_count(2, n);
353 25           return (UV) (RiemannR( (long double) n ) + 0.5 );
354             }
355              
356             /* See http://numbers.computation.free.fr/Constants/Primes/twin.pdf, page 5 */
357             /* Upper limit is in Wu, Acta Arith 114 (2004). 4.48857*x/(log(x)*log(x) */
358 623           UV twin_prime_count_approx(UV n)
359             {
360             /* Best would be another estimate for n < ~ 5000 */
361 623 100         if (n < 2000) return twin_prime_count(3,n);
362             {
363             /* Sebah and Gourdon 2002 */
364 216           const long double two_C2 = 1.32032363169373914785562422L;
365 216           const long double two_over_log_two = 2.8853900817779268147198494L;
366 216           long double ln = (long double) n;
367 216           long double logn = logl(ln);
368 216           long double li2 = Ei(logn) + two_over_log_two-ln/logn;
369             /* try to minimize MSE */
370 216 100         if (n < 32000000) {
371             long double fm;
372 86 50         if (n < 4000) fm = 0.2952;
373 86 100         else if (n < 8000) fm = 0.3152;
374 85 50         else if (n < 16000) fm = 0.3090;
375 85 100         else if (n < 32000) fm = 0.3096;
376 70 100         else if (n < 64000) fm = 0.3100;
377 56 100         else if (n < 128000) fm = 0.3089;
378 41 50         else if (n < 256000) fm = 0.3099;
379 41 100         else if (n < 600000) fm = .3091 + (n-256000) * (.3056-.3091) / (600000-256000);
380 22 50         else if (n < 1000000) fm = .3062 + (n-600000) * (.3042-.3062) / (1000000-600000);
381 22 50         else if (n < 4000000) fm = .3067 + (n-1000000) * (.3041-.3067) / (4000000-1000000);
382 22 50         else if (n <16000000) fm = .3033 + (n-4000000) * (.2983-.3033) / (16000000-4000000);
383 0           else fm = .2980 + (n-16000000) * (.2965-.2980) / (32000000-16000000);
384 86           li2 *= fm * logl(12+logn);
385             }
386 216           return (UV) (two_C2 * li2 + 0.5L);
387             }
388             }
389              
390 16449           UV prime_count_lower(UV n)
391             {
392             long double fn, fl1, fl2, lower, a;
393              
394 16449 100         if (n < 33000) return segment_prime_count(2, n);
395              
396 979           fn = (long double) n;
397 979           fl1 = logl(n);
398 979           fl2 = fl1 * fl1;
399              
400 979 100         if (n <= 300000) { /* Quite accurate and avoids calling Li for speed. */
401 138 100         a = (n < 70200) ? 947 : (n < 176000) ? 904 : 829;
    100          
402 138           lower = fn / (fl1 - 1 - 1/fl1 - 2.85/fl2 - 13.15/(fl1*fl2) + a/(fl2*fl2));
403 841 100         } else if (n < UVCONST(4000000000)) {
404             /* Loose enough that FP differences in Li(n) should be ok. */
405 822           a = (n < 88783) ? 4.0L
406 1644 50         : (n < 300000) ? -3.0L
407 1644 50         : (n < 303000) ? 5.0L
408 1644 50         : (n < 1100000) ? -7.0L
409 1369 100         : (n < 4500000) ? -37.0L
410 581 100         : (n < 10200000) ? -70.0L
411 42 100         : (n < 36900000) ? -53.0L
412 15 100         : (n < 38100000) ? -29.0L
413 7 50         : -84.0L;
414 822           lower = Li(fn) - (sqrtl(fn)/fl1) * (1.94L + 2.50L/fl1 + a/fl2);
415 19 100         } else if (fn < 1e19) { /* Büthe 2015 1.9 1511.02032v1.pdf */
416 17           lower = Li(fn) - (sqrtl(fn)/fl1) * (1.94L + 3.88L/fl1 + 27.57L/fl2);
417             } else { /* Büthe 2014 v3 7.2 1410.7015v3.pdf */
418 2           lower = Li(fn) - fl1*sqrtl(fn)/25.132741228718345907701147L;
419             }
420 979           return (UV) ceill(lower);
421             }
422              
423             typedef struct {
424             UV thresh;
425             float aval;
426             } thresh_t;
427              
428             static const thresh_t _upper_thresh[] = {
429             { 59000, 2.48 },
430             { 355991, 2.54 },
431             { 3550000, 2.51 },
432             { 3560000, 2.49 },
433             { 5000000, 2.48 },
434             { 8000000, 2.47 },
435             { 13000000, 2.46 },
436             { 18000000, 2.45 },
437             { 31000000, 2.44 },
438             { 41000000, 2.43 },
439             { 48000000, 2.42 },
440             { 119000000, 2.41 },
441             { 182000000, 2.40 },
442             { 192000000, 2.395 },
443             { 213000000, 2.390 },
444             { 271000000, 2.385 },
445             { 322000000, 2.380 },
446             { 400000000, 2.375 },
447             { 510000000, 2.370 },
448             { 682000000, 2.367 },
449             { UVCONST(2953652287), 2.362 }
450             };
451             #define NUPPER_THRESH (sizeof(_upper_thresh)/sizeof(_upper_thresh[0]))
452              
453 6144           UV prime_count_upper(UV n)
454             {
455             int i;
456             long double fn, fl1, fl2, upper, a;
457              
458 6144 100         if (n < 33000) return segment_prime_count(2, n);
459              
460 600           fn = (long double) n;
461 600           fl1 = logl(n);
462 600           fl2 = fl1 * fl1;
463              
464 600 100         if (BITS_PER_WORD == 32 || fn <= 821800000.0) { /* Dusart 2010, page 2 */
465 1866 50         for (i = 0; i < (int)NUPPER_THRESH; i++)
466 1866 100         if (n < _upper_thresh[i].thresh)
467 579           break;
468 579 50         a = (i < (int)NUPPER_THRESH) ? _upper_thresh[i].aval : 2.334L;
469 579           upper = fn/fl1 * (1.0L + 1.0L/fl1 + a/fl2);
470 21 100         } else if (fn < 1e19) { /* Büthe 2015 1.10 Skewes number lower limit */
471 19 100         a = (fn < 1100000000.0) ? 0.032 /* Empirical */
    100          
    100          
472             : (fn < 10010000000.0) ? 0.027 /* Empirical */
473             : (fn < 101260000000.0) ? 0.021 /* Empirical */
474             : 0.0;
475 19           upper = Li(fn) - a * fl1*sqrtl(fn)/25.132741228718345907701147L;
476             } else { /* Büthe 2014 7.4 */
477 2           upper = Li(fn) + fl1*sqrtl(fn)/25.132741228718345907701147L;
478             }
479 600           return (UV) floorl(upper);
480             }
481              
482 3004           static void simple_nth_limits(UV *lo, UV *hi, long double n, long double logn, long double loglogn) {
483 3004 100         const long double a = (n < 228) ? .6483 : (n < 948) ? .8032 : (n < 2195) ? .8800 : (n < 39017) ? .9019 : .9484;
    100          
    100          
    100          
484 3004           *lo = n * (logn + loglogn - 1.0 + ((loglogn-2.10)/logn));
485 3004           *hi = n * (logn + loglogn - a);
486 3004 50         if (*hi < *lo) *hi = MPU_MAX_PRIME;
487 3004           }
488              
489             /* The nth prime will be less or equal to this number */
490 2927           UV nth_prime_upper(UV n)
491             {
492             long double fn, flogn, flog2n, upper;
493              
494 2927 100         if (n < NPRIMES_SMALL)
495 429           return primes_small[n];
496              
497 2498           fn = (long double) n;
498 2498           flogn = logl(n);
499 2498           flog2n = logl(flogn); /* Note distinction between log_2(n) and log^2(n) */
500              
501 2498 100         if (n < 688383) {
502             UV lo,hi;
503 2406           simple_nth_limits(&lo, &hi, fn, flogn, flog2n);
504 18352 100         while (lo < hi) {
505 15946           UV mid = lo + (hi-lo)/2;
506 15946 100         if (prime_count_lower(mid) < n) lo = mid+1;
507 8131           else hi = mid;
508             }
509 2406           return lo;
510             }
511              
512             /* Dusart 2010 page 2 */
513 92           upper = fn * (flogn + flog2n - 1.0 + ((flog2n-2.00)/flogn));
514 92 100         if (n >= 46254381) {
515             /* Axler 2017 http//arxiv.org/pdf/1706.03651.pdf Corollary 1.2 */
516 12           upper -= fn * ((flog2n*flog2n-6*flog2n+10.667)/(2*flogn*flogn));
517 80 100         } else if (n >= 8009824) {
518             /* Axler 2013 page viii Korollar G */
519 61           upper -= fn * ((flog2n*flog2n-6*flog2n+10.273)/(2*flogn*flogn));
520             }
521              
522 92 50         if (upper >= (long double)UV_MAX) {
523 0 0         if (n <= MPU_MAX_PRIME_IDX) return MPU_MAX_PRIME;
524 0           croak("nth_prime_upper(%"UVuf") overflow", n);
525             }
526              
527 92           return (UV) floorl(upper);
528             }
529              
530             /* The nth prime will be greater than or equal to this number */
531 1313           UV nth_prime_lower(UV n)
532             {
533             double fn, flogn, flog2n, lower;
534              
535 1313 100         if (n < NPRIMES_SMALL)
536 638           return primes_small[n];
537              
538 675           fn = (double) n;
539 675           flogn = log(n);
540 675           flog2n = log(flogn);
541              
542             /* For small values, do a binary search on the inverse prime count */
543 675 100         if (n < 2000000) {
544             UV lo,hi;
545 598           simple_nth_limits(&lo, &hi, fn, flogn, flog2n);
546 4392 100         while (lo < hi) {
547 3794           UV mid = lo + (hi-lo)/2;
548 3794 100         if (prime_count_upper(mid) < n) lo = mid+1;
549 1870           else hi = mid;
550             }
551 598           return lo;
552             }
553              
554             { /* Axler 2017 http//arxiv.org/pdf/1706.03651.pdf Corollary 1.4 */
555 77 100         double b1 = (n < 56000000) ? 11.200 : 11.508;
556 77           lower = fn * (flogn + flog2n-1.0 + ((flog2n-2.00)/flogn) - ((flog2n*flog2n-6*flog2n+b1)/(2*flogn*flogn)));
557             }
558              
559 77           return (UV) ceill(lower);
560             }
561              
562 25           UV nth_prime_approx(UV n)
563             {
564 25 100         return (n < NPRIMES_SMALL) ? primes_small[n] : inverse_R(n);
565             }
566              
567              
568 7006           UV nth_prime(UV n)
569             {
570             const unsigned char* cache_sieve;
571             unsigned char* segment;
572             UV upper_limit, segbase, segment_size, p, count, target;
573              
574             /* If very small, return the table entry */
575 7006 100         if (n < NPRIMES_SMALL)
576 5917           return primes_small[n];
577              
578             /* Determine a bound on the nth prime. We know it comes before this. */
579 1089           upper_limit = nth_prime_upper(n);
580 1089 50         MPUassert(upper_limit > 0, "nth_prime got an upper limit of 0");
581 1089           p = count = 0;
582 1089           target = n-3;
583              
584             /* For relatively small values, generate a sieve and count the results.
585             *
586             * For larger values, compute an approximate low estimate, use our fast
587             * prime count, then segment sieve forwards or backwards for the rest.
588             */
589 1089 100         if (upper_limit <= get_prime_cache(0, 0) || upper_limit <= 32*1024*30) {
    100          
590             /* Generate a sieve and count. */
591 1067           segment_size = get_prime_cache(upper_limit, &cache_sieve) / 30;
592             /* Count up everything in the cached sieve. */
593 2134 50         if (segment_size > 0)
594 1067           count += count_segment_maxcount(cache_sieve, 0, segment_size, target, &p);
595             release_prime_cache(cache_sieve);
596             } else {
597             /* A binary search on RiemannR is nice, but ends up either often being
598             * being higher (requiring going backwards) or biased and then far too
599             * low. Using the inverse Li is easier and more consistent. */
600 22           UV lower_limit = inverse_li(n);
601             /* For even better performance, add in half the usual correction, which
602             * will get us even closer, so even less sieving required. However, it
603             * is now possible to get a result higher than the value, so we'll need
604             * to handle that case. It still ends up being a better deal than R,
605             * given that we don't have a fast backward sieve. */
606 22           lower_limit += inverse_li(isqrt(n))/4;
607 22           segment_size = lower_limit / 30;
608 22           lower_limit = 30 * segment_size - 1;
609 22           count = prime_count(2,lower_limit);
610              
611             /* printf("We've estimated %lu too %s.\n", (count>n)?count-n:n-count, (count>n)?"FAR":"little"); */
612             /* printf("Our limit %lu %s a prime\n", lower_limit, is_prime(lower_limit) ? "is" : "is not"); */
613              
614 22 100         if (count >= n) { /* Too far. Walk backwards */
615 3 50         if (is_prime(lower_limit)) count--;
616 7 100         for (p = 0; p <= (count-n); p++)
617 4           lower_limit = prev_prime(lower_limit);
618 3           return lower_limit;
619             }
620 19           count -= 3;
621              
622             /* Make sure the segment siever won't have to keep resieving. */
623 19           prime_precalc(isqrt(upper_limit));
624             }
625              
626 1086 100         if (count == target)
627 1067           return p;
628              
629             /* Start segment sieving. Get memory to sieve into. */
630 19           segbase = segment_size;
631 19           segment = get_prime_segment(&segment_size);
632              
633 38 100         while (count < target) {
634             /* Limit the segment size if we know the answer comes earlier */
635 19 100         if ( (30*(segbase+segment_size)+29) > upper_limit )
636 18           segment_size = (upper_limit - segbase*30 + 30) / 30;
637              
638             /* Do the actual sieving in the range */
639 19           sieve_segment(segment, segbase, segbase + segment_size-1);
640              
641             /* Count up everything in this segment */
642 19           count += count_segment_maxcount(segment, 30*segbase, segment_size, target-count, &p);
643              
644 19 50         if (count < target)
645 0           segbase += segment_size;
646             }
647 19           release_prime_segment(segment);
648 19 50         MPUassert(count == target, "nth_prime got incorrect count");
649 7006           return ( (segbase*30) + p );
650             }
651              
652             /******************************************************************************/
653             /* TWIN PRIMES */
654             /******************************************************************************/
655              
656             #if BITS_PER_WORD < 64
657             static const UV twin_steps[] =
658             {58980,48427,45485,43861,42348,41457,40908,39984,39640,39222,
659             373059,353109,341253,332437,326131,320567,315883,312511,309244,
660             2963535,2822103,2734294,2673728,
661             };
662             static const unsigned int twin_num_exponents = 3;
663             static const unsigned int twin_last_mult = 4; /* 4000M */
664             #else
665             static const UV twin_steps[] =
666             {58980,48427,45485,43861,42348,41457,40908,39984,39640,39222,
667             373059,353109,341253,332437,326131,320567,315883,312511,309244,
668             2963535,2822103,2734294,2673728,2626243,2585752,2554015,2527034,2501469,
669             24096420,23046519,22401089,21946975,21590715,21300632,21060884,20854501,20665634,
670             199708605,191801047,186932018,183404596,180694619,178477447,176604059,174989299,173597482,
671             1682185723,1620989842,1583071291,1555660927,1534349481,1517031854,1502382532,1489745250, 1478662752,
672             14364197903,13879821868,13578563641,13361034187,13191416949,13053013447,12936030624,12835090276, 12746487898,
673             124078078589,120182602778,117753842540,115995331742,114622738809,113499818125,112551549250,111732637241,111012321565,
674             1082549061370,1050759497170,1030883829367,1016473645857,1005206830409,995980796683,988183329733,981441437376,975508027029,
675             9527651328494, 9264843314051, 9100153493509, 8980561036751, 8886953365929, 8810223086411, 8745329823109, 8689179566509, 8639748641098,
676             84499489470819, 82302056642520, 80922166953330, 79918799449753, 79132610984280, 78487688897426, 77941865286827, 77469296499217, 77053075040105,
677             754527610498466, 735967887462370, 724291736697048,
678             };
679             static const unsigned int twin_num_exponents = 12;
680             static const unsigned int twin_last_mult = 4; /* 4e19 */
681             #endif
682              
683 418           UV twin_prime_count(UV beg, UV end)
684             {
685             unsigned char* segment;
686 418           UV sum = 0;
687              
688             /* First use the tables of #e# from 1e7 to 2e16. */
689 418 100         if (beg <= 3 && end >= 10000000) {
    100          
690 6           UV mult, exp, step = 0, base = 10000000;
691 40 50         for (exp = 0; exp < twin_num_exponents && end >= base; exp++) {
    100          
692 312 100         for (mult = 1; mult < 10 && end >= mult*base; mult++) {
    100          
693 278           sum += twin_steps[step++];
694 278           beg = mult*base;
695 278 50         if (exp == twin_num_exponents-1 && mult >= twin_last_mult) break;
    0          
696             }
697 34           base *= 10;
698             }
699             }
700 418 100         if (beg <= 3 && end >= 3) sum++;
    50          
701 418 100         if (beg <= 5 && end >= 5) sum++;
    50          
702 418 100         if (beg < 11) beg = 7;
703 418 50         if (beg <= end) {
704             /* Make end points odd */
705 418           beg |= 1;
706 418           end = (end-1) | 1;
707             /* Cheesy way of counting the partial-byte edges */
708 5392 100         while ((beg % 30) != 1) {
709 4974 100         if (is_prime(beg) && is_prime(beg+2) && beg <= end) sum++;
    100          
    50          
710 4974           beg += 2;
711             }
712 3264 100         while ((end % 30) != 29) {
713 2861 100         if (is_prime(end) && is_prime(end+2) && beg <= end) sum++;
    100          
    50          
714 2861 100         end -= 2; if (beg > end) break;
715             }
716             }
717 418 100         if (beg <= end) {
718             UV seg_base, seg_low, seg_high;
719 403           void* ctx = start_segment_primes(beg, end, &segment);
720 806 100         while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
721 403           UV bytes = seg_high/30 - seg_low/30 + 1;
722             unsigned char s;
723 403           const unsigned char* sp = segment;
724 403           const unsigned char* const spend = segment + bytes - 1;
725 53716 100         while (sp < spend) {
726 53313           s = *sp++;
727 53313 100         if (!(s & 0x0C)) sum++;
728 53313 100         if (!(s & 0x30)) sum++;
729 53313 100         if (!(s & 0x80) && !(*sp & 0x01)) sum++;
    100          
730             }
731 403           s = *sp;
732 403 100         if (!(s & 0x0C)) sum++;
733 403 100         if (!(s & 0x30)) sum++;
734 403 100         if (!(s & 0x80) && is_prime(seg_high+2)) sum++;
    100          
735             }
736 403           end_segment_primes(ctx);
737             }
738 418           return sum;
739             }
740              
741 54           UV nth_twin_prime(UV n)
742             {
743             unsigned char* segment;
744             double dend;
745 54           UV nth = 0;
746             UV beg, end;
747              
748 54 100         if (n < 6) {
749 6           switch (n) {
750 0           case 0: nth = 0; break;
751 1           case 1: nth = 3; break;
752 1           case 2: nth = 5; break;
753 1           case 3: nth =11; break;
754 1           case 4: nth =17; break;
755             case 5:
756 2           default: nth =29; break;
757             }
758 6           return nth;
759             }
760              
761 48           end = UV_MAX - 16;
762 48           dend = 800.0 + 1.01L * (double)nth_twin_prime_approx(n);
763 48 50         if (dend < (double)end) end = (UV) dend;
764              
765 48           beg = 2;
766 48 50         if (n > 58980) { /* Use twin_prime_count tables to accelerate if possible */
767 0           UV mult, exp, step = 0, base = 10000000;
768 0 0         for (exp = 0; exp < twin_num_exponents && end >= base; exp++) {
    0          
769 0 0         for (mult = 1; mult < 10 && n > twin_steps[step]; mult++) {
    0          
770 0           n -= twin_steps[step++];
771 0           beg = mult*base;
772 0 0         if (exp == twin_num_exponents-1 && mult >= twin_last_mult) break;
    0          
773             }
774 0           base *= 10;
775             }
776             }
777 48 50         if (beg == 2) { beg = 31; n -= 5; }
778              
779             {
780             UV seg_base, seg_low, seg_high;
781 48           void* ctx = start_segment_primes(beg, end, &segment);
782 96 100         while (n && next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
    50          
783 48           UV p, bytes = seg_high/30 - seg_low/30 + 1;
784 48           UV s = ((UV)segment[0]) << 8;
785 4428 50         for (p = 0; p < bytes; p++) {
786 4428           s >>= 8;
787 4428 50         if (p+1 < bytes) s |= (((UV)segment[p+1]) << 8);
788 0 0         else if (!is_prime(seg_high+2)) s |= 0xFF00;
789 4428 100         if (!(s & 0x000C) && !--n) { nth=seg_base+p*30+11; break; }
    100          
790 4409 100         if (!(s & 0x0030) && !--n) { nth=seg_base+p*30+17; break; }
    100          
791 4396 100         if (!(s & 0x0180) && !--n) { nth=seg_base+p*30+29; break; }
    100          
792             }
793             }
794 48           end_segment_primes(ctx);
795             }
796 54           return nth;
797             }
798              
799 58           UV nth_twin_prime_approx(UV n)
800             {
801 58           long double fn = (long double) n;
802 58           long double flogn = logl(n);
803 58           long double fnlog2n = fn * flogn * flogn;
804             UV lo, hi;
805              
806 58 100         if (n < 6)
807 1           return nth_twin_prime(n);
808              
809             /* Binary search on the TPC estimate.
810             * Good results require that the TPC estimate is both fast and accurate.
811             * These bounds are good for the actual nth_twin_prime values.
812             */
813 57           lo = (UV) (0.9 * fnlog2n);
814 57 50         hi = (UV) ( (n >= 1e16) ? (1.04 * fnlog2n) :
    50          
    100          
    100          
815 57           (n >= 1e13) ? (1.10 * fnlog2n) :
816 59           (n >= 1e7 ) ? (1.31 * fnlog2n) :
817 5           (n >= 1200) ? (1.70 * fnlog2n) :
818 50           (2.3 * fnlog2n + 5) );
819 57 50         if (hi <= lo) hi = UV_MAX;
820 673 100         while (lo < hi) {
821 616           UV mid = lo + (hi-lo)/2;
822 616 100         if (twin_prime_count_approx(mid) < fn) lo = mid+1;
823 287           else hi = mid;
824             }
825 57           return lo;
826             }
827              
828             /******************************************************************************/
829             /* SUMS */
830             /******************************************************************************/
831              
832             /* The fastest way to compute the sum of primes is using a combinatorial
833             * algorithm such as Deleglise 2012. Since this code is purely native,
834             * it will overflow a 64-bit result quite quickly. Hence a relatively small
835             * table plus sum over sieved primes works quite well.
836             *
837             * The following info is useful if we ever return 128-bit results or for a
838             * GMP implementation.
839             *
840             * Combinatorial sum of primes < n. Call with phisum(n, isqrt(n)).
841             * Needs optimization, either caching, Lehmer, or LMO.
842             * http://mathoverflow.net/questions/81443/fastest-algorithm-to-compute-the-sum-of-primes
843             * http://www.ams.org/journals/mcom/2009-78-268/S0025-5718-09-02249-2/S0025-5718-09-02249-2.pdf
844             * http://mathematica.stackexchange.com/questions/80291/efficient-way-to-sum-all-the-primes-below-n-million-in-mathematica
845             * Deleglise 2012, page 27, simple Meissel:
846             * y = x^1/3
847             * a = Pi(y)
848             * Pi_f(x) = phisum(x,a) + Pi_f(y) - 1 - P_2(x,a)
849             * P_2(x,a) = sum prime p : y < p <= sqrt(x) of f(p) * Pi_f(x/p) -
850             * sum prime p : y < p <= sqrt(x) of f(p) * Pi_f(p-1)
851             */
852              
853             static const unsigned char byte_sum[256] =
854             {120,119,113,112,109,108,102,101,107,106,100,99,96,95,89,88,103,102,96,95,92,
855             91,85,84,90,89,83,82,79,78,72,71,101,100,94,93,90,89,83,82,88,87,81,80,77,
856             76,70,69,84,83,77,76,73,72,66,65,71,70,64,63,60,59,53,52,97,96,90,89,86,85,
857             79,78,84,83,77,76,73,72,66,65,80,79,73,72,69,68,62,61,67,66,60,59,56,55,49,
858             48,78,77,71,70,67,66,60,59,65,64,58,57,54,53,47,46,61,60,54,53,50,49,43,42,
859             48,47,41,40,37,36,30,29,91,90,84,83,80,79,73,72,78,77,71,70,67,66,60,59,74,
860             73,67,66,63,62,56,55,61,60,54,53,50,49,43,42,72,71,65,64,61,60,54,53,59,58,
861             52,51,48,47,41,40,55,54,48,47,44,43,37,36,42,41,35,34,31,30,24,23,68,67,61,
862             60,57,56,50,49,55,54,48,47,44,43,37,36,51,50,44,43,40,39,33,32,38,37,31,30,
863             27,26,20,19,49,48,42,41,38,37,31,30,36,35,29,28,25,24,18,17,32,31,25,24,21,
864             20,14,13,19,18,12,11,8,7,1,0};
865              
866             #if BITS_PER_WORD == 64
867             /* We have a much more limited range, so use a fixed interval. We should be
868             * able to get any 64-bit sum in under a half-second. */
869             static const UV sum_table_2e8[] =
870             {1075207199997324,3071230303170813,4990865886639877,6872723092050268,8729485610396243,10566436676784677,12388862798895708,14198556341669206,15997206121881531,17783028661796383,19566685687136351,21339485298848693,23108856419719148,
871             24861364231151903,26619321031799321,28368484289421890,30110050320271201,31856321671656548,33592089385327108,35316546074029522,37040262208390735,38774260466286299,40490125006181147,42207686658844380,43915802985817228,45635106002281013,
872             47337822860157465,49047713696453759,50750666660265584,52449748364487290,54152689180758005,55832433395290183,57540651847418233,59224867245128289,60907462954737468,62597192477315868,64283665223856098,65961576139329367,67641982565760928,
873             69339211720915217,71006044680007261,72690896543747616,74358564592509127,76016548794894677,77694517638354266,79351385193517953,81053240048141953,82698120948724835,84380724263091726,86028655116421543,87679091888973563,89348007111430334,
874             90995902774878695,92678527127292212,94313220293410120,95988730932107432,97603162494502485,99310622699836698,100935243057337310,102572075478649557,104236362884241550,105885045921116836,107546170993472638,109163445284201278,
875             110835950755374921,112461991135144669,114116351921245042,115740770232532531,117408250788520189,119007914428335965,120652479429703269,122317415246500401,123951466213858688,125596789655927842,127204379051939418,128867944265073217,
876             130480037123800711,132121840147764197,133752985360747726,135365954823762234,137014594650995101,138614165689305879,140269121741383097,141915099618762647,143529289083557618,145146413750649432,146751434858695468,148397902396643807,
877             149990139346918801,151661665434334577,153236861034424304,154885985064643097,156500983286383741,158120868946747299,159735201435796748,161399264792716319,162999489977602579,164566400448130092,166219688860475191,167836981098849796,
878             169447127305804401,171078187147848898,172678849082290997,174284436375728242,175918609754056455,177525046501311788,179125593738290153,180765176633753371,182338473848291683,183966529541155489,185585792988238475,187131988176321434,
879             188797837140841381,190397649440649965,191981841583560122,193609739194967419,195166830650558070,196865965063113041,198400070713177440,200057161591648721,201621899486413406,203238279253414934,204790684829891896,206407676204061001,
880             208061050481364659,209641606658938873,211192088300183855,212855420483750498,214394145510853736,216036806225784861,217628995137940563,219277567478725189,220833877268454872,222430818525363309,224007307616922530,225640739533952807,
881             227213096159236934,228853318075566255,230401824696558125,231961445347821085,233593317860593895,235124654760954338,236777716068869769,238431514923528303,239965003913481640,241515977959535845,243129874530821395};
882             #define N_SUM_TABLE (sizeof(sum_table_2e8)/sizeof(sum_table_2e8[0]))
883             #endif
884              
885             /* Add n to the double-word hi,lo */
886             #define ADD_128(hi, lo, n) \
887             do { UV _n = n; \
888             if (_n > (UV_MAX-lo)) { hi++; if (hi == 0) overflow = 1; } \
889             lo += _n; } while (0)
890             #define SET_128(hi, lo, n) \
891             do { hi = (UV) (((n) >> 64) & UV_MAX); \
892             lo = (UV) (((n) ) & UV_MAX); } while (0)
893              
894             /* Legendre method for prime sum */
895 0           int sum_primes128(UV n, UV *hi_sum, UV *lo_sum) {
896             #if BITS_PER_WORD == 64 && HAVE_UINT128
897             uint128_t *V, *S;
898 0           UV j, k, r = isqrt(n), r2 = r + n/(r+1);
899              
900 0 0         New(0, V, r2+1, uint128_t);
901 0 0         New(0, S, r2+1, uint128_t);
902 0 0         for (k = 0; k <= r2; k++) {
903 0 0         uint128_t v = (k <= r) ? k : n/(r2-k+1);
904 0           V[k] = v;
905 0           S[k] = (v*(v+1))/2 - 1;
906             }
907              
908 0 0         START_DO_FOR_EACH_PRIME(2, r) {
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
909 0           uint128_t a, b, sp = S[p-1], p2 = ((uint128_t)p) * p;
910 0 0         for (j = k-1; j > 1 && V[j] >= p2; j--) {
    0          
911 0           a = V[j], b = a/p;
912 0 0         if (a > r) a = r2 - n/a + 1;
913 0 0         if (b > r) b = r2 - n/b + 1;
914 0           S[a] -= p * (S[b] - sp); /* sp = sum of primes less than p */
915             }
916 0           } END_DO_FOR_EACH_PRIME;
917 0           SET_128(*hi_sum, *lo_sum, S[r2]);
918 0           Safefree(V);
919 0           Safefree(S);
920 0           return 1;
921             #else
922             return 0;
923             #endif
924             }
925              
926 1003           int sum_primes(UV low, UV high, UV *return_sum) {
927 1003           UV sum = 0;
928 1003           int overflow = 0;
929              
930 1003 100         if ((low <= 2) && (high >= 2)) sum += 2;
    50          
931 1003 100         if ((low <= 3) && (high >= 3)) sum += 3;
    100          
932 1003 100         if ((low <= 5) && (high >= 5)) sum += 5;
    100          
933 1003 100         if (low < 7) low = 7;
934              
935             /* If we know the range will overflow, return now */
936             #if BITS_PER_WORD == 64
937 1003 100         if (low == 7 && high >= 29505444491) return 0;
    50          
938 1003 50         if (low >= 1e10 && (high-low) >= 32e9) return 0;
    0          
939 1003 50         if (low >= 1e13 && (high-low) >= 5e7) return 0;
    0          
940             #else
941             if (low == 7 && high >= 323381) return 0;
942             #endif
943              
944             #if 1 && BITS_PER_WORD == 64 /* Tables */
945 1003 100         if (low == 7 && high >= 2e8) {
    50          
946             UV step;
947 0 0         for (step = 1; high >= (step * 2e8) && step < N_SUM_TABLE; step++) {
    0          
948 0           sum += sum_table_2e8[step-1];
949 0           low = step * 2e8;
950             }
951             }
952             #endif
953              
954 1003 100         if (low <= high) {
955             unsigned char* segment;
956             UV seg_base, seg_low, seg_high;
957 998           void* ctx = start_segment_primes(low, high, &segment);
958 1996 50         while (!overflow && next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
    100          
959 998           UV bytes = seg_high/30 - seg_low/30 + 1;
960             unsigned char s;
961 998           unsigned char* sp = segment;
962 998           unsigned char* const spend = segment + bytes - 1;
963 998           UV i, p, pbase = 30*(seg_low/30);
964              
965             /* Clear primes before and after our range */
966 998           p = pbase;
967 2005 50         for (i = 0; i < 8 && p+wheel30[i] < low; i++)
    100          
968 1007 100         if ( (*sp & (1<
969 3           *sp |= (1 << i);
970              
971 998           p = 30*(seg_high/30);
972 8982 100         for (i = 0; i < 8; i++)
973 7984 100         if ( (*spend & (1< high )
    100          
974 2459           *spend |= (1 << i);
975              
976 29639 100         while (sp <= spend) {
977 28641           s = *sp++;
978 28641 50         if (sum < (UV_MAX >> 3) && pbase < (UV_MAX >> 5)) {
    50          
979             /* sum block of 8 all at once */
980 28641           sum += pbase * byte_zeros[s] + byte_sum[s];
981             } else {
982             /* sum block of 8, checking for overflow at each step */
983 0 0         for (i = 0; i < byte_zeros[s]; i++) {
984 0 0         if (sum+pbase < sum) overflow = 1;
985 0           sum += pbase;
986             }
987 0 0         if (sum+byte_sum[s] < sum) overflow = 1;
988 0           sum += byte_sum[s];
989 0 0         if (overflow) break;
990             }
991 28641           pbase += 30;
992             }
993             }
994 998           end_segment_primes(ctx);
995             }
996 1003 50         if (!overflow && return_sum != 0) *return_sum = sum;
    50          
997 1003           return !overflow;
998             }
999              
1000 526           double ramanujan_sa_gn(UV un)
1001             {
1002 526           long double n = (long double) un;
1003 526           long double logn = logl(n);
1004 526           long double log2 = logl(2);
1005              
1006 526           return (double)( (logn + logl(logn) - log2 - 0.5) / (log2 + 0.5) );
1007             }