No Description
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

Miller-Rabin primality test (big).b93 3.6KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134
  1. v works for 23 < n < 3,825,123,056,546,413,051
  2. v <<
  3. 7 ,
  4. 4 > > > >$0v +
  5. * > > $1v 5
  6. > >1+: :2\`#^_:2-!#^_:2%!#^_:9\`#^_:3%!#^_:5%!#^_1 :v 5
  7. v_v#+%3\`\4::_v#+%2\`\3::p13:+1g13<:_v#-+*< 3 .
  8. :>> > $^ 2 1
  9. :| 9 p
  10. 6>\:11p1-0\>:2% !#v_v v ++!!+1-g< # ^5< < :
  11. \# ^\+1\/2< \ >3-#v_$$ 1>31g\ !|> |
  12. `! vp01p03p04 g11p12< >:*11v1 >$1 #$^ >^
  13. \ >120pv v%g04*<v-1\ %g<^ \!!-1:<$0
  14. 5 vg030< v-1\ < >10g^ >\:::*11 g%1-!\^>^
  15. % >$1\> :#v_ $ 21g >:#^_$1-!! ^
  16. + >:!#^_\1+\2v\ ^_^#!%2/\g03p<
  17. >^ ^p02*2g02/ <>:*40g%20g2/:20^
  18. [1, 0] modExp :: a
  19. [2, 0] modExp :: bi
  20. [3, 0] modExp :: b
  21. [4, 0] modExp :: n
  22. [1, 1] witness :: n
  23. [2, 1] witness :: t
  24. [3, 1] isMRPrime :: a
  25. // http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
  26. // ModularExp [a, b, n] -> result
  27. >>> 40p30p10pv v%g04*<
  28. v030p021<v-1\ < >10g^
  29. g >$1\> :#v_ $ >>>
  30. >>:!#^_\1+\2v\ ^_^#!%2/\g03p<
  31. ^p02*2g02/ <>:*40g%20g2/:20^
  32. // Witness [a, n] -> result
  33. >>> :11p1-0\>:2% !#v_vv1-g11\!-1:\!-<
  34. ^\+1\/2< \>+!++3-#v_$$ 1> >>>
  35. vp01p03p04 g11p12< >:*11v%
  36. >120pv v%g04*<v-1\ %g<g
  37. vg030< v-1\ < >10g^ >\:::*11^
  38. >$1\> :#v_ $ 21g >:#^_$1-!! ^
  39. >:!#^_\1+\2v\ ^_^#!%2/\g03p<
  40. ^p02*2g02/ <>:*40g%20g2/:20^
  41. //#############################################################################
  42. bool isMRPrime(int n)
  43. {
  44. if (n <= 1) return false;
  45. if (n == 2) return true;
  46. if (n % 2 == 0) return false;
  47. if (n < 9) return true;
  48. if (n % 3 == 0) return false;
  49. if (n % 5 == 0) return false;
  50. // correct for n < 3,825,123,056,546,413,051
  51. if (Witness(2, n)) return false;
  52. if (Witness(3, n)) return false;
  53. if (Witness(5, n)) return false;
  54. if (Witness(7, n)) return false;
  55. if (Witness(11, n)) return false;
  56. if (Witness(13, n)) return false;
  57. if (Witness(17, n)) return false;
  58. if (Witness(19, n)) return false;
  59. if (Witness(23, n)) return false;
  60. return true;
  61. }
  62. bool Witness(int a, int n)
  63. {
  64. int t = 0;
  65. int u = n - 1;
  66. while (u % 2 == 0)
  67. {
  68. t++;
  69. u /= 2;
  70. }
  71. long xi1 = ModularExp(a, u, n);
  72. for (; t > 0; t--)
  73. {
  74. if (((xi1 * xi1) % n == 1) && (xi1 != 1) && (xi1 != (n - 1))) return true;
  75. xi1 = (xi1 * xi1) % n;
  76. }
  77. if (xi1 != 1) return true;
  78. return false;
  79. }
  80. long ModularExp(int a, int b, int n) // (a^b) % n
  81. {
  82. int k = 0;
  83. int bt = b;
  84. int bi = 1;
  85. while (bt > 0)
  86. {
  87. k++;
  88. bt /= 2;
  89. bi *= 2;
  90. }
  91. long d = 1;
  92. for (k--; k >= 0; k--)
  93. {
  94. d = (d * d) % n;
  95. bi /= 2;
  96. if (((b / bi) % 2) > 0) d = (d * a) % n;
  97. }
  98. return d;
  99. }