有谁能看懂这段梅森数的lucas-lehmer判别法的代码？

mathematica

1. LucasLehmer[n_] :=
   
2. Module[{Mp, s, k}, Mp = 2^n - 1; s = 4;
   
3.   Do[s = Mod[s^2 - 2, Mp], {k, 1, n - 2}];
   
4.   If[s == 0, Return[True], Return[False]]]
   
5. 
   
6. LucasLehmer2[p_] :=
   
7. Module[{Mp, s, k}, Mp = 2^p - 1; s = 4;
   
8.   Do[s = s^2 - 2; If[s < 0, s += Mp];
   
9.    s = BitAnd[s, Mp] + BitShiftRight[s, p];
   
10.    If[s >= Mp, s -= Mp], {k, 1, p - 2}];
    
11.   If[s == 0, Return[True], Return[False]]]

复制代码

https://baike.baidu.com/item/%E5 ... /3757240?fr=aladdin
卢卡斯-莱默检验法 ，用来判别梅森数是不是素数的算法

s = Mod[s^2 - 2, Mp]
这行代码被替换成了下面的
s = s^2 - 2;
If[s < 0, s += Mp];
s = BitAnd[s, Mp] + BitShiftRight[s, p];
If[s >= Mp, s -= Mp]
这代码替换了一下，速度快了
执行
Timing[LucasLehmer2[86243]]
运行结果
{10.9825, True}
执行
Timing[LucasLehmer[86243]]
运行结果
{49.4367, True}
仅仅换了一下代码，运行效率有五倍之差！！！！！！！！！！
效率低的代码是我写的，效率高的代码是别人写的！
我最近研究mathematica软件与prime95的时候，对比发现判别梅森数方面，mathematica比prime95慢了很多！
我还怀疑prime95，真的是越来越发现自己很无知！
这段取模的置换的代码，是不是只对梅森数有效？

mathematica

1. LucasLehmer[n_] := Module[{Mp, s, k}, Mp = 2^n - 1; s = 4;
   
2.   Do[s = Mod[s^2 - 2, Mp]; Print[s], {k, 1, n - 2}];
   
3.   If[s == 0, Return[True], Return[False]]]
   
4. 
   
5. LucasLehmer2[p_] := Module[{Mp, s, k}, Mp = 2^p - 1; s = 4;
   
6.   Do[s = s^2 - 2; If[s < 0, s += Mp];
   
7.    s = BitAnd[s, Mp] + BitShiftRight[s, p];
   
8.    If[s >= Mp, s -= Mp]; Print[s], {k, 1, p - 2}];
   
9.   If[s == 0, Return[True], Return[False]]]

复制代码

我当然也怀疑代码的正确性，就把代码的中间结果也打印出来，每个函数都增加了

1. Print[s]

复制代码

这句，
结果用LucasLehmer[11]与LucasLehmer2[11]比较中间结果，发现结果完全一样！

问题来了，是不是这种算法只有对梅森数有效果？这段代码的原理又是什么？
真把我打蒙了！

.·.·.

mathematica 发表于 2019-3-4 12:15
我当然也怀疑代码的正确性，就把代码的中间结果也打印出来，每个函数都增加了这句，
结果用LucasLehmer[ ...

取模运算比位运算慢是很显然的事情
你写的是一般的取模运算
人家的代码把特定的取模运算（模Mp）简化成了位运算
所以人家的算法快
这是很正常的事情

原理差不多是
a=2^p*(a-(a mod 2^p))/2^p+(a mod 2^p)
两边mod2^p-1
a mod 2^p-1=(a-(a mod 2^p))/2^p+(a mod 2^p)
/2^p被bitshift替代,mod 2^p被bitand(#,2^p-1)&替代