RE: Datatypes for sieve algorithms

> if you want to save more, then you can pack 64 bit in each vecsmall
> components by using bitand/bitor/bittest

The C program which uses char vectors can sieve 10^8 in 2 seconds,

while this needs 45 seconds for 10^7, because of the modulo operations

in tst() and clr() which need work and/or gp2c.

_______________________________________________

\\ Sieve of Erastothenes with bit vectors in GP

NB = 64; /* Bits per limb */

ALL = 2^NB-1; /* All limb bits set */

tst(n) = my([k,j]=divrem(n+NB,NB)); bittest(P[k],j);

clr(n) = my([k,j]=divrem(n+NB,NB)); P[k] = bitand(P[k],ALL-2^j);

sieve(N) = { P = vector(1+ceil(N/NB),k,ALL); clr(1);

for(k=2,N/2, if(tst(k), for(j=2,N/k, clr(j*k) )));}

/* The test */

ck(N=2^20) = sieve(N); for(n=1,N, \

if(tst(n) != isprime(n), error("sieve: n = ",n)));

_______________________________________________

Jacques Gélinas

De : Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>
Envoyé : 30 mars 2018 15:33
À : pari-users@pari.math.u-bordeaux.fr
Objet : Re: Datatypes for sieve algorithms

On Sat, Mar 24, 2018 at 10:49:00PM +0000, Jacques Gélinas wrote:
> { The algorithm in question is the sieve of S.P. Sundaram, an Indian
> mathematician who discovered it in 1934. It is proven correct
> in http://en.wikipedia.org/wiki/Sieve_of_Sundaram }
>
> Still, the sieve implementation in PARI/GP raises some
> questions for me, particularly concerning memory.
> While PARI/GP only needs 16MB for 1,048,576 primes,
>
> NP = 2^20;
> GP = primes_zv(NP);
>
> the working vector P has 8,501,423 elements
> and the stack now grows to 128MB.
>
> N = (GP[NP]-1)/2;
> P = vectorsmall(N,n,2*n+1);
>
> However, only one bit per integer is in fact needed to test
> the Sundaram algorithm, which means only ~1MB.
> The sieve operations should involve offset calculations
> and bitset/bittest operations only.
>
> my(l=1,m); for(k=1,(N-1)/3, l+=2; m=k; for(j=1,(N-k)/l, m+=l; P[m] = 0 ));
> MP = N + 1 - vecsum(apply(n->!n,Vec(P))); \\ 256MB !!!

you can replace this line by
MP = N + 1 - sum(i=1,#P,!P[i]);
to save memory.

if you want to save more, then you can pack 64 bit in each vecsmall
components by using bitand/bitor/bittest, but at this point it is
simpler to write the code in C.

Cheers,
Bill.