找回密码
 欢迎注册
楼主: mathematica

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

[复制链接]
 楼主| 发表于 2019-1-26 12:11:47 | 显示全部楼层
本帖最后由 mathematica 于 2019-1-26 12:13 编辑
mathematica 发表于 2008-11-30 18:47
上面的那个数是一个素数。这个例子给出一个强伪素数,它是前46个素数
的强伪素数。可以看出,最后的结果判 ...


观察下面的结果
{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}
提取
{2,1,-1,prime}
{7,1,-1,prime}
这个意味着
\(x^2\equiv-1\pmod{n}\)
\(y^2\equiv-1\pmod{n}\)
那么
\((x*y)^2\equiv1^2\pmod{n}\)
并且x与y互质,
所以可以利用这个同余等式来分解质因数
分别计算
x*y-1与n的最大公约数
x*y+1与n的最大公约数
然后就能得到质因子
4009582166394996054183064520845468530051881660411325087745062047380032\
1707011962427162231915972197335821631650853581669691452338139171692875\
27980445796800452592031836601

代码如下:
  1. Clear["Global`*"];(*Clear all variables*)
  2. n=803837457453639491257079614341942108138837688287558145837488917522\
  3. 2974273765333652186502336163960045457915042023603208766569966760987284\
  4. 0439654082329287387918508691668573282677617710293896977394701670823042\
  5. 8687109997439976544144845341155872450633409279022275296229414984230688\
  6. 1685404326457534018329786111298960644845216191652872597534901;
  7. m=n-1;s=0;
  8. While[Mod[m,2]==0,m=m/2;s=s+1];
  9. kk=PowerMod[2*7,m,n]
  10. ga=GCD[kk-1,n]
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-2-24 09:33:56 | 显示全部楼层
  1. Clear["Global`*"];
  2. MR[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
  3.     s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
  4.     t1=PowerMod[a,m,n];
  5.     If[t1==1,Return[True]];
  6.     k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
  7.     If[t1==n-1,Return[True],Return[False]]
  8. ]
复制代码


这是我重新写的miller rabin算法的代码,又简化了!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-2-24 13:58:54 | 显示全部楼层
mathematica 发表于 2019-2-24 09:33
这是我重新写的miller rabin算法的代码,又简化了!
  1. Clear["Global`*"];
  2. n=2887148238050771212671429597130393991977609459279722700926516024197432\
  3. 3037991527331163289831446392259419778031109293496555784189494417409338\
  4. 0561511397999942154241693397290542371100275104208013496673175515285922\
  5. 6962916775325475044445856101949404200039904432116776619949629539250452\
  6. 6987193290703735640322737012784538991261203092448414947289768854060249\
  7. 76768122077071687938121709811322297802059565867;
  8. (*miller rabin子函数*)
  9. MR[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
  10.     s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
  11.     t1=PowerMod[a,m,n];
  12.     If[t1==1,Return[True]];
  13.     k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
  14.     If[t1==n-1,Return[True],Return[False]]
  15. ]
  16. Do[Print[{k,MR[n,k]}],{k,1,307}]
复制代码

运行结果
  1. {1,True}
  2. {2,True}
  3. {3,True}
  4. {4,True}
  5. {5,True}
  6. {6,True}
  7. {7,True}
  8. {8,True}
  9. {9,True}
  10. {10,True}
  11. {11,True}
  12. {12,True}
  13. {13,True}
  14. {14,True}
  15. {15,True}
  16. {16,True}
  17. {17,True}
  18. {18,True}
  19. {19,True}
  20. {20,True}
  21. {21,True}
  22. {22,True}
  23. {23,True}
  24. {24,True}
  25. {25,True}
  26. {26,True}
  27. {27,True}
  28. {28,True}
  29. {29,True}
  30. {30,True}
  31. {31,True}
  32. {32,True}
  33. {33,True}
  34. {34,True}
  35. {35,True}
  36. {36,True}
  37. {37,True}
  38. {38,True}
  39. {39,True}
  40. {40,True}
  41. {41,True}
  42. {42,True}
  43. {43,True}
  44. {44,True}
  45. {45,True}
  46. {46,True}
  47. {47,True}
  48. {48,True}
  49. {49,True}
  50. {50,True}
  51. {51,True}
  52. {52,True}
  53. {53,True}
  54. {54,True}
  55. {55,True}
  56. {56,True}
  57. {57,True}
  58. {58,True}
  59. {59,True}
  60. {60,True}
  61. {61,True}
  62. {62,True}
  63. {63,True}
  64. {64,True}
  65. {65,True}
  66. {66,True}
  67. {67,True}
  68. {68,True}
  69. {69,True}
  70. {70,True}
  71. {71,True}
  72. {72,True}
  73. {73,True}
  74. {74,True}
  75. {75,True}
  76. {76,True}
  77. {77,True}
  78. {78,True}
  79. {79,True}
  80. {80,True}
  81. {81,True}
  82. {82,True}
  83. {83,True}
  84. {84,True}
  85. {85,True}
  86. {86,True}
  87. {87,True}
  88. {88,True}
  89. {89,True}
  90. {90,True}
  91. {91,True}
  92. {92,True}
  93. {93,True}
  94. {94,True}
  95. {95,True}
  96. {96,True}
  97. {97,True}
  98. {98,True}
  99. {99,True}
  100. {100,True}
  101. {101,True}
  102. {102,True}
  103. {103,True}
  104. {104,True}
  105. {105,True}
  106. {106,True}
  107. {107,True}
  108. {108,True}
  109. {109,True}
  110. {110,True}
  111. {111,True}
  112. {112,True}
  113. {113,True}
  114. {114,True}
  115. {115,True}
  116. {116,True}
  117. {117,True}
  118. {118,True}
  119. {119,True}
  120. {120,True}
  121. {121,True}
  122. {122,True}
  123. {123,True}
  124. {124,True}
  125. {125,True}
  126. {126,True}
  127. {127,True}
  128. {128,True}
  129. {129,True}
  130. {130,True}
  131. {131,True}
  132. {132,True}
  133. {133,True}
  134. {134,True}
  135. {135,True}
  136. {136,True}
  137. {137,True}
  138. {138,True}
  139. {139,True}
  140. {140,True}
  141. {141,True}
  142. {142,True}
  143. {143,True}
  144. {144,True}
  145. {145,True}
  146. {146,True}
  147. {147,True}
  148. {148,True}
  149. {149,True}
  150. {150,True}
  151. {151,True}
  152. {152,True}
  153. {153,True}
  154. {154,True}
  155. {155,True}
  156. {156,True}
  157. {157,True}
  158. {158,True}
  159. {159,True}
  160. {160,True}
  161. {161,True}
  162. {162,True}
  163. {163,True}
  164. {164,True}
  165. {165,True}
  166. {166,True}
  167. {167,True}
  168. {168,True}
  169. {169,True}
  170. {170,True}
  171. {171,True}
  172. {172,True}
  173. {173,True}
  174. {174,True}
  175. {175,True}
  176. {176,True}
  177. {177,True}
  178. {178,True}
  179. {179,True}
  180. {180,True}
  181. {181,True}
  182. {182,True}
  183. {183,True}
  184. {184,True}
  185. {185,True}
  186. {186,True}
  187. {187,True}
  188. {188,True}
  189. {189,True}
  190. {190,True}
  191. {191,True}
  192. {192,True}
  193. {193,True}
  194. {194,True}
  195. {195,True}
  196. {196,True}
  197. {197,True}
  198. {198,True}
  199. {199,True}
  200. {200,True}
  201. {201,True}
  202. {202,True}
  203. {203,True}
  204. {204,True}
  205. {205,True}
  206. {206,True}
  207. {207,True}
  208. {208,True}
  209. {209,True}
  210. {210,True}
  211. {211,True}
  212. {212,True}
  213. {213,True}
  214. {214,True}
  215. {215,True}
  216. {216,True}
  217. {217,True}
  218. {218,True}
  219. {219,True}
  220. {220,True}
  221. {221,True}
  222. {222,True}
  223. {223,True}
  224. {224,True}
  225. {225,True}
  226. {226,True}
  227. {227,True}
  228. {228,True}
  229. {229,True}
  230. {230,True}
  231. {231,True}
  232. {232,True}
  233. {233,True}
  234. {234,True}
  235. {235,True}
  236. {236,True}
  237. {237,True}
  238. {238,True}
  239. {239,True}
  240. {240,True}
  241. {241,True}
  242. {242,True}
  243. {243,True}
  244. {244,True}
  245. {245,True}
  246. {246,True}
  247. {247,True}
  248. {248,True}
  249. {249,True}
  250. {250,True}
  251. {251,True}
  252. {252,True}
  253. {253,True}
  254. {254,True}
  255. {255,True}
  256. {256,True}
  257. {257,True}
  258. {258,True}
  259. {259,True}
  260. {260,True}
  261. {261,True}
  262. {262,True}
  263. {263,True}
  264. {264,True}
  265. {265,True}
  266. {266,True}
  267. {267,True}
  268. {268,True}
  269. {269,True}
  270. {270,True}
  271. {271,True}
  272. {272,True}
  273. {273,True}
  274. {274,True}
  275. {275,True}
  276. {276,True}
  277. {277,True}
  278. {278,True}
  279. {279,True}
  280. {280,True}
  281. {281,True}
  282. {282,True}
  283. {283,True}
  284. {284,True}
  285. {285,True}
  286. {286,True}
  287. {287,True}
  288. {288,True}
  289. {289,True}
  290. {290,True}
  291. {291,True}
  292. {292,True}
  293. {293,True}
  294. {294,True}
  295. {295,True}
  296. {296,True}
  297. {297,True}
  298. {298,True}
  299. {299,True}
  300. {300,True}
  301. {301,True}
  302. {302,True}
  303. {303,True}
  304. {304,True}
  305. {305,True}
  306. {306,True}
  307. {307,False}
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-2-24 14:07:09 | 显示全部楼层
https://bbs.emath.ac.cn/forum.ph ... 734&fromuid=865
用此处的lucas判定办法,就回得到false的答案!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2022-9-19 14:21:20 | 显示全部楼层
好地方 发表于 2008-11-30 20:31
这个伪素数够强,不过越强的伪素数其实越脆弱,略施小计,只需一眨眼的时间就将它分解为下面两个素数的乘积 ...

能说说你是怎么分解的吗?我需要的是方法
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-3-29 15:57 , Processed in 0.055415 second(s), 14 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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