-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathmiller rabin preimality test.cpp
67 lines (58 loc) · 1.43 KB
/
miller rabin preimality test.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
/// ***Miller-Rabin Primality test [120]
/** usigned long long is not working sometimes
* left shift == mult by 2 right shift == divide bye 2
* multiplying two numbers (a*b)%c avoiding overflow
*/
ll mulmod(ll a, ll b, ll mod){
ll x = 0,y = a % mod;
while (b > 0){
if (b&1) x = (x + y) % mod;
y = (y<<1) % mod;
b >>= 1;
}
return x % mod;
}
//Bigmod
ll modulo(ll n, ll r, ll mod){
ll x = 1;
while (r > 0){
if (r&1) x = mulmod(x,n,mod);
n = mulmod(n,n,mod);
r >>= 1;
}
return x % mod;
}
/// higher value of "it" ensure higher percision [recomendation 7]
bool MillerRabin(ll p,int it)
{
if (p < 2) return 0;
if (p != 2 && p % 2==0) return 0;
ll i,a,tmp,mod,s=p-1;
while(s%2==0){
s>>=1;
}
for(i=0;i<it;i++){
a = rand()%(p-1)+1;
tmp = s;
mod = modulo(a, tmp, p);
if(mod==1 || mod == p-1)continue;
while (tmp!=p-1&&mod!=1&&mod!=p-1){
mod = mulmod(mod, mod, p);
tmp <<= 1;
}
if(mod!=p-1) return 0;
}
return 1;
}
/** this body of miller rabin giving correct ans upto 10^9
* with lower time complexity
*/
for(i=1;i<=it;i++){
a = rand() % (p - 1) + 1, tmp = s;
mod = mod_pow(a, tmp, p);
while(tmp != p - 1 && mod != 1 && mod != p - 1) {
mod = mod_mul(mod, mod, p);
tmp *= 2;
}
if (mod != p - 1 && tmp % 2 == 0) return false;
}