#include "gmp.h"
#include "gmp-impl.h"
#define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve) \
__max_i = (end); \
\
do { \
++__i; \
if (((sieve)[__index] & __mask) == 0) \
{ \
(prime) = id_to_n(__i)
#define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve) \
do { \
mp_limb_t __mask, __index, __max_i, __i; \
\
__i = (start)-(off); \
__index = __i / GMP_LIMB_BITS; \
__mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS); \
__i += (off); \
\
LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)
#define LOOP_ON_SIEVE_STOP \
} \
__mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1); \
__index += __mask & 1; \
} while (__i <= __max_i) \
#define LOOP_ON_SIEVE_END \
LOOP_ON_SIEVE_STOP; \
} while (0)
#if 0#endif
static mp_limb_t
id_to_n (mp_limb_t id) { return id*3+1+(id&1); }
static mp_limb_t
n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
#if 0#endif
#if GMP_LIMB_BITS > 61
#define SIEVE_SEED CNST_LIMB(0x3294C9E069128480)
#define SEED_LIMIT 202
#else
#if GMP_LIMB_BITS > 30
#define SIEVE_SEED CNST_LIMB(0x69128480)
#define SEED_LIMIT 114
#else
#if GMP_LIMB_BITS > 15
#define SIEVE_SEED CNST_LIMB(0x8480)
#define SEED_LIMIT 54
#else
#if GMP_LIMB_BITS > 7
#define SIEVE_SEED CNST_LIMB(0x80)
#define SEED_LIMIT 34
#else
#define SIEVE_SEED CNST_LIMB(0x0)
#define SEED_LIMIT 24
#endif
#endif
#endif
#endif
static void
first_block_primesieve (mp_ptr bit_array, mp_limb_t n)
{
mp_size_t bits, limbs;
ASSERT (n > 4);
bits = n_to_bit(n);
limbs = bits / GMP_LIMB_BITS + 1;
MPN_ZERO (bit_array, limbs);
bit_array[0] = SIEVE_SEED;
if ((bits + 1) % GMP_LIMB_BITS != 0)
bit_array[limbs-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
if (n > SEED_LIMIT) {
mp_limb_t mask, index, i;
ASSERT (n > 49);
mask = 1;
index = 0;
i = 1;
do {
if ((bit_array[index] & mask) == 0)
{
mp_size_t step, lindex;
mp_limb_t lmask;
unsigned maskrot;
step = id_to_n(i);
lindex = i*(step+1)-1+(-(i&1)&(i+1));
if (lindex > bits)
break;
step <<= 1;
maskrot = step % GMP_LIMB_BITS;
lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
do {
bit_array[lindex / GMP_LIMB_BITS] |= lmask;
lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
lindex += step;
} while (lindex <= bits);
lindex = i*(i*3+6)+(i&1);
lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
for ( ; lindex <= bits; lindex += step) {
bit_array[lindex / GMP_LIMB_BITS] |= lmask;
lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
};
}
mask = mask << 1 | mask >> (GMP_LIMB_BITS-1);
index += mask & 1;
i++;
} while (1);
}
}
static void
block_resieve (mp_ptr bit_array, mp_size_t limbs, mp_limb_t offset,
mp_srcptr sieve, mp_limb_t sieve_bits)
{
mp_size_t bits, step;
ASSERT (limbs > 0);
bits = limbs * GMP_LIMB_BITS - 1;
MPN_ZERO (bit_array, limbs);
LOOP_ON_SIEVE_BEGIN(step,0,sieve_bits,0,sieve);
{
mp_size_t lindex;
mp_limb_t lmask;
unsigned maskrot;
lindex = __i*(step+1)-1+(-(__i&1)&(__i+1));
if (lindex > bits + offset)
break;
step <<= 1;
maskrot = step % GMP_LIMB_BITS;
if (lindex < offset)
lindex += step * ((offset - lindex - 1) / step + 1);
lindex -= offset;
lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
for ( ; lindex <= bits; lindex += step) {
bit_array[lindex / GMP_LIMB_BITS] |= lmask;
lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
};
lindex = __i*(__i*3+6)+(__i&1);
if (lindex > bits + offset)
continue;
if (lindex < offset)
lindex += step * ((offset - lindex - 1) / step + 1);
lindex -= offset;
lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
for ( ; lindex <= bits; lindex += step) {
bit_array[lindex / GMP_LIMB_BITS] |= lmask;
lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
};
}
LOOP_ON_SIEVE_END;
}
#define BLOCK_SIZE 2048
mp_limb_t
gmp_primesieve (mp_ptr bit_array, mp_limb_t n)
{
mp_size_t size;
mp_limb_t bits;
ASSERT (n > 4);
bits = n_to_bit(n);
size = bits / GMP_LIMB_BITS + 1;
if (size > BLOCK_SIZE * 2) {
mp_size_t off;
off = BLOCK_SIZE + (size % BLOCK_SIZE);
first_block_primesieve (bit_array, id_to_n (off * GMP_LIMB_BITS));
for ( ; off < size; off += BLOCK_SIZE)
block_resieve (bit_array + off, BLOCK_SIZE, off * GMP_LIMB_BITS, bit_array, off * GMP_LIMB_BITS - 1);
} else {
first_block_primesieve (bit_array, n);
}
if ((bits + 1) % GMP_LIMB_BITS != 0)
bit_array[size-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
return size * GMP_LIMB_BITS - mpn_popcount (bit_array, size);
}
#undef BLOCK_SIZE
#undef SEED_LIMIT
#undef SIEVE_SEED
#undef LOOP_ON_SIEVE_END
#undef LOOP_ON_SIEVE_STOP
#undef LOOP_ON_SIEVE_BEGIN
#undef LOOP_ON_SIEVE_CONTINUE