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.b93 3.3KB

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