Test this out and see if you can find the flaw in the program
Clear[a, b, p]; k = 2^3217 - 1; p =
Solve[{(2 (1 + 2 k) (-k + a (1 + k - b)) b)/((-1 + a b) (2 k +
a b) (1 + 2 k + a b)) == 0, a > 0, b > 0}, {a, b},
Integers]; p = FromDigits[Differences[Flatten[{a, b} /. p]]]; If[
p == 0, Print[k, " <-Is Prime"], Print[k, " <- Is Not Prime"]]
try some values for k.
The equation is needlessly complicated, but that’s not the main issue. Try it with k=2^327-1 and you will see it is on the slow side. The reason is Solve is basically factoring your k. That’s not a great way to show primality…
Yes Daniel I realise it’s not the way, a slight variation of the program above was written to solve a combination type problem, it was because it was behaving like it was testing the primality of numbers that made me think a little. It wasn’t until I had investigated further that I found out how it was doing it. If with a slight alteration you Reduce the equation above you will get.
k == (-a + a b)/(-1 + a)
which is identical to just finding the divisors of k, as you point out. It was just a realisation on my part that is how Solve must work in certain situations.