Searching Mersenne Primes with Mathematica

[dcheuk]

Originally Posted by mattprim

Mathematica has the native prime Number Theory routines: PrimeQ, NextPrime, PrimePi etc. Just was mainly curious how PrimeQ[2^n-1] giving the answer if 2^n-1 is prime. You can actually link stuff like PariGP to Mathematica using Mathlink so you can include own faster routines for Mersenne.

For science I picked 3 small M_p, of which one is prime, and run it with Mathematica and Prime95 with LL.

Unsurprising, Prime95 is exponentially faster since it is optimized to specifically test M_p for primality.

[mattprim]

Lucas-Lehmer test in native Mathematica:

p = 2976221;
s = 4;
Mp = 2^p - 1;
Do[s = Mod[(s*s - 2) , Mp];
If[s == 0, Print["PRIME"]], {i, 0, p - 2}];

puts me in top prime guys fast at least in 60-ties (below GHz*Day).

[Uncwilly]

Originally Posted by mattprim

p = 2976221;

That is 30 times smaller than the current wave front.
It got double checked more than 20 years ago.

puts me in top prime guys fast at least in 60-ties (below GHz*Day).

Top prime guys? What do you mean by that statement?

[LaurV]

Originally Posted by Uncwilly

Top prime guys? What do you mean by that statement?

He means the year. His test is as fast as some supercomputer ~55 years ago . Albeit I have to admit I don't know how to pronounce the the thing he wrote... sixty ties? sixty tits?

[LaurV]

Originally Posted by mattprim

Lucas-Lehmer test in native Mathematica:

p = 2976221;
s = 4;
Mp = 2^p - 1;
Do[s = Mod[(s*s - 2) , Mp];
If[s == 0, Print["PRIME"]], {i, 0, p - 2}];

puts me in top prime guys fast at least in 60-ties (below GHz*Day).

Can you prove this test is correct? I mean, [excluding the fact that you could make it much faster by skipping the initialization (see the syntax of the "do" in Wolfram), and by taking the "if" out of the loop], can you prove that for some composite, you don't hit zero on the way? Can you prove? The "do" will execute the "if" test at every iteration, and if it displays "prime" you will have no idea after how many iterations it did that. You may have a composite that runs into the following residue chain: r0, r1, r2, .... rN, 0, -2 (mod Mp), 2, 2, 2, 2, 2, ..., this will still display "prime", for any N smaller than p-2.
Homework: prove this is not possible.

[mattprim]

Originally Posted by mattprim

Lucas-Lehmer test in native Mathematica:

p = 2976221;
s = 4;
Mp = 2^p - 1;
Do[s = Mod[(s*s - 2) , Mp];
If[s == 0, Print["PRIME"]], {i, 0, p - 2}];

puts me in top prime guys fast at least in 60-ties (below GHz*Day).

While PariGP in 70-ties ?: LL(p)={
my(m=Mod(4,1<<p-1));
for(i=3,p,m=m^2-2);
m==0
};
?
? search()={
print("2^2-1");
forprime(p=3,43112609,
if(LL(p), print("2^"p"-1"))
)
};

? search()
2^2-1
2^3-1
2^5-1
2^7-1
2^13-1
2^17-1
2^19-1
2^31-1
2^61-1
2^89-1
2^107-1
2^127-1
2^521-1
2^607-1
2^1279-1
2^2203-1
2^2281-1
2^3217-1
2^4253-1
2^4423-1
2^9689-1
2^9941-1
2^11213-1
2^19937-1
2^21701-1
2^23209-1

[mattprim]

Originally Posted by LaurV

Can you prove this test is correct? I mean, [excluding the fact that you could make it much faster by skipping the initialization (see the syntax of the "do" in Wolfram), and by taking the "if" out of the loop], can you prove that for some composite, you don't hit zero on the way? Can you prove? The "do" will execute the "if" test at every iteration, and if it displays "prime" you will have no idea after how many iterations it did that. You may have a composite that runs into the following residue chain: r0, r1, r2, .... rN, 0, -2 (mod Mp), 2, 2, 2, 2, 2, ..., this will still display "prime", for any N smaller than p-2.
Homework: prove this is not possible.

Here is the proof of LL https://rosettacode.org/wiki/Lucas-Lehmer_test test https://planetmath.org/proofoflucaslehmerprimalitytest . It's kind of tricky because it kicks out of the integer space to irrational roots even if all stays in integer space and is a sort of expanded Little Fermat Theorem https://planetmath.org/inductiveproo...letheoremproof (s_i's with s_(n-2) divisible by M_n if and only if it is the prime are generating quasi-binomial in x=4=2^2 but more complicated). While s_i's are growing giant fast Prime95 for example is using the FFT multiplication of giant integers: http://numbers.computation.free.fr/C...ithms/fft.html using the "electronics signal" fact that the Fourier transform of the convolution f*g is the product of the transforms of two functions.

[kriesel]

Mods: relocate thread to Misc Math?

[Uncwilly]

Originally Posted by kriesel

Mods: relocate thread to Misc Math?

Why? Is the math bad? Is this not about using different software to confirm MP's?

[dcheuk]

Originally Posted by mattprim

Here is the proof of LL https://rosettacode.org/wiki/Lucas-Lehmer_test test https://planetmath.org/proofoflucaslehmerprimalitytest . It's kind of tricky because it kicks out of the integer space to irrational roots even if all stays in integer space and is a sort of expanded Little Fermat Theorem https://planetmath.org/inductiveproo...letheoremproof (s_i's with s_(n-2) divisible by M_n if and only if it is the prime are generating quasi-binomial in x=4=2^2 but more complicated). While s_i's are growing giant fast Prime95 for example is using the FFT multiplication of giant integers: http://numbers.computation.free.fr/C...ithms/fft.html using the "electronics signal" fact that the Fourier transform of the convolution f*g is the product of the transforms of two functions.

I don’t think that’s what LaurV is getting at

Also Fermats littles theorem could be proven using Lagrange in one sentence. And what the heck is integer space to irrational roots you mean extension fields?

It’s not necessary to tell people here how Prime95 works.. what have we done to deserve this?

[mattprim]

Also PariGP LL test for n=41766037 https://www.mersenne.ca/exponent/41766037 on above 1GHz Machine some 12 seconds per step at 240-th step and slowing down. Extimated check time 5800 days above 10 years so native Mathematica above 100 years.

? default(parisize,"64M")
*** Warning: new stack size = 64000000 (61.035 Mbytes).
? LL(p)={
my(m=Mod(4,1<<p-1));print(p);
for(i=3,p,m=m^2-2;print(i));
m==0
};
?
? search()={
print("2^2-1");
p=41766037;
if(LL(p), print("2^"p"-1")
)
};
? search()
2^2-1
41766037
3
4
5
6
.
.
.
239
240