[QUOTE=3.14159;222408].. I can't code it altogether. It has to be done in separate tests for the original.[/QUOTE]

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

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

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

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

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

Here you go:
[code]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
};[/code]

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

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
[code]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).");[/code]

[QUOTE=3.14159;222414]Here you go:
[code]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
};[/code][/QUOTE]

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.

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

[QUOTE=3.14159;222418]The test keeps in all the damn pseudoprimes :no: [/QUOTE]

Sure. IF you don't want that, you'll need some re-coding. Split isSPRP into its own function, then do
[code]if(!isSPRP(n), return(0));
if(isprime(n), return(n), return(0))[/code]
after the Fermat test.

Also, as you already know but probably don't realize (:smile:), the last result of the program is simply returned. So my code snippet can be simplified to
[code]if(!isSPRP(n), return(0));
if(isprime(n), n, 0)[/code]
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
[code]if(!isSPRP(n), return(0));
if(isprime(n), n)[/code]

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.

pseudo(n) =!isprime(n)&isSPRP[COLOR="Red"](n)[/COLOR], surely?

I'll just leave this link to Sloane's [url=http://oeis.org/classic/A001262]A001262[/url] here in case you need it.

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

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.