Re: Datatypes for sieve algorithms

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Fri, Mar 23, 2018 at 08:44:55PM +0000, Jacques Gélinas wrote:
> A prime sieve algorithm is proposed in http://vixra.org/pdf/1803.0493v1.pdf
> 
> From the list of even integers,                          [2,4,6,8,10,12,14,16,18,20,22,...]
> eliminate 4+2r + n(2+4r) where r,n=1,2,3,...  [12,18,24,30,...],  [18,28,38,...], ... 
> then subtract 3 from those remaining to get  [_,1,3,5, 7, _,11,13, _,17,19,_,...]
> The list of primes below 4(r^2+r+1) is said to be complete as this number is eliminated.
> 
> In order to use PARI/GP to test this (unproven++) algorithm, what kind of datatypes/structures
> are available and efficient ? Vectors of GP integers seem to me to be wasteful here !!
> What is needed is an index for primes that could be used in vecextract([1..n]).

You could try to use Vecsmall/vectorsmall. At least this would save some
memory.

But honestly, sieves are simpler to implement in C in general.

> for(k=1,N/3-1, for(j=1,(N-2-k)/(1+2*k), P[k+j*(1+2*k)] = 0 ));

in GP this will not make much of a difference, but in a sieve in
general, multiplications needed to be avoided and replaced by additions.
In fact this is the source of the efficiency of a sieve.
You can do:

my(l=1,m);for(k=1,N/3-1, l+=2; m=k; for(j=1,(N-2-k)/l, m+=l; P[m] = 0));

Please try the file in attachment.

Cheers,
Bill.

install(primes_zv,L)
test(e)=
{
NP = 2^e;               /* number of primes  */
GP = primes_zv(NP);
N  = (GP[NP]+3)/2;  /* largest index needed in P */

P  = vectorsmall(N-2,n,2*n+1);
my(l=1,m);for(k=1,N/3-1, l+=2; m=k; for(j=1,(N-2-k)/l, m+=l; P[m] = 0 ));
MP = N - 1 - vecsum(apply(n->!n,Vec(P)));

Q  = vectorsmall(MP); Q[1]=2;     /* drop zeros from P */
m=1; for(n=1,N-2,if(P[n],Q[m++]=P[n]));

Q == GP                 /* the test */
}
test(20);
##