dipBRUR

Miller Robin Primality Test

Jan 31st, 2018
95
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 2.39 KB | None | 0 0
  1. // Miller Robin Primality Test
  2.  
  3. // mulMod(a, b, p) : returns (a*b)%p,
  4. // When a, b and p are both long long int
  5. // If we multiply a and b then it will overflow
  6. ll mulMod(ll a, ll b, ll p)
  7. {
  8.     ll x = 0;
  9.     ll y = a%p;    
  10.     while (b)
  11.     {
  12.         if (b & 1)
  13.             x = (x+y)%p;
  14.         y = (y*2)%p;
  15.         b /= 2;
  16.     }
  17.     return x%p;
  18. }
  19.  
  20. // Mod(a, b, p) : return (a^b)%p
  21. // When a, b and p are both long long int
  22. ll Mod(ll a, ll b, ll p)
  23. {
  24.     ll res = 1;
  25.     ll x = a % p;
  26.     while (b)
  27.     {
  28.         if (b & 1)
  29.             res = mulMod(res, x, p);
  30.         x = mulMod(x, x, p);
  31.         b /= 2;
  32.     }
  33.     return res%p;
  34. }
  35.  
  36. // This function is called for all k trials. It returns
  37. // false if n is composite and returns true if n is
  38. // probably prime.
  39. // s is an odd number such that n-1 = 2^d * s, d >= 0
  40. bool millerTest(ll d, ll s, ll n)
  41. {
  42.     // Pick a random number in [2..n-2]
  43.     // Corner cases make sure that n > 4
  44.     ll a = rand() % (n-4) + 2;
  45.    
  46.     // Compute a^s % n
  47.     ll x = Mod(a, s, n);
  48.    
  49.     if (x == 1 || x == n-1)
  50.         return true;
  51.     // Keep squaring x while one of the following doesn't happen
  52.     //   i. s doesn't reach n-1
  53.     //  ii. (x^2)%n is not 1
  54.     // iii. (x^2)%n is not n-1
  55.  
  56.     // i.e. check all r, where 0 < r < d
  57.     //                   where n-1 = 2^d * s
  58.  
  59.     for (int r = 1; r < d; r++)
  60.     {
  61.         x = mulMod(x, x, n);
  62.         if (x == 1)     return false;
  63.         if (x == n-1)   return true;
  64.     }
  65.     // Return composite
  66.     return false;
  67. }
  68.  
  69. // It returns false if 'n' is composite and returns true
  70. // if n is probably prime. 'k' is an input parameter that
  71. // determines accuaracy level. Higher value of 'k' indicates
  72. // more accuaracy.
  73. bool isPrime(ll n, int k = 20)
  74. {
  75.     // Corner cases
  76.     if (n <= 1 || n == 4)  return false;
  77.     if (n <= 3)            return true;
  78.     if (n%2 == 0)          return false;
  79.  
  80.     // find d and s such that, n-1 = 2^d * s, where s = odd, d >= 0
  81.     ll s = n-1, d = 0;
  82.     while (s%2 == 0)
  83.     {
  84.         s /= 2;
  85.         d++;
  86.     }
  87.     // Iterative given number of 'k' times
  88.     for (int i = 0; i < k; i++)
  89.     {
  90.         if (millerTest(d, s, n) == false)
  91.             return false;
  92.     }
  93.     return true;
  94. }
  95.  
  96. int main()
  97. {
  98.     ll n;
  99.     cin >> n;
  100.     if (isPrime(n))    cout << "Prime" << endl;
  101.     else               cout << "Not Prime" << endl;
  102. }
Add Comment
Please, Sign In to add comment