\\ Generation of primes in pari

\\ The first numprimes primes
numprimes=20000;

\\ List of primes
primelist=vector(numprimes);

\\ Last found multiple of corresponding prime
multiples=vector(numprimes);

\\ First candidate prime (priming the prime list!)
candidate=2;

found=0;

\\ Where the squares in the list start
\\ To begin with all are "virtual" squares
\\ (i.e. will be initialised to be squares)
squarestart=1;

\\ Using Eratosthenes sieve
sieve(n)=
{
	local(i,c);

	i=1;c=0;

	while(i<squarestart,
		if(n==multiples[i],
			c=1;
			multiples[i]=multiples[i]+primelist[i];
		);
		i=i+1;
	);

	\\ Now check whether we have a square
	if(n==multiples[squarestart],
		c=1;
		multiples[squarestart]=multiples[squarestart]+primelist[squarestart];
		squarestart=squarestart+1;
	);
	
	\\ Add to the list
	if( c==0,
		found=found+1;
		primelist[found]=n;
		multiples[found]=n*n;
	);
}

while(found<numprimes, sieve(candidate); candidate=candidate+1);

halfprimediffs=vector(numprimes);

i=0;lastp=0;while(i<numprimes,i=i+1;halfprimediffs[i]=(primelist[i]-lastp)/2;lastp=primelist[i]);

write(primesave,printp(halfprimediffs));