#11
by かずま » 8年前
clear さんが書きました:
algorithm miller-rabin-test(n,l)
begin
Compute odd number u and positive integer k such that n = u*2^k;
for i ← 1 to l do
if M-R-T(n, u, k) = NO then return(NO);
return(YES);
end
function M-R-T(n, u, k)
begin
Pick a ∈ {2, 3, ..., n-2} at random;
b ← a^u mod n;
if b ∈ {1, n-1} then return(YES);
for i ← to k-1 do begin
b ← b^2 mod n;
if b = n-1 then return(YES);
if b = 1 then return(YES);
end;
return(NO);
end
n = u*2
k だったら、n は 2の倍数なので素数ではありません。
Wikipedia を見たら、n-1 = u*2
k でした。
if b ∈ {1, n-1} then return(YES);
と
if b = n-1 then return(YES);
if b = 1 then return(YES);
とは、同じ意味ですよね。C で書くと、
if (b == 1 || b == n-1) return YES;
多倍長計算を実装するのは面倒なので unsigned long long で、
9桁までなら、アルゴリズムの検証ができるかやってみました。
コード:
#include <stdio.h> // printf
#include <stdlib.h> // srand, rand
#include <time.h> // time
#define YES 1
#define NO 0
typedef unsigned long long TYPE;
TYPE mod_pow(TYPE x, TYPE y, TYPE n)
{
TYPE z = x;
while (--y > 0) z = z * x % n;
return z;
}
TYPE MRT(TYPE n, TYPE u, TYPE k)
{
if (n <= 3) return YES;
TYPE a = rand() % (n - 3) + 2;
TYPE b = mod_pow(a, u, n);
//printf("n=%llu, u=%llu, k=%llu, a=%llu, b=%llu\n", n, u, k, a, b);
if (b == 1 || b == n-1) return YES;
for (TYPE i = 0; i < k; i++) {
b = mod_pow(b, 2, n);
if (b == 1 || b == n-1) return YES;
}
return NO;
}
TYPE miller_rabin_test(TYPE n, TYPE l)
{
TYPE u = n - 1, k = 0;
while ((u & 1) == 0) u >>= 1, k++;
for (TYPE i = 0; i < l; i++)
if (MRT(n, u, k) == NO) return NO;
return YES;
}
int main(void)
{
srand(time(0));
printf("start n: ");
TYPE n;
if (scanf("%llu", &n) != 1) return 1;
n |= 1; // 奇数にする
printf("times k: ");
int k;
if (scanf("%d", &k) != 1) return 1;
k /= 2; // 奇数だけなので回数は半分
for (int i = 0; i < k; i++, n += 2)
if (miller_rabin_test(n, 5) == YES) {
printf(" %llu", n);
fflush(stdout);
}
printf("\n");
return 0;
}
実行結果
コード:
start n: 99999900
times k: 100
99999931 99999941 99999959 99999971 99999989
すごく時間がかかるので、多倍長計算でも 100桁やるのは大変だと思います。
[quote="clear" id=3,19259,145823]
[b]algorithm[/b] miller-rabin-test(n,l)
[b]begin[/b]
Compute odd number [i]u[/i] and positive integer [i]k[/i] such that [i]n = u*2^k[/i];
[b]for[/b] [i]i[/i] ← 1 [b]to[/b] [i]l[/i] [b]do[/b]
if M-R-T([i]n, u, k[/i]) = NO [b]then return[/b](NO);
[b]return[/b](YES);
[b]end[/b]
[b]function[/b] M-R-T([i]n, u, k[/i])
[b]begin[/b]
Pick [i]a[/i] ∈ {2, 3, ..., n-2} at random;
[i]b[/i] ← [i]a^u[/i] mod [i]n[/i];
[b]if[/b] [i]b[/i] ∈ {1, n-1} [b]then return[/b](YES);
[b]for[/b] [i]i[/i] ← [b]to[/b] k-1 [b]do begin[/b]
[i]b[/i] ← [i]b^2[/i] mod [i]n[/i];
[b]if[/b] [i]b[/i] = n-1 [b]then return[/b](YES);
[b]if[/b] [i]b[/i] = 1 [b]then return[/b](YES);
[b]end[/b];
[b]return[/b](NO);
[b]end[/b][/quote]
n = u*2[sup]k[/sup] だったら、n は 2の倍数なので素数ではありません。
Wikipedia を見たら、n-1 = u*2[sup]k[/sup] でした。
if b ∈ {1, n-1} then return(YES);
と
if b = n-1 then return(YES);
if b = 1 then return(YES);
とは、同じ意味ですよね。C で書くと、
if (b == 1 || b == n-1) return YES;
多倍長計算を実装するのは面倒なので unsigned long long で、
9桁までなら、アルゴリズムの検証ができるかやってみました。
[code=c]
#include <stdio.h> // printf
#include <stdlib.h> // srand, rand
#include <time.h> // time
#define YES 1
#define NO 0
typedef unsigned long long TYPE;
TYPE mod_pow(TYPE x, TYPE y, TYPE n)
{
TYPE z = x;
while (--y > 0) z = z * x % n;
return z;
}
TYPE MRT(TYPE n, TYPE u, TYPE k)
{
if (n <= 3) return YES;
TYPE a = rand() % (n - 3) + 2;
TYPE b = mod_pow(a, u, n);
//printf("n=%llu, u=%llu, k=%llu, a=%llu, b=%llu\n", n, u, k, a, b);
if (b == 1 || b == n-1) return YES;
for (TYPE i = 0; i < k; i++) {
b = mod_pow(b, 2, n);
if (b == 1 || b == n-1) return YES;
}
return NO;
}
TYPE miller_rabin_test(TYPE n, TYPE l)
{
TYPE u = n - 1, k = 0;
while ((u & 1) == 0) u >>= 1, k++;
for (TYPE i = 0; i < l; i++)
if (MRT(n, u, k) == NO) return NO;
return YES;
}
int main(void)
{
srand(time(0));
printf("start n: ");
TYPE n;
if (scanf("%llu", &n) != 1) return 1;
n |= 1; // 奇数にする
printf("times k: ");
int k;
if (scanf("%d", &k) != 1) return 1;
k /= 2; // 奇数だけなので回数は半分
for (int i = 0; i < k; i++, n += 2)
if (miller_rabin_test(n, 5) == YES) {
printf(" %llu", n);
fflush(stdout);
}
printf("\n");
return 0;
}
[/code]
実行結果
[code=text]
start n: 99999900
times k: 100
99999931 99999941 99999959 99999971 99999989
[/code]
すごく時間がかかるので、多倍長計算でも 100桁やるのは大変だと思います。