Wheel factorization: Efficient?

[CRGreathouse]

Originally Posted by 3.14159

.. I can't code it altogether. It has to be done in separate tests for the original.

OK. Can you make the fixes I (and sm88) suggested to your test?

[3.14159]

OK. Can you make the fixes I (and sm88) suggested to your test?

Okay. Fixing code, then running some tests. It works.

Results: It lets all the pseudoprimes in. (Tested the well-known example, 1373653)

[CRGreathouse]

Let's see the code! I'll post my version and we can compare.

[3.14159]

Here you go:

primetest(n, b=2) = {
if(ispower(n) == 1,return(0));
forprime(p=2,nextprime(500),
if(n%p==0,return(p))
);
if(Mod(b, n)^(n-1) == 1,return(n));
my(s=valuation(n-1,2),d=n>>s);
n = Mod(b, n)^d;
if (n == 1, return(n-1));
while(s--,
if (n == -1, return(n));
n = n^2;
);
n == -1
};

I forgot to add the base component. It fails for pseudoprimes that pass trial division.

[CRGreathouse]

I haven't run yours, but it will fail on any prime larger than nextprime(501), since you haven't defined your base.

I have

primetest(n:int,b:int=2) = {
	my(s,d);
	if(ispower(n),return(0));
	forprime(p=2,500,
		if(n%p==0,return(p))
	);
	if(Mod(b, n)^(n-1) == 1, return(n));
	s=valuation(n-1,2);
	d=n>>s;
	n = Mod(b, n)^d;
	if (n == 1, return(1));
	while(s--,
		if (n == -1, return(1));
		n = n^2;
	);
	n == -1
};
addhelp(primetest, "primetest(n, {b=2}): Returns the smallest prime factor of n, if one is found. Otherwise, returns 0 if the number is composite and 1 if the number passes all tests (in which case it may be prime or pseudoprime).");

[science_man_88]

Originally Posted by 3.14159

Here you go:

primetest(n) = {
if(ispower(n) == 1,return(0));
forprime(p=2,nextprime(500),
if(n%p==0,return(p))
);
if(Mod(b, n)^(n-1) == 1,return(n));
my(s=valuation(n-1,2),d=n>>s);
n = Mod(b, n)^d;
if (n == 1, return(n-1));
while(s--,
if (n == -1, return(n));
n = n^2;
);
n == -1
};

not sure you need the ==1 one part for ispower as if(ispower(n) means if it returns 1 anyways if(!ispower(n) would return 1 if it was false but if !ispower is true when declaring ispower(n) would declare it flase and stop the execution anyways.

[3.14159]

The test keeps in all the damn pseudoprimes Also: Code fixed, == 1 was removed from nth power check. New code = if(ispower(n), return(0));

[CRGreathouse]

Originally Posted by 3.14159

The test keeps in all the damn pseudoprimes

Sure. IF you don't want that, you'll need some re-coding. Split isSPRP into its own function, then do

if(!isSPRP(n), return(0));
if(isprime(n), return(n), return(0))

after the Fermat test.

Also, as you already know but probably don't realize (), the last result of the program is simply returned. So my code snippet can be simplified to

if(!isSPRP(n), return(0));
if(isprime(n), n, 0)

since you won't be doing anything after that. And actually, since an omitted isFalse part of the if statement returns 0 by default, it could even be

if(!isSPRP(n), return(0));
if(isprime(n), n)

[3.14159]

I'm going to check if the SPRP test I implemented as any pseudoprime numbers:

pseudo(n) =!isprime(n)&isSPRP

forstep(n=3,1e8,2,if(pseudo(n),print(n))

It just prints me every odd number. Damned thing never does anything right.

Let's try:

if(pseudo(n),print(n))

Meh. Manual search.

[CRGreathouse]

pseudo(n) =!isprime(n)&isSPRP(n), surely?

I'll just leave this link to Sloane's A001262 here in case you need it.

[3.14159]

I'll just leave this link to Sloane's A001262 here in case you need it.

No, dammit. Those all fail the TD test. A pseudoprime to my test = Passes TD to 1000003 + Passes SPRP test (No base specified)

Pseudoprime number found: 58398839626753 = 5403649 * 10807297. (Base 2 only)

Wait.. To weed out pseudoprimes
I should change b to (random(n)).

Success: Now, random bases are tested.

Now, to test the same integer length range I did yesterday, to compare timings:
≈ 1000 digits: 312 ms.
≈ 2000 digits: 1.28 seconds
≈ 3000 digits: 3.41 seconds
≈ 13000 digits (second largest prime I've ever found): 134.97 seconds.