| hermann on Tue, 06 Jan 2026 21:00:02 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| Re: isprime() runtime jump between 10^18 and 10^20? |
On 2026-01-06 11:59, Aurel Page wrote:
Thanks for the explanation. But long is a signed type, shouldn't the transition from long to t_INT happen at +2^63?Dear Hermann, There are two prime testing functions in pari: - ispseudoprime: a fast pseudo-primality test, which in particular guarantees primality of numbers up to 2^64. - isprime: a certified primality test, which only runs ispseudoprime for numbers up to 2^64, but performs a much more costly test otherwise. In addition, as you know you should expect operations on integers that do not fit in a long to be much slower, as we immediately switch to multiprecision integers (t_INT).I did a quick test on my laptop (don't consider these numbers as universal):- ispseudoprime was about 17 times slower for numbers between 2^63 and 2^64 than for numbers between 2^64 and 2^65 (long -> t_INT). - isprime was about 12 times slower than ispseudoprime for numbers between 2^64 and 2^65.Together, these make a jump of x200 for isprime from [2^63,2^64] to [2^64,2^65].Note that 2^64 ~ 10^19, corresponding to your observed jump.
I don't know why I did not think on parallel execution of GP script before, now I did. And because of the results I changed sum() in new oeis sequence to parsum() because that works for single as well as pthread gp engine.
I found the 7 parallel "par*" GP statements and searched for their use on oeis.org. No hit for 5, one (false) hit for "parfor" (Maple code) and only two hits for "parsum".
Runtime reduced from 5min to 5 seconds, and from 23h to 19 minutes using parsum()!
hermann@x3950-X6:~$ nproc 384 hermann@x3950-X6:~$ numactl -C 0-191 gp -q ? a(n)=parsum(k=1, 10^n, isprime(k^2+(k+1)^2)); ? # timer = 1 (on) ? a(9) cpu time = 14min, 51,150 ms, real time = 4,819 ms. 68588950 ? a(10) cpu time = 59h, 57min, 42,370 ms, real time = 18min, 46,398 ms. 614983774 ?In addition I did the parsum() count for 10^10 < k <= 2*10^10, and now the 19 minutes with lots of ispseudoprime() for processing 10^10 ks became 1h (not too bad) for processing the same amount of ks with zero ispseudoprime():
? parsum(k=10^10+1,2*10^10,isprime(k^2+(k+1)^2)); cpu time = 191h, 45min, 10,976 ms, real time = 1h, 3,472 ms. ?In my C++ code I was responsible to make sure that OpenMP creates work on all cores. Using the "par*" GP functions is so much easier, thanks to the dev team!
Regards, Hermann.