用mathematica写的miller rabin算法

SAY

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是要求判定的奇数，一般很大，
关于循环的范围可以自己选取，这么一个简单的程序我就不
过多地解释了！

evuojad

上面的那个数是一个素数。这个例子给出一个强伪素数，它是前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代码

gocehuzana

这个伪素数够强，不过越强的伪素数其实越脆弱，略施小计，只需一眨眼的时间就将它分解为下面两个素数的乘积。

4009582166394996054183064520845468530051881660411325087745062047380032170701196242716223191597219733582163165085358166969145233813917169287527980445796800452592031836601

2004791083197498027091532260422734265025940830205662543872531023690016085350598121358111595798609866791081582542679083484572616906958584643763990222898400226296015918301

而那些无法通过miller－rabin测试的大合数，反而难以分解。是不是数也和人一样，有得必有失，呵呵。

Rickykek

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

ogemama

我发过一个素性检测Miller-Rabin基选择的帖子
可惜没得到结论
你有兴趣做下
还有个过2，3，5，7的测试的帖子
你可以参考下

另外，GxQ说过，单一测试方法并不是好的
混合测试才好吧
Baillie-PSW应该属于一个混合测试的例子
先过Miller-Rabin强伪素数测试，再Lucas测试

ugejoleyaz

原帖由 无心人 于 2008-11-30 20:56 发表
我发过一个素性检测Miller-Rabin基选择的帖子
可惜没得到结论
你有兴趣做下
还有个过2，3，5，7的测试的帖子
你可以参考下

另外，GxQ说过，单一测试方法并不是好的
混合测试才好吧
Baillie-PSW应该属于一个 ...

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

uajoegecp

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

BristSes

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

mirakhatun

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

复制代码

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

saluozu

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

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

复制代码

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