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:
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.
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?
f(n) = #{n is a k-strong pseudoprime, 1 < k < n}
where n is a composite integer.
The naive Pari code:
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.
Code:
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
};'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?