incorrect isprime results for large numbers

[waldyrious]

julia> isprime(1000000007) # 10 digits
true

julia> isprime(10000000019) # 11 digits
false

[StefanKarpinski]

This doesn't seem to be a matter of witness choices so something is up with the algorithm itself. Good catch.

[StefanKarpinski]

The powermod function is the culprit. According to Wolfram Alpha, 2^5000000009 = 10000000018 = -1 (mod 10000000019), but powermod(2,5000000009,10000000019) is giving 7653772409 instead.

[JeffBezanson]

One of the intermediate multiplications is overflowing.

[JeffBezanson]

One fix is to do the multiplications with widemul.

[StefanKarpinski]

Hmm. I'd like to see if I can avoid that. Resorting to 128-bit multiply might be slow.

[StefanKarpinski]

Decided that it was more important to fix the broken isprime behavior than to have the absolute fastest possible powermod algorithm. We can always improve it at some point in the future if it becomes a bottleneck for someone.