用mathematica写的miller rabin算法

mathematica

n = 33940771162697620079218535056615004152380624474841627811237; m = n - 1; s = 0; While[Mod[m, 2] == 0, m = m/2; s = s + 1]; Do[t1 = PowerMod[a, m, n]; k = 0; If[t1 == 1, Print[{a, k, t1, prime}]; Continue[]]; t2 = t1; While[k < s - 1 && t2 != n - 1, k = k + 1; t2 = Mod[t2*t2, n]]; If[t2 == n - 1, Print[{a, k, t2 - n, prime}], Print[{a, composite}]; Break[]], {a, 2, 100}] 其中n是要求判定的奇数，一般很大， 关于循环的范围可以自己选取，这么一个简单的程序我就不 过多地解释了！

mathematica

上面的那个数是一个素数。这个例子给出一个强伪素数，它是前46个素数 的强伪素数。可以看出，最后的结果判定为合数，这个简单的程序给出的 判定结果有些过于详细了！ n = 803837457453639491257079614341942108138837688287558145837488917522\ 2974273765333652186502336163960045457915042023603208766569966760987284\ 0439654082329287387918508691668573282677617710293896977394701670823042\ 8687109997439976544144845341155872450633409279022275296229414984230688\ 1685404326457534018329786111298960644845216191652872597534901; m = n - 1; s = 0; While[Mod[m, 2] == 0, m = m/2; s = s + 1]; Do[t1 = PowerMod[a, m, n]; k = 0; If[t1 == 1, Print[{a, k, t1, prime}]; Continue[]]; t2 = t1; While[k < s - 1 && t2 != n - 1, k = k + 1; t2 = Mod[t2*t2, n]]; If[t2 == n - 1, Print[{a, k, t2 - n, prime}], Print[{a, composite}]; Break[]], {a, 2, 100}] 结果如下 {2,1,-1,prime} {3,0,-1,prime} {4,0,-1,prime} {5,0,1,prime} {6,1,-1,prime} {7,1,-1,prime} {8,1,-1,prime} {9,0,1,prime} {10,1,-1,prime} {11,1,-1,prime} {12,0,1,prime} {13,1,-1,prime} {14,composite} 第14个结果显示该数是合数。 我们也可以分解这个合数。 其中的{a，2，100}，范围可以由用户自己选择。 郭先强的程序也提供了随机的millerrabin测试， 但是我们并不知道其背后产生的基，所以我发一个 输出结果比较详细的mathematica代码

好地方

这个伪素数够强，不过越强的伪素数其实越脆弱，略施小计，只需一眨眼的时间就将它分解为下面两个素数的乘积。 4009582166394996054183064520845468530051881660411325087745062047380032170701196242716223191597219733582163165085358166969145233813917169287527980445796800452592031836601 2004791083197498027091532260422734265025940830205662543872531023690016085350598121358111595798609866791081582542679083484572616906958584643763990222898400226296015918301 而那些无法通过miller－rabin测试的大合数，反而难以分解。是不是数也和人一样，有得必有失，呵呵。

mathematica

其实强伪素数才容易分解了！至少因为目前知道的强伪素数都是 特殊形状的，p-1=k(p-1),n=p*q,目前知道的强伪素数大多是这种 形状的，强伪素数很容易被lucas判定排除掉的，一次不行两次， 两次不行三次！

无心人

我发过一个素性检测Miller-Rabin基选择的帖子 可惜没得到结论 你有兴趣做下 还有个过2，3，5，7的测试的帖子 你可以参考下 另外，GxQ说过，单一测试方法并不是好的 混合测试才好吧 Baillie-PSW应该属于一个混合测试的例子 先过Miller-Rabin强伪素数测试，再Lucas测试

mathematica

原帖由 无心人 于 2008-11-30 20:56 发表 我发过一个素性检测Miller-Rabin基选择的帖子 可惜没得到结论 你有兴趣做下 还有个过2，3，5，7的测试的帖子 你可以参考下 另外，GxQ说过，单一测试方法并不是好的 混合测试才好吧 Baillie-PSW应该属于一个 ...

不知道为什么我看不到以你打头的帖子，虽然你说你发了很多类似之类的。 先吃饭了。我走人了！！！！！！！！！！！！！！！！！！！！！！！ 天已经很晚了！！！！！！！！！！！！！！

无心人

深度隐藏在论坛里呢 你帮论坛晒下旧货 就能找到很多的 好多东西我都忘记自己曾经说过了 呵呵

gxqcn

在算法交流版里，查版主推荐或精华，很容易就找到了：能通过2，3，5，7的检验的合数

mathematica

好地方 发表于 2008-11-30 20:31
这个伪素数够强，不过越强的伪素数其实越脆弱，略施小计，只需一眨眼的时间就将它分解为下面两个素数的乘积 ...

1. Clear["Global`*"];(*Clear all variables*)
   
2. n=803837457453639491257079614341942108138837688287558145837488917522\
   
3. 2974273765333652186502336163960045457915042023603208766569966760987284\
   
4. 0439654082329287387918508691668573282677617710293896977394701670823042\
   
5. 8687109997439976544144845341155872450633409279022275296229414984230688\
   
6. 1685404326457534018329786111298960644845216191652872597534901;
   
7. Solve[p-1==2*(q-1)&&n==p*q&&p>0&&q>0,{p,q}]
   

复制代码

{{p -> 400958216639499605418306452084546853005188166041132508774506204\
7380032170701196242716223191597219733582163165085358166969145233813917\
169287527980445796800452592031836601,
  q -> 200479108319749802709153226042273426502594083020566254387253102\
3690016085350598121358111595798609866791081582542679083484572616906958\
584643763990222898400226296015918301}}

mathematica

mathematica 发表于 2019-1-25 12:59
{{p -> 400958216639499605418306452084546853005188166041132508774506204\
73800321707011962 ...

上面的2是猜出来的,可以一个一个地试

1. N[Sqrt[n/2], 400]

复制代码

看输出结果那么多9,肯定很接近整数了
2.00479108319749802709153226042273426502594083020566254387253102369001\
6085350598121358111595798609866791081582542679083484572616906958584643\
7639902228984002262960159183007499999999999999999999999999999999999999\
9999999999999999999999999999999999999999999999999999999999999999999999\
9999999999999999999999999999999999999999999999999999999999984412340885\
835101213896864153365058221324929248602762249812835*10^168