What is millerRabin algorithm
The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm which determines whether a given number is prime, similar to the Fermat primality test and the Solovay–Strassen primality test. It was first discovered by Russian mathematician M. M. Artjuhov in 1967.[1] Gary L. Miller rediscovered it in 1976; Miller’s version of the test is deterministic, but its correctness relies on the unproven extended Riemann hypothesis.[2] Michael O. Rabin modified it to obtain an unconditional probabilistic algorithm in 1980.[3]
PowerMod
What problem it solve
When you want to mod a number which is the product of two big numbers, this algorithm is necessary.
Like a^n mod m.
Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
| // a^n mod m
uint64_t powermod(uint64_t a, uint64_t n, uint64_t m){
uint64_t prod = 1;
// a^11 11 = 8 + 2 + 1 = 01011 // a^11 = a^(1+2+8) = a^1 * a^2 * a^8
// the basic idea is shift the bit. It will more understandable than odd/even number.
// if the last bit is 0, we make a become a^2.
// if the last bit is 1, prod * a.
while(n > 0){
if(n % 2 != 0) // n&1 // n is odd, n&1 is true
prod = prod * a % m;
a = a * a % m; // can not be prod = prod * prod % m
// Because prod * prod will multipul extra a inside.
// For n = 11, n(binary)=01011, (((prod * a)^2 * a)^2^2 * a)^2. this is wrong
// However, prod should be prod * a * a^2 * a^8.
n = n / 2; // compiler n>>=1
}
return prod;
}//O(logn)
|
How Miller-Rabin works
Example code
Here is the sample code
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
67
68
| #include <iostream>
#include <cmath>
using namespace std;
uint64_t powermod(uint64_t a, uint64_t n, uint64_t m){
uint64_t prod = 1;
// a^11 11 = 8 + 2 + 1 = 01011 // a^11 = a^(1+2+8) = a^1 * a^2 * a^8
// the basic idea is shift the bit. It will more understandable than odd/even number.
// if the last bit is 0, we make a become a^2.
// if the last bit is 1, prod * a.
while(n > 0){
if(n % 2 != 0) // n&1 // n is odd, n&1 is true
prod = prod * a % m;
a = a * a % m; // can not be prod = prod * prod % m
// Because prod * prod will multiple extra a inside.
// For n = 11, n(binary)=01011, (((prod * a)^2 * a)^2^2 * a)^2. this is wrong
// However, prod should be prod * a * a^2 * a^8.
n = n / 2; // compiler n>>=1
}
return prod;
}//O(logn)
// https://www.cnblogs.com/Norlan/p/5350243.html
bool MillerRabin(uint64_t p, int k){
if(p < 2) return false;
for (int i = 0; i < k ;i++){
uint64_t a = random(2,p-2);
uint64_t d = p-1;
uint64_t s = 0;
// cout << "MillerRabin a= "<< a << " p="<< p<< endl;
while (d % 2 == 0) {// d & 1 == 0
s++;
d /= 2; //d >> =1
}
// cout <<"d=" << d << " s=" << s << endl;
// d*2^s = p-1 (d is odd number)
// d contains high-order, non-zero bits(stripedd low zreo bits off)
// s = #of bits that were stripp off
uint64_t x = powermod(a,d,p);
cout << "x= " << x << endl;
if(x == 1 || x == p-1)
continue;
// kurger version //
for(int j = 0; j < s-1; j++){
x = x * x % p;
if (x == 1)
return false;
if (x == p-1)
goto nextTry;
}
return false;
nextTry:;
// ************** //
// my version //
// for (int j = 0; j < s-1; ++j)
// {
// x = x * x % p
// if(x == 1)
// return false;
// if(x == p-1)
// break;
// }
// if(x != p-1)
// return false;
// ************** //
}
return true; // probably?
}
|
