- Feb 8 9:24 PMOn Wednesday 09 February 2005 01:59, you wrote:
> Can you combine this fast caching method with using the prime factorization

I coded the sieve (code reproduced below). Despite the fact that it segfaults

> 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.

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] - << Previous post in topic Next post in topic >>