What techniques can speed up counting pseudoprimes?

[CRGreathouse]

I'm calculating the high-water marks for the following function:
f(n) = #{n is a k-strong pseudoprime, 1 < k < n}
where n is a composite integer.

The naive Pari code:

ff(n)=sum(k=2,n-1,isSPRP(n, k))

record=0;forstep(n=3, 1e6, 2, if(isprime(n), next); k = ff(n); if (k > record, record = k; print (k" "n)))

Note: I'm using the following straightforward code for checking if a number is a strong pseudoprime; here's the code if you want to duplicate my above functions. I'm not looking for improvements to this part, though if you have any I'd love to hear about it.

isSPRP(n,b)={
	my(d = n, s = 0);
	until(bitand(d,1), d >>= 1; s++);
	d = Mod(b, n)^d;
	if (d == 1, return(1));
	for(i=1,s-1,
		if (d == -1, return(1));
		d = d^2;
	);
	d == -1
};

But this is very slow: checking ff for a number around a million takes about 10 seconds, so I'm looking at weeks or months of calculation even to check to a million (and I'd like to check more than that).

'Everyone' knows that a composite number can't be a strong pseudoprime to more than a quarter of the k, but more is true: ff(n) is at most phi(n)/4. This simple test significantly speeds results, allowing me to skip between 20% and 90% of the numbers in the range (the relatively smooth/abundant ones) -- the further closer together the records are, the more I can skip. Are there any other ideas?

 

[fresh_42]

Maybe this would help:
http://mathworld.wolfram.com/Pseudoprime.htmlhttps://www.maths.lancs.ac.uk/~jameson//psp.pdf