数学研发论坛

 找回密码
 欢迎注册
查看: 7682|回复: 15

[原创] 用mathematica写的miller rabin算法

[复制链接]
发表于 2008-11-30 18:38:34 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
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是要求判定的奇数,一般很大,
关于循环的范围可以自己选取,这么一个简单的程序我就不
过多地解释了!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-30 18:47:36 | 显示全部楼层

再给出一个样板

上面的那个数是一个素数。这个例子给出一个强伪素数,它是前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代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-30 20:31:12 | 显示全部楼层
这个伪素数够强,不过越强的伪素数其实越脆弱,略施小计,只需一眨眼的时间就将它分解为下面两个素数的乘积。

4009582166394996054183064520845468530051881660411325087745062047380032170701196242716223191597219733582163165085358166969145233813917169287527980445796800452592031836601

2004791083197498027091532260422734265025940830205662543872531023690016085350598121358111595798609866791081582542679083484572616906958584643763990222898400226296015918301

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

点评

Clear["Global`*"];(*Clear all variables*) n=803837457453639491257079614341942108138837688287558145837488917522\ 2974273765333652186502336163960045457915042023603208766569966760987284\ 043965408   发表于 2019-1-25 12:59
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-30 20:45:17 | 显示全部楼层

你的分解算不了什么

其实强伪素数才容易分解了!至少因为目前知道的强伪素数都是
特殊形状的,p-1=k(p-1),n=p*q,目前知道的强伪素数大多是这种
形状的,强伪素数很容易被lucas判定排除掉的,一次不行两次,
两次不行三次!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-30 20:56:25 | 显示全部楼层
我发过一个素性检测Miller-Rabin基选择的帖子
可惜没得到结论
你有兴趣做下
还有个过2,3,5,7的测试的帖子
你可以参考下

另外,GxQ说过,单一测试方法并不是好的
混合测试才好吧
Baillie-PSW应该属于一个混合测试的例子
先过Miller-Rabin强伪素数测试,再Lucas测试
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-11-30 21:06:06 | 显示全部楼层
原帖由 无心人 于 2008-11-30 20:56 发表
我发过一个素性检测Miller-Rabin基选择的帖子
可惜没得到结论
你有兴趣做下
还有个过2,3,5,7的测试的帖子
你可以参考下

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


不知道为什么我看不到以你打头的帖子,虽然你说你发了很多类似之类的。
先吃饭了。我走人了!!!!!!!!!!!!!!!!!!!!!!!
天已经很晚了!!!!!!!!!!!!!!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-11-30 21:20:14 | 显示全部楼层


深度隐藏在论坛里呢
你帮论坛晒下旧货
就能找到很多的
好多东西我都忘记自己曾经说过了
呵呵
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-12-1 07:39:50 | 显示全部楼层

回复 6# mathematica 的帖子

在算法交流版里,查版主推荐或精华,很容易就找到了:能通过2,3,5,7的检验的合数
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-1-25 12:59:49 | 显示全部楼层
好地方 发表于 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}}
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-1-25 13:03:22 | 显示全部楼层
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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2022-10-3 06:58 , Processed in 0.074480 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表