CODEWITHSUNDEEP
c++
Maddy
Maddy

Reputation: 155

Miller Rabin Primality test in C++ bug?

I'm currently solving Project Euler Problem 216. First, I implemented Miller-Rabin primality test in Python:

alist=[2,3,5,7,11,13,17,19,23]
def isPrime(n):
    s=0
    d=n-1
    while not d&1:
        s+=1
        d/=2
    for a in alist:
        if a>=n:continue
        compo=True
        val=pow(a,d,n)
        if val==1 or val==n-1:continue
        for r in range(s-1):
            val=val*val%n
            if val==n-1:
                compo=False
                break
        if(compo):return False
    return True

N=10000
cnt=0
for i in range(2,N+1):
    if isPrime(2*i*i-1):cnt+=1
print cnt

It looks fine because example in PE for N=10000 matches. But python is slower than C++, so I translated this code to C++. In Problem, N=5e7, so we should use long long 64-bit int, and for exponent, we should use 128-bit int.

#include <cstdio>
long long list[9]={2,3,5,7,11,13,17,19,23};
long long exp(int a,long long d,long long n){
    if (d==0)return 1LL;
    if (d&1){
        __int128_t tmp=exp(a,d-1,n);
        return tmp*a%n;
        //return exp(a,d-1,n)*a%n
    }
    __int128_t tmp=exp(a,d/2,n);
    tmp=tmp*tmp%n;
    return tmp;
}
bool isPrime(long long n){
    int s=0;
    long long d=n-1;
    while(!d&1){
        s++;
        d/=2;
    }
    for(int i=0;i<9;i++){
        int a=list[i];
        if(a>=n)continue;
        bool com=true;
        long long val=exp(a,d,n);
        if (val==1||val==n-1)continue;
        for (int r=0;r<s-1;r++){
            __int128_t tmp=val;
            tmp=tmp*tmp%n;
            val=tmp;
            if (val==n-1){
                com=false;
                break;
            }
        }
        if(com)return false;
    }
    return true;
}
int main(){
    long long N=10000;
    int cnt=0;
    for(long long i=2;i<=N;i++){
        if (isPrime(2LL*i*i-1))cnt++;
    }
    printf("%d \n",cnt);
    return 0;
}

This is deterministic code because, if n < 3,825,123,056,546,413,051, it has been shown it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, and 23.

But surprisingly, it prints 2203, while python prints 2202. And I tested primality for small numbers (<1e7), there is no problem. I'm guessing this little difference does not mean my code has serious error such as integer overflow, but some error in using 128-bit integers.

Also i tried to determine is there any integer overflow evidence, I edited my code if val in isPrime method gets <0, then assert, but it didn't assert for N=1e6.

Why did this error occur? I used g++ 4.6.3 on Windows.

After some debugging, I discovered that C++ says 2*1939*1939-1 is a prime, but it is actually not.

Upvotes: 2

Views: 691

Answers (1)

Jakube
Jakube

Reputation: 3565

The error is in the line: 

while(!d&1)

If you look at C++'s Operator Precedence, you see that ! has the priority 3 and & has only the priority 10. That means, that the condition !d&1 gets parsed as (!d)&1. But you actually want !(d&1).

For the future, how can you find such errors? Simple run both codes in parallel and compare some of the values. Since you wrote that the i = 1939 is the bad case. I simple made a breakpoint after the while loop and compared the values of s and d and noticed, that they are different in the Python version and the C++ version. And if you don't want to use a debugger, you can simple insert a line in both codes, that prints the values of s and d.

Upvotes: 6

Related Questions