Is the function isprime deterministic?

I was playing with primes, and I was happy to verify that

primes(uint32(2^32))

works. However, I wanted to test a huge prime, and found that (2ˆ64)-59 is prime according to SageMath and Wolfram Alpha. However MATLAB function is prime returns false.

isprime(uint64(2^64)-59)

Which one is right? The documentation of the function isprime seems to suggest that it is deterministic, but it would be nice to say it explicitly. Anyone would, please, confirm it?

[Rafael S.T. Vieira]

Deterministic means that it tests all viable numbers up to x to know if x is prime or not. A non-deterministic algorithm might use some heuristic as the Miller-Rabin primality test.

Answer 1

[John D'Errico]

Your problem is that number is NOT represented in MATLAB as you think it is. MATLAB CANNOT represent a number as a double exactly when the integer is larger than 2^53-1. Double precision is limited in that respect. And yes, I know you THINK you created the number 2^64-59 as a uint64. You did not.

So, YES, I believe it can be shown that isprime is deterministic within the limits of double precision. However, NO, if you throw numbers at it that are too large, then it will fail. The biggest reason why isprime failed for you was that you created it as a double precision number, AND only then converted to uint64. 2^64 is NOT a number you can trust, because it is created as a double.

There is another problem, even more serious. Can you even compute 2^64, and store it in a uint64 variable? NO!

uint64(2)^64
ans =
  uint64
   18446744073709551615
>> intmax('uint64')
ans =
  uint64
   18446744073709551615

Should 2^64 be an odd number? Of course not, yet as you see, 2^64 overflows even a uint64 variable.

Instead, use sym, which her will be quite efficient.

isprime(sym(2)^64-59)
ans =
  logical
   1

See that I created 2 as a symbolic number, only then raising it to a power.

If I did not have the symbolic toolbox, I could have done this - subtly avoiding any overflows, because I recognize that intmax('uint64') is just 2^64-1.

isprime(intmax('uint64') - 58)
ans =
  logical
   1

Now it works. the symbolic toolbox isprime is MUCH faster, although I expect it is really best called isprobable prime for large integers, since it will be based on a multiple set of calls to a MIller-Rabin test, and possibly other tricks. I recall looking up the basis of the mathematica isprime test, which is similar. A true test for primality that cannot be fooled is more difficult to achieve, and will be MUCH slower running most of the time.

By the way, what number did you test to see if it was prime?

uint64(2^64) - 59
ans =
  uint64
   18446744073709551556
   
intmax('uint64') - 58
ans =
  uint64
   18446744073709551557

You should see clearly why MATLAB thought the former version was not prime. After all, most even numbers are not prime. Almost all of them, in fact.

[Rafael S.T. Vieira]

Thanks, John D'Errico. Great explanation and I understand what I did wrong; I'd forgot about adding the zero (I've tried other constructions of 2^64-59 for is prime, but then they had probably failed because of the double limit as you explained). And I did not know I could use the Symbolic toolbox like that, thanks again.

PS: uint64(2^64) == intmax('uint64').

[Steven Lord]

Your PS is correct but slightly misleading.

2^64 is greater than intmax('uint64'). When you convert a value that is outside the range of values an integer type can hold into that integer type, MATLAB saturates at the boundary of that range as stated on this documentation page.

int8(Inf)
int16(-1e6)

I wouldn't say that Inf is 127, but the closest int8 number to Inf is 127.

[Rafael S.T. Vieira]

Thanks, Steven Lord. I realize I am causing an overflow, but for MATLAB that works and it is shorter than writting intmax('uint64') or 0xffffffffffffffff. That's why I've added the PS.

In accordance with the documentation link that you've shared "For binary operations involving a 64-bit integer array and a scalar double, MATLAB computes the operation as if 80-bit extended-precision arithmetic were used, to prevent loss of precision".

PS: If you could, please, look at the MATLAB's code and just confirm if the isprime test is deterministic or not, I would really appreciate it. :)

[John D'Errico]

In general, MathWorks does not tell some algorithmic details as you are asking to learn. Could they be doing 2^32 divides in this case to test that number to see if it is prime?

tic,isprime(intmax('uint64') - 58),toc
ans =
  logical
   1
Elapsed time is 33.309201 seconds.

tic,isprime(sym(2)^64 - 59),toc
ans =
  logical
   1
Elapsed time is 0.027687 seconds.

It seems entirely possible that the uint64/isprime call is essentially doing trial divides below sqrt(N). With some care you can reduce the number of actual tests needed anyway. I can imagine a variety of tricks I might employ involving a partial sieve and rough numbers to minimize the divides needed. In that 33 seconds of time, only one core was running, so the test was not farmed out to all 8 cores on my machine. However, 33 seconds is a relatively LONG time.

But I'm sorry, you probably won't learn what the code does internally if they do not tell you already.

Your other comment is about using syms. DO THAT!!!!!! The sym/isprime utility is pretty fast. I've used it as far out as 1e150000 so far, though it is getting a little slow out there (as much as an hour per test.) However, it does use a Miller-Rabin test out thatfar. I don't know of any explicit place where MATLAB describes their exact scheme for testing for primality under the hood, the reasoning being that it may change at some point.

By way of example though, I just did a search for what Mathematica does.

"The Wolfram Language implements the multiple Rabin-Miller test in bases 2 and 3 combined with a Lucas pseudoprime test as the primality test used by the function PrimeQ[n]. Like many such algorithms, it is a probabilistic test using pseudoprimes. In order to guarantee primality, a much slower deterministic algorithm must be used. However, no numbers are actually known that pass advanced probabilistic tests (such as Rabin-Miller) yet are actually composite."

So it looks like Mathematica's isprime uses Miller-Rabin with two bases, as well as a Lucas-pseudoprime test. This is not uncommon, to use several tests to best insure the result.

Another comparison I might make is to Python, where I also do not know the exact underlying algorithm. But I will state that Python's isprime call is roughly 3x as fast as the sym/isprime. (Sigh. That is important when I have tested numbers on the order of 1e150000.)

So any such test for primality beyond that point will be based most likely on something like MIller-Rabin. There are other things you can do too, in the form of a Lucas based test if you know the full factorization of N+1 to make a test that is fully valid, not just an isprobable prime. Luckily, in my own work that works perfectly. So I can use isprime to tell me when to look a little harder.

[Walter Roberson]

And yes, I know you THINK you created the number 2^64-59 as a uint64. You did not.

>> uint64(2^64)-59
ans =
  uint64
   18446744073709551556

The number is created as uint64, not as double. However because uint64 cannot represent 2^64 it saturates as 2^64-1 and 2^64-1 - 59 is an even number and so not prime.

But you can

isprime(intmax('uint64') - 58)

and it will say true.

[Steven Lord]

isprime is implemented as a MATLAB function file. You can read it and see that there's no calls to any random number generator in it.

[John D'Errico]

@Walter - while A number is created as uint64, THE number that was created was not 2^64-59. I even showed how to create that number, as you did, using intmax. So I don't see the point of your response. I did not claim the number was not uint64, but that it was not the number that was intended to be created. Perhaps you perceived an ambiguity in my wording.

@Steve - I should have checked the code for isprime. I thought I had done so, but that was not the case, since it is provided. (I know that I had looked at the sym version of isprime, and seen that it redirects into symobj::isprime.)

However, in the case of a normal integer/floating point type, isprime does compute all primes less than sqrt(N), and does trial divides, actually trial rems.

[Walter Roberson]

John, you talked specifically about representing as a double, but that was pretty much irrelevant.

uint64(2^64) is calculated as double(2)^double(64) resulting in a exact double 0x43f0000000000000 . The uint64() operator is then applied to that, but the exact double is more than the maximum uint64 so it saturates.

You talk about "Double precision is limited in that respect. And yes, I know you THINK you created the number 2^64-59 as a uint64. You did not. " . Most anyone is going to understand you as indicating that the result of the user's uint64(2^64)-59 is a double instead of uint64. The reason why the result is not 2^64-59 has nothing to do with double precision .

[Rafael S.T. Vieira]

Thanks, Walter Roberson and Steven Lord. The function is indeed deterministic, and I was mistaken in my calculations.