ミラー-ラビン素数判定法のC言語での実装

[かびるん]

こんにちは。
ミラー-ラビン素数判定法を用いて決定的素数判定をするため、Cでコードを書きました。
しかし確認してみると、50000を超えたあたりから素数も通らなくなりました。
10000以下はすべて正しく判定出来ているようなのですが…
どこかおかしいところはないか、ご教示願います。

コード:

#include<stdio.h>
#include<math.h>
// 冪剰余(a^b %d の計算)
int modpow(int a, long long b, int d) {
  if (b == 0)     return 1;
  if (b % 2 == 0) return modpow(a*a % d, b/2, d);
  else            return a*modpow(a, b-1, d) % d;
}

bool Miller_Rabin_Test(int n){
    int prime_list[4] = {2, 3, 61};

    if(n == 2 || n == 3 || n == 61){  // 2,3,61は素数
        return true;
    }
    if(n == 1 || n%2 == 0){  // 2以外の偶数と1は素数でない
        return false;
    }
    // n = 2^s*d + 1 (dは奇数)の形に変形
    int d = n-1;
    int s = 0;
    while(d%2 == 0){  // dが奇数になるまで変形
        d = d/2;
        s++;
    }
    for(int i=0; i<3; i++){  // 判定する底は3つ
        bool isProbPrime = false;
        int a = prime_list[i];
        int r = modpow(a, d, n);
        if(r == 1 || r == n-1){
            isProbPrime = true;
        }
        for(int j=0; j<s; j++){
            r = modpow(a, pow(2,j)*d, n);
            if(r == n-1){
                isProbPrime = true;
            }
        }
        if(!isProbPrime){
            return false;
        }
    }
    return true;
}

int main(){
    int n = 0;
    scanf_s("%d",&n);
    if(Miller_Rabin_Test(n)){
        printf("%d is prime\n",n);
    }else{
        printf("%d is not prime\n",n);
    }
    return 0;
}

[かずま]

6行目を何度も通ると、a*a のせいで、a の値が
どんどん大きくなって、int の最大値を超えています。

int modpow(long long a, long long b, int d) にしてみては？

[かずま]

もう少しちゃんと説明すると、6行目の a*a % d ですが、入力である d が
50000 を超えると、余りも 50000 になることがあって、それを新しい a と
すると、次の a*a が 25億になってしまいます。
これは int の最大値 2147483647を超えるということです。
だから、long long a にしましょう。

[かびるん]

なるほど。ありがとうございます。
unsigned long long型を双方に適用して、型の最大値の平方根 4294967295 までの素数判定が決定的になることがわかりました。
(とりあえず1億まで数え上げさせたら、個数が一致しました)

それと、判定の底が間違ってたみたいなので直しました。

コード:

// 底を2,7,61にしてMiller–Rabin素数判定法を用いると、
// 2^32までの素数判定は決定的であることが保証されている.
#include<stdio.h>
#include<math.h>
#define ullint unsigned long long int
// 冪剰余(a^b %d の計算)
int modpow(ullint a, ullint b, int d) {
  if (b == 0)     return 1;
  if (b % 2 == 0) return modpow(a*a % d, b/2, d);
  else            return a*modpow(a, b-1, d) % d;
}

bool Miller_Rabin_Test(int n){
    int prime_list[4] = {2, 7, 61};

    if(n == 2 || n == 7 || n == 61){  // 2,7,61は素数
        return true;
    }
    if(n == 1 || n%2 == 0){  // 2以外の偶数と1は素数でない
        return false;
    }
    // n = 2^s*d + 1 (dは奇数)の形に変形
    int d = n-1;
    int s = 0;
    while(d%2 == 0){  // dが奇数になるまで変形
        d = d/2;
        s++;
    }
    for(int i=0; i<3; i++){  // 判定する底は3つ
        bool isProbPrime = false;
        int a = prime_list[i];
        int r = modpow(a, d, n);
        if(r == 1 || r == n-1){
            isProbPrime = true;
        }
        for(int j=0; j<s; j++){
            r = modpow(a, pow(2,j)*d, n);
            if(r == n-1){
                isProbPrime = true;
            }
        }
        if(!isProbPrime){
            return false;
        }
    }
    return true;
}
int main(){
    int n = 0;
    scanf_s("%d",&n);
    if(Miller_Rabin_Test(n)){
        printf("%d is prime\n",n);
    }else{
        printf("%d is not prime\n",n);
    }
    return 0;
}

[かびるん]

あっ、breakとcontinue挟み忘れてました…