合計 昨日 今日

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

[このトピックは解決済みです]

フォーラムルール
フォーラムルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
Name: かびるん
[URL]
Date: 2017年4月17日(月) 00:11
No: 1
(OFFLINE)

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

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

コード[C]: 全て選択
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
#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;
}

Name: かずま
[URL]
Date: 2017年4月17日(月) 00:49
No: 2
(OFFLINE)

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

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

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

Name: かずま
[URL]
Date: 2017年4月17日(月) 00:58
No: 3
(OFFLINE)

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

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

Name: かびるん
[URL]
Date: 2017年4月17日(月) 01:32
No: 4
(OFFLINE)

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

[解決!]

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

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

コード[C++]: 全て選択
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
// 底を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;
}

Name: かびるん
[URL]
Date: 2017年4月17日(月) 01:42
No: 5
(OFFLINE)

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

[解決!]

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


Return to C言語何でも質問掲示板

オンラインデータ

このフォーラムを閲覧中のユーザー: なし & ゲスト[9人]