Setting up Miller-Rabin in C
source link: https://www.codesd.com/item/setting-up-miller-rabin-in-c.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
Setting up Miller-Rabin in C
I'm trying to implement the Miller-Rabin primality test in C99, but I'm coming across some problems getting it to work. I crafted a small test-set to verify whether or not the implementation works, here's how I'm checking for primes
int main() {
int foo[11] = {0, 1, 2, 3, 4, 7, 28, 73, 125, 991, 1000};
for (int i = 0; i < 11; i++) {
printf("%s; ", isprime(foo[i], 5000) ? "Yes" : "No");
}
return 0;
}
From the numbers listed, the expected output would be
No; No; Yes; Yes; No; Yes; No; Yes; No; Yes; No;
However, as implemented , the output I get is the following:
No; No; Yes; Yes; No; Yes; No; No; No; No; No;
Here's how I wrote the algorithm
int randint (int low, int up){
return rand() % (++up - low)+low;
}
int modpow(int a, int b, int m) {
int c = 1;
while (b) {
if (b & 1) {
c *= a;
}
b >>= 1;
a *= a;
}
return c % m;
}
bool witness(int a, int s, int d, int n) {
int x = modpow(a,d,n);
if(x == 1) return true;
for(int i = 0; i< s-1; i++){
if(x == n-1) return true;
x = modpow(x,2,n);
}
return (x == n- 1);
}
bool isprime(int x, int j) {
if (x == 2) {
return true;
}
if (!(x & 1) || x <= 1) {
return false;
}
int a = 0;
int s = 0;
int d = x - 1;
while (!d&1){
d >>=1;
s+=1;
}
for(int i = 0; i < j; i++){
a = randint(2, x-1);
if(!witness(a,s,d,x)){
return false;
}
}
return true;
}
What am I doing wrong? Why is the test failing for "large" primes, but working for very small ones? How may I fix this?
With Visual Studio 2015 Community edition I found two problems. First the line:
while (!d&1){
needs to be:
while (!(d&1)){
Secondly, as mentioned in the comments your modpow function is overflowing. Try:
int modpow(int a, int d, int m) {
int c = a;
for (int i = 1; i < d; i++)
c = (c*a) % m;
return c % m;
}
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK