#restart; read `d:/SPprime.txt`;

#AllSPprime
#input: positive number N and base b
#output: set of positive numer N that is 
# a strong pseudoprime to the base b.

AllSPprime := proc(N,b) local i,S;

S := {};

for i from 4 to N do

	if isprime(i) = false and i mod 2 = 1 
	and SPprime(i,b) = true then
		S := S union {i};
	fi:
od:

return(S);
end:

#SPprime
# input: odd composite integer n, 
# output: true/false

SPprime := proc(n,b) local t,r;

t := n-1;

while type(t,odd)=false do
	t := t/2;
od:


if b^t mod n = 1 then
	return(true);
fi:

for r from 0 to log[2]((n-1)/t)-1 do 

	if b^((2^r)*t) mod n = n-1 then 
		return(true);
	fi:
od:

return(false);
end: