#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
/******************** functions to support sequences ********************/
#include "sequences.h"
static UV isqrt(UV n) {
UV root;
#if BITS_PER_WORD == 32
if (n >= W_CONST(4294836225)) return W_CONST(65535);
#else
if (n >= W_CONST(18446744065119617025)) return W_CONST(4294967295);
#endif
root = (UV) sqrt((double)n);
while (root*root > n) root--;
while ((root+1)*(root+1) <= n) root++;
return root;
}
unsigned char* sieve_erat30(WTYPE end);
/* Used for moving between primes */
static unsigned char nextwheel30[30] = {
1, 7, 7, 7, 7, 7, 7, 11, 11, 11, 11, 13, 13, 17, 17,
17, 17, 19, 19, 23, 23, 23, 23, 29, 29, 29, 29, 29, 29, 1 };
static unsigned char prevwheel30[30] = {
29, 29, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 11, 11, 13,
13, 13, 13, 17, 17, 19, 19, 19, 19, 23, 23, 23, 23, 23, 23 };
/* The bit mask within a byte */
static unsigned char masktab30[30] = {
0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4, 0, 8, 0,
0, 0, 16, 0, 32, 0, 0, 0, 64, 0, 0, 0, 0, 0,128 };
/* Add this to a number and you'll ensure you're on a wheel location */
static unsigned char distancewheel30[30] = {
1, 0, 5, 4, 3, 2, 1, 0, 3, 2, 1, 0, 1, 0, 3,
2, 1, 0, 1, 0, 3, 2, 1, 0, 5, 4, 3, 2, 1, 0 };
#if 0
static int is_prime_in_sieve(const unsigned char* sieve, WTYPE p) {
WTYPE d = p/30;
WTYPE m = p - d*30;
/* If m isn't part of the wheel, we return 0 */
return ( (masktab30[m] != 0) && ((sieve[d] & masktab30[m]) == 0) );
}
#endif
/* Warning -- can go off the end of the sieve */
static WTYPE next_prime_in_sieve(const unsigned char* sieve, WTYPE p) {
WTYPE d, m;
if (p < 7)
return (p < 2) ? 2 : (p < 3) ? 3 : (p < 5) ? 5 : 7;
d = p/30;
m = p - d*30;
do {
if (m==29) { d++; m = 1; while (sieve[d] == 0xFF) d++; }
else { m = nextwheel30[m]; }
} while (sieve[d] & masktab30[m]);
return(d*30+m);
}
static WTYPE prev_prime_in_sieve(const unsigned char* sieve, WTYPE p) {
WTYPE d, m;
if (p <= 7)
return (p <= 2) ? 0 : (p <= 3) ? 2 : (p <= 5) ? 3 : 5;
d = p/30;
m = p - d*30;
do {
m = prevwheel30[m]; if (m==29) { if (d == 0) return 0; d--; }
} while (sieve[d] & masktab30[m]);
return(d*30+m);
}
/* Useful macros for the wheel-30 sieve array */
#define START_DO_FOR_EACH_SIEVE_PRIME(sieve, a, b) \
{ \
WTYPE p = a; \
WTYPE l_ = b; \
WTYPE d_ = p/30; \
WTYPE m_ = p-d_*30; \
m_ += distancewheel30[m_]; \
p = d_*30 + m_; \
while ( p <= l_ ) { \
if ((sieve[d_] & masktab30[m_]) == 0)
#define END_DO_FOR_EACH_SIEVE_PRIME \
m_ = nextwheel30[m_]; if (m_ == 1) { d_++; } \
p = d_*30+m_; \
} \
}
#if 0
static __inline__ uint64_t rdtsc(void)
{
unsigned a, d;
asm volatile("rdtsc" : "=a" (a), "=d" (d));
return ((uint64_t)a) | (((uint64_t)d) << 32);
}
/* uint64_t ts = rdtsc(); .... te = tdtsc(); tot += te-ts; */
#endif
/******************** primes ********************/
/* GCC 3.4 - 4.1 has broken 64-bit popcount.
* GCC 4.2+ can generate awful code when it doesn't have asm (GCC bug 36041).
* When the asm is present (e.g. compile with -march=native on a platform that
* has them, like Nahelem+), then it is almost as fast as the direct asm. */
#if BITS_PER_WORD == 64
#if defined(__POPCNT__) && defined(__GNUC__) && (__GNUC__> 4 || (__GNUC__== 4 && __GNUC_MINOR__> 1))
#define popcnt(b) __builtin_popcountll(b)
#else
static UV popcnt(UV b) {
b -= (b >> 1) & 0x5555555555555555;
b = (b & 0x3333333333333333) + ((b >> 2) & 0x3333333333333333);
b = (b + (b >> 4)) & 0x0f0f0f0f0f0f0f0f;
return (b * 0x0101010101010101) >> 56;
}
#endif
#endif
#if defined(__GNUC__)
#define word_unaligned(m,wordsize) ((uintptr_t)m & (wordsize-1))
#else /* uintptr_t is part of C99 */
#define word_unaligned(m,wordsize) ((unsigned int)m & (wordsize-1))
#endif
static const unsigned char byte_zeros[256] =
{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,
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,
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,
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,
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,
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,
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,
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};
static WTYPE count_zero_bits(const unsigned char* m, WTYPE nbytes)
{
WTYPE count = 0;
#if BITS_PER_WORD == 64
if (nbytes >= 16) {
while ( word_unaligned(m,sizeof(UV)) && nbytes--)
count += byte_zeros[*m++];
if (nbytes >= 8) {
UV* wordptr = (UV*)m;
UV nwords = nbytes / 8;
UV nzeros = nwords * 64;
m += nwords * 8;
nbytes %= 8;
while (nwords--)
nzeros -= popcnt(*wordptr++);
count += nzeros;
}
}
#endif
while (nbytes--)
count += byte_zeros[*m++];
return count;
}
static unsigned char* prime_cache_sieve = 0;
static WTYPE prime_cache_size = 0;
/*
* Get the size and a pointer to the cached prime sieve.
* Returns the maximum sieved value in the sieve.
* Allocates and sieves if needed.
*
* The sieve holds 30 numbers per byte, using a mod-30 wheel.
*/
static WTYPE get_prime_cache(WTYPE n, const unsigned char** sieve)
{
if (prime_cache_size < n) {
if (prime_cache_sieve != 0)
Safefree(prime_cache_sieve);
prime_cache_size = 0;
/* Sieve a bit more than asked, to mitigate thrashing */
if (n < (W_FFFF-3840))
n += 3840;
/* TODO: testing near 2^32-1 */
prime_cache_sieve = sieve_erat30(n);
if (prime_cache_sieve != 0)
prime_cache_size = n;
}
if (sieve != 0)
*sieve = prime_cache_sieve;
return prime_cache_size;
}
/* Marked bits for each n, indicating if the number is prime */
static const unsigned char prime_is_small[] =
{0xac,0x28,0x8a,0xa0,0x20,0x8a,0x20,0x28,0x88,0x82,0x08,0x02,0xa2,0x28,0x02,
0x80,0x08,0x0a,0xa0,0x20,0x88,0x20,0x28,0x80,0xa2,0x00,0x08,0x80,0x28,0x82,
0x02,0x08,0x82,0xa0,0x20,0x0a,0x20,0x00,0x88,0x22,0x00,0x08,0x02,0x28,0x82,
0x80,0x20,0x88,0x20,0x20,0x02,0x02,0x28,0x80,0x82,0x08,0x02,0xa2,0x08,0x80,
0x80,0x08,0x88,0x20,0x00,0x0a,0x00,0x20,0x08,0x20,0x08,0x0a,0x02,0x08,0x82,
0x82,0x20,0x0a,0x80,0x00,0x8a,0x20,0x28,0x00,0x22,0x08,0x08,0x20,0x20,0x80,
0x80,0x20,0x88,0x80,0x20,0x02,0x22,0x00,0x08,0x20,0x00,0x0a,0xa0,0x28,0x80,
0x00,0x20,0x8a,0x00,0x20,0x8a,0x00,0x00,0x88,0x80,0x00,0x02,0x22,0x08,0x02};
#define NPRIME_IS_SMALL (sizeof(prime_is_small)/sizeof(prime_is_small[0]))
int is_prime(WTYPE n)
{
WTYPE d;
unsigned char mtab;
if ( n < (NPRIME_IS_SMALL*8))
return ((prime_is_small[n/8] >> (n%8)) & 1);
d = n/30;
mtab = masktab30[ n - d*30 ]; /* Bitmask in mod30 wheel */
if (mtab == 0) return 0; /* Return 0 if a multiple of 2, 3, or 5 */
if (n <= prime_cache_size)
return ((prime_cache_sieve[d] & mtab) == 0);
if (!(n%7) || !(n%11) || !(n%13) || !(n%17) || !(n%23) || !(n%29) || !(n%31))
return 0;
{
UV limit = isqrt(n);
UV i = 37;
while (1) { /* trial division, skipping multiples of 2/3/5 */
if (i > limit) break; if ((n % i) == 0) return 0; i += 4;
if (i > limit) break; if ((n % i) == 0) return 0; i += 2;
if (i > limit) break; if ((n % i) == 0) return 0; i += 4;
if (i > limit) break; if ((n % i) == 0) return 0; i += 2;
if (i > limit) break; if ((n % i) == 0) return 0; i += 4;
if (i > limit) break; if ((n % i) == 0) return 0; i += 6;
if (i > limit) break; if ((n % i) == 0) return 0; i += 2;
if (i > limit) break; if ((n % i) == 0) return 0; i += 6;
}
}
return 1;
}
static const unsigned char prime_count_small[] =
{0,0,1,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,8,8,8,8,9,9,9,9,9,9,10,10,
11,11,11,11,11,11,12,12,12,12,13,13,14,14,14,14,15,15,15,15,15,15,
16,16,16,16,16,16,17,17,18,18,18,18,18,18,19};
#define NPRIME_COUNT_SMALL (sizeof(prime_count_small)/sizeof(prime_count_small[0]))
static const unsigned char primes_small[] =
{0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71};
#define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))
/* The nth prime will be less than or equal to this number */
static UV nth_prime_upper(WTYPE n)
{
double fn, flogn, flog2n, upper;
if (n < NPRIMES_SMALL)
return primes_small[n];
fn = (double) n;
flogn = log(n);
flog2n = log(flogn); /* Note distinction between log_2(n) and log^2(n) */
if (n >= 688383) /* Dusart 2010 page 2 */
upper = fn * (flogn + flog2n - 1.0 + ((flog2n-2.00)/flogn));
else if (n >= 178974) /* Dusart 2010 page 7 */
upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.95)/flogn));
else if (n >= 39017) /* Dusart 1999 page 14 */
upper = fn * (flogn + flog2n - 0.9484);
else if (n >= 6) /* Modified from Robin 1983 for 6-39016 _only_ */
upper = fn * ( flogn + 0.6000 * flog2n );
else
upper = fn * ( flogn + flog2n );
/* Watch out for overflow */
if (upper >= (double)UV_MAX) {
#if BITS_PER_WORD == 32
if (n <= W_CONST(203280221)) return W_CONST(4294967291);
#else
if (n <= W_CONST(425656284035217743)) return W_CONST(18446744073709551557);
#endif
croak("nth_prime_upper(%"UVuf") overflow", n);
}
return (WTYPE) ceil(upper);
}
UV nth_prime(WTYPE n)
{
const unsigned char* sieve;
UV upper_limit, start, count, s, bytes_left;
if (n < NPRIMES_SMALL)
return primes_small[n];
upper_limit = nth_prime_upper(n);
if (upper_limit == 0) {
croak("nth_prime(%lu) would overflow", (unsigned long)n);
return 0;
}
/* The nth prime is guaranteed to be within this range */
if (get_prime_cache(upper_limit, &sieve) < upper_limit) {
croak("Couldn't generate sieve for nth(%lu) [sieve to %lu]", (unsigned long)n, (unsigned long)upper_limit);
return 0;
}
count = 3;
start = 7;
s = 0;
bytes_left = (n-count) / ((n<24000)?8:(n<3000000)?4:3);
while ( bytes_left > 0 ) {
/* There is at minimum one byte we can count (and probably many more) */
count += count_zero_bits(sieve+s, bytes_left);
assert(count <= n);
s += bytes_left;
bytes_left = (n-count) / 8;
}
if (s > 0)
start = s * 30;
START_DO_FOR_EACH_SIEVE_PRIME(sieve, start, upper_limit)
if (++count == n) return p;
END_DO_FOR_EACH_SIEVE_PRIME;
croak("nth_prime failed for %lu, not found in range %lu - %lu", (unsigned long)n, (unsigned long) start, (unsigned long)upper_limit);
return 0;
}
void prime_init(WTYPE n)
{
if ( (n == 0) && (prime_cache_sieve == 0) ) {
/* On init, make a few primes (2-30k using 1k memory) */
size_t initial_primes_to = 30 * (1024-8);
prime_cache_sieve = sieve_erat30(initial_primes_to);
if (prime_cache_sieve != 0)
prime_cache_size = initial_primes_to;
return;
}
get_prime_cache(n, 0); /* Sieve to n */
}
UV prime_count(WTYPE n)
{
const unsigned char* sieve;
static WTYPE last_bytes = 0;
static UV last_count = 3;
WTYPE s, bytes;
UV count = 3;
if (n < NPRIME_COUNT_SMALL)
return prime_count_small[n];
/* Get the cached sieve. */
if (get_prime_cache(n, &sieve) < n) {
croak("Couldn't generate sieve for prime_count");
return 0;
}
#if 0
/* The really simple way -- walk the sieve */
START_DO_FOR_EACH_SIEVE_PRIME(sieve, 7, n)
count++;
END_DO_FOR_EACH_SIEVE_PRIME;
#else
bytes = n / 30;
s = 0;
/* Start from last word position if we can. This is a big speedup when
* calling prime_count many times with successively larger numbers. */
if (bytes >= last_bytes) {
s = last_bytes;
count = last_count;
}
count += count_zero_bits(sieve+s, bytes-s);
last_bytes = bytes;
last_count = count;
START_DO_FOR_EACH_SIEVE_PRIME(sieve, 30*bytes, n)
count++;
END_DO_FOR_EACH_SIEVE_PRIME;
#endif
return count;
}
/* Crude way to get this for 7 or 11:
perl -E 'my $n = "0" x 210; for ($s=7; $s<210; $s+=7) { substr($n,$s,1,"1"); } for $s (0..length($n)-1) { $b .= substr($n,$s,1) if $s%2 && $s%3 && $s%5 } say join ",", map { sprintf "0x%02x", oct("0b".reverse(substr($b,$_*8,8))); } 0..6'
perl -E 'my $n = "0" x 2310; for ($s=7; $s<2310; $s+=7) { substr($n,$s,1,"1"); } for ($s=11; $s<2310; $s+=11) { substr($n,$s,1,"1"); } for $s (0..length($n)-1) { $b .= substr($n,$s,1) if $s%2 && $s%3 && $s%5 } say join ",", map { sprintf "0x%02x", oct("0b".reverse(substr($b,$_*8,8))); } 0..(7*11-1)'
*/
#define PRESIEVE_SIZE (7*11)
static const unsigned char presieve11[PRESIEVE_SIZE] =
{ 0x06,0x20,0x10,0x81,0x49,0x04,0xc2,0x02,0x28,0x10,0xa1,0x0c,0x04,0x50,0x02,0x61,0x10,0x83,0x08,0x0c,0x40,0x22,0x24,0x10,0x91,0x08,0x45,0x40,0x82,0x20,0x18,0x81,0x28,0x04,0x40,0x12,0x20,0x51,0x81,0x8a,0x04,0x48,0x02,0x20,0x14,0x81,0x18,0x04,0x41,0x02,0xa2,0x10,0x89,0x08,0x24,0x44,0x02,0x30,0x10,0xc1,0x08,0x86,0x40,0x0a,0x20,0x30,0x85,0x08,0x14,0x40,0x43,0x20,0x92,0x81,0x08,0x04,0x60 };
static void memtile(unsigned char* src, UV from, UV to) {
while (from < to) {
UV bytes = (2*from > to) ? to-from : from;
memcpy(src+from, src, bytes);
from += bytes;
}
}
static UV sieve_prefill(unsigned char* mem, UV startd, UV endd)
{
UV nbytes = endd - startd + 1;
if (nbytes > 0) {
memcpy(mem, presieve11, (nbytes < PRESIEVE_SIZE) ? nbytes : PRESIEVE_SIZE);
memtile(mem, PRESIEVE_SIZE, nbytes);
if (startd == 0) mem[0] = 0x01; /* Correct first byte */
}
return 13;
}
/* Wheel 30 sieve. Ideas from Terje Mathisen and Quesada / Van Pelt. */
unsigned char* sieve_erat30(WTYPE end)
{
unsigned char* mem;
WTYPE max_buf, limit, prime;
max_buf = (end/30) + ((end%30) != 0);
/* Round up to a word */
max_buf = ((max_buf + sizeof(UV) - 1) / sizeof(UV)) * sizeof(UV);
New(0, mem, max_buf, unsigned char);
/* Fill buffer marked with small primes 7+ */
prime = sieve_prefill(mem, 0, max_buf-1);
limit = isqrt(end); /* prime*prime can overflow */
for ( ; prime <= limit; prime = next_prime_in_sieve(mem,prime)) {
WTYPE d = (prime*prime)/30;
WTYPE m = (prime*prime) - d*30;
WTYPE dinc = (2*prime)/30;
WTYPE minc = (2*prime) - dinc*30;
WTYPE wdinc[8];
unsigned char wmask[8];
int i;
/* Find the positions of the next composites we will mark */
for (i = 1; i <= 8; i++) {
WTYPE dlast = d;
do {
d += dinc;
m += minc;
if (m >= 30) { d++; m -= 30; }
} while ( masktab30[m] == 0 );
wdinc[i-1] = d - dlast;
wmask[i%8] = masktab30[m];
}
d -= prime;
#if 0
assert(d == ((prime*prime)/30));
assert(d < max_buf);
assert(prime = (wdinc[0]+wdinc[1]+wdinc[2]+wdinc[3]+wdinc[4]+wdinc[5]+wdinc[6]+wdinc[7]));
#endif
i = 0; /* Mark the composites */
do {
mem[d] |= wmask[i];
d += wdinc[i];
i = (i+1) & 7;
} while (d < max_buf);
}
return mem;
}
/******************** best pair (for Additive) ********************/
static int gamma_length(WTYPE n)
{
#if defined(__GNUC__) && (__GNUC__ >= 4 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4))
#if BITS_PER_WORD == 64
WTYPE log2 = 63 - __builtin_clzll(n+1);
#else
WTYPE log2 = 31 - __builtin_clzl(n+1);
#endif
#else
WTYPE log2 = 0;
while (n >= ((2 << log2)-1)) log2++;
#endif
return ((2*log2)+1);
}
/* adder is used to modify the stored indices. A function would be better. */
int find_best_pair(WTYPE* basis, int basislen, WTYPE val, int adder, int* a, int* b)
{
int maxbasis;
int bestlen = INT_MAX;
int i, j;
assert( (basis != 0) && (a != 0) && (b != 0) && (basislen >= 1) );
/* Find how far in basis to look */
if ((basislen > 15) && (val > basis[15])) {
/* Binary search for large values */
i = 0;
j = basislen-1;
while (i < j) {
int mid = (i+j)/2;
if (basis[mid] < val) i = mid+1;
else j = mid;
}
maxbasis = i-1;
} else {
/* Iteration for small values */
maxbasis = 0;
while ( ((maxbasis+1) < basislen) && (basis[maxbasis+1] < val) )
maxbasis++;
}
assert(maxbasis < basislen);
assert(basis[maxbasis] <= val);
assert( ((maxbasis+1) == basislen) || (basis[maxbasis+1] >= val) );
i = 0;
j = maxbasis;
while (i <= j) {
WTYPE sum = basis[i] + basis[j];
if (sum > val) {
j--;
} else {
if (sum == val) {
int p1 = i + adder;
int p2 = j - i + adder;
int glen = gamma_length(p1) + gamma_length(p2);
/* printf("found %llu+%llu=%llu pair %d,%d (%d,%d) with length %d\n", basis[i], basis[j], sum, i, j, p1, p2, glen); */
if (glen < bestlen) {
*a = p1;
*b = p2;
bestlen = glen;
}
}
i++;
}
}
return (bestlen < INT_MAX);
}
/* If you roll your own prev_prime and next_prime, you can make this
* about 35% faster. I decided it wasn't worth the obfuscation. E.g.
*
* if (i <= 3) { pim = pi = (i==1) ? 3 : (i==2) ? 5 : 7;
* } else { do { pim = nextwheel30[pim]; if (pim == 1) pid++;
* } while (sieve[pid] & masktab30[pim]);
* pi = pid*30+pim; }
*/
int find_best_prime_pair(WTYPE val, int adder, int* a, int* b)
{
int bestlen = INT_MAX;
int i, j;
WTYPE pi, pj;
const unsigned char* sieve;
assert( (a != 0) && (b != 0) );
if (get_prime_cache(val, &sieve) < val) {
croak("Couldn't generate sieve for find_best_prime_pair");
return 0;
}
pi = 1;
pj = prev_prime_in_sieve(sieve,val+1);
i = 0;
j = (val <= 2) ? 1 : prime_count(pj)-1;
while (i <= j) {
WTYPE sum = pi + pj;
if (sum > val) {
j--;
pj = (j == 0) ? 1 : prev_prime_in_sieve(sieve,pj);
} else {
if (sum == val) {
int p1 = i + adder;
int p2 = j - i + adder;
int glen = gamma_length(p1) + gamma_length(p2);
if (glen <= bestlen) { /* Prefer a smaller j */
*a = p1;
*b = p2;
bestlen = glen;
}
}
i++;
pi = (i == 1) ? 3 : next_prime_in_sieve(sieve,pi);
}
}
return (bestlen < INT_MAX);
}