Expand Messages
• Feb 8 9:24 PM
• 0 Attachment
On Wednesday 09 February 2005 01:59, you wrote:
> Can you combine this fast caching method with using the prime factorization
> of n as you suggest (except for initialization if you don't start at n=1)?
> I considered the factorization method but thought my suggestion would be
> much faster.

I coded the sieve (code reproduced below). Despite the fact that it segfaults
with gp2c, so that I had to run it interpreted instead of compiled, it's
still a speed demon. With 10 minutes of CPU time (P4 2.6 GHz here) I
eliminated nearly 4000 composites. I shudder to think that this could
probably have cost me somewhere in the vicinity of 10 CPU-days.

I'm sure I commited a few egregious errors down there, not to mention some
ugly hacks, so please refrain from flaming me (:

I wonder, at what point should I switch back to trial factoring?

Décio

---
primorial(x) =
{
local(a);
a = 1;
forprime(p=2,x,a*=p);
return(a);
}

\\ Reverse indexes the primes -- the inverse of PARI's prime() function
init_invprime(nmax) =
{
local(last,k);
global(invprime);
invprime = vector(nmax);
last = 1;
k = 0;
forprime(p=2,nmax,
k++;
for(i=last,p-1,
invprime[i] = k-1;
);
last = p;
);
}

\\ Computes the prime factorization of n!/n#
legendre_factorial(n) =
{
local(v,i,k);
\\ v[] is initialized with -1 to simulate division by n#
v = vector(invprime[n],i,-1);
i = 0;
forprime(p=2,n,
i++;
k = 1;
while((cont = floor(n/p^k)) > 0,
v[i] += cont;
k++;
);
);
return(v);
}

sieve(nmin,nmax,pmax) =
{
local(v,i,j,a,b,apb,amb);
\\ Integers of the form n! +/- n# +/- 1 are not divisible by any
\\ prime up to n, and since the minimum value of n that we're
\\ considering is nmin, then we'll start there. Sure, we'll waste
\\ a little time checking whether the primes from nmin to a given n
\\ divide n! +/- n# +/- 1, for n > nmin, but that's not too costly
\\ anyway.
j = 0;
forprime(p=nmin,pmax,
if(bitand(j++,255)==0,
print("Sieving "p"...");
);
v = legendre_factorial(nmin);

\\ a == nmin!/nmin# (mod p)
a = Mod(1,p);
i = 0;
forprime(q=2,floor(nmin/2),
a *= Mod(q,p)^v[i++];
);

\\ b == nmin# (mod p)
b = Mod(1,p);
forprime(q=2,nmin,
b *= Mod(q,p);
);

\\ Now a == nmin! (mod p)
a *= b;

for(i=nmin,nmax,
apb = a + b;
amb = a - b;

if(apb == Mod(p-1,p),
write("composite_"nmin"_"nmax, i"! + "i"# + 1, "p);
);
if(apb == Mod(1,p),
write("composite_"nmin"_"nmax, i"! + "i"# - 1, "p);
);
if(amb == Mod(p-1,p),
write("composite_"nmin"_"nmax, i"! - "i"# + 1, "p);
);
if(amb == Mod(1,p),
write("composite_"nmin"_"nmax, i"! - "i"# - 1, "p);
);

a *= Mod(i+1,p);
if(isprime(i+1),
b *= Mod(i+1,p);
);
);
);
}

[Non-text portions of this message have been removed]
• Show all 9 messages in this topic