在 C 中实现 Miller-Rabin
Implementing Miller-Rabin in C
我正在尝试在 C99 中实现 Miller-Rabin primality test,但在使其工作时遇到了一些问题。我制作了一个小型测试集来验证实现是否有效,下面是我检查素数的方式
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;
}
根据列出的数字,预期输出为
No; No; Yes; Yes; No; Yes; No; Yes; No; Yes; No;
然而,作为实施,我得到的输出如下:
No; No; Yes; Yes; No; Yes; No; No; No; No; No;
下面是我写的算法
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;
}
我做错了什么?为什么测试对 "large" 个素数失败,但对非常小的素数有效?我该如何解决这个问题?
使用 Visual Studio 2015 Community Edition 我发现了两个问题。第一行:
while (!d&1){
需要:
while (!(d&1)){
其次,如评论中所述,您的 modpow 函数溢出。尝试:
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;
}
您的 modpow()
功能有问题。您可能希望对参数和结果使用无符号类型(无论如何,负数 m
是什么意思?)其次,它会溢出,因为它会尝试计算 a^b
,然后再对它取模-m
.您需要减少 a
和 c
。
处理这个问题的最好方法是编写一些测试:
unsigned modpow(unsigned a, unsigned b, unsigned m) {
unsigned c = 1;
while (b) {
if (b & 1) {
c *= a;
}
b >>= 1;
a *= a;
}
return c % m;
}
#include <stdio.h>
unsigned test(unsigned a, unsigned b, unsigned m, unsigned expected)
{
unsigned actual = modpow(a, b, m);
if (actual == expected)
return 0;
printf("modpow(%u, %u, %u) returned %u; expected %u\n",
a, b, m, actual, expected);
return 1;
}
int main()
{
return test(0, 0, 2, 1)
+ test(1005, 16, 100, 25)
;
}
第一个测试通过(但提出了一个问题——m < 2
时你想要什么结果?);第二个失败:
modpow(1005, 16, 100) returned 49; expected 25
让我们修改 modpow()
以减少每一步的结果:
unsigned modpow(unsigned a, unsigned b, unsigned m) {
unsigned c = 1;
while (b) {
if (b & 1) {
c *= a;
c %= m;
}
b >>= 1;
a *= a;
a %= m;
}
return c;
}
现在通过了!我们可以再做一个失败的测试:
int main()
{
return test(0, 0, 2, 1)
+ test(1005, 16, 100, 25)
+ test(100000005, 16, 1000000000, 587890625)
;
}
modpow(100000005, 16, 1000000000) returned 919214529; expected 587890625
现在我们需要使用更大的类型来计算乘法:
unsigned modpow(unsigned a, unsigned b, unsigned m) {
unsigned long long c = 1;
while (b) {
if (b & 1) {
c = (unsigned long long)c * a % m;
}
b >>= 1;
a = (unsigned long long)a * a % m;
}
return c;
}
一旦我们对 modpow()
函数有足够的信心,我们就可以调试算法的其余部分了。
请注意,如果您的整数大小与我的不同,您需要在测试中使用不同的值来复制结果。我选择了以 005
结尾的大数字,因为我们知道最后两位数字是不变的 25
,无论其幂如何。您可能会发现 dc -e '???|p'
有助于生成测试用例(在其标准输入中提供三个参数,它将打印预期值)。
我正在尝试在 C99 中实现 Miller-Rabin primality test,但在使其工作时遇到了一些问题。我制作了一个小型测试集来验证实现是否有效,下面是我检查素数的方式
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;
}
根据列出的数字,预期输出为
No; No; Yes; Yes; No; Yes; No; Yes; No; Yes; No;
然而,作为实施,我得到的输出如下:
No; No; Yes; Yes; No; Yes; No; No; No; No; No;
下面是我写的算法
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;
}
我做错了什么?为什么测试对 "large" 个素数失败,但对非常小的素数有效?我该如何解决这个问题?
使用 Visual Studio 2015 Community Edition 我发现了两个问题。第一行:
while (!d&1){
需要:
while (!(d&1)){
其次,如评论中所述,您的 modpow 函数溢出。尝试:
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;
}
您的 modpow()
功能有问题。您可能希望对参数和结果使用无符号类型(无论如何,负数 m
是什么意思?)其次,它会溢出,因为它会尝试计算 a^b
,然后再对它取模-m
.您需要减少 a
和 c
。
处理这个问题的最好方法是编写一些测试:
unsigned modpow(unsigned a, unsigned b, unsigned m) {
unsigned c = 1;
while (b) {
if (b & 1) {
c *= a;
}
b >>= 1;
a *= a;
}
return c % m;
}
#include <stdio.h>
unsigned test(unsigned a, unsigned b, unsigned m, unsigned expected)
{
unsigned actual = modpow(a, b, m);
if (actual == expected)
return 0;
printf("modpow(%u, %u, %u) returned %u; expected %u\n",
a, b, m, actual, expected);
return 1;
}
int main()
{
return test(0, 0, 2, 1)
+ test(1005, 16, 100, 25)
;
}
第一个测试通过(但提出了一个问题——m < 2
时你想要什么结果?);第二个失败:
modpow(1005, 16, 100) returned 49; expected 25
让我们修改 modpow()
以减少每一步的结果:
unsigned modpow(unsigned a, unsigned b, unsigned m) {
unsigned c = 1;
while (b) {
if (b & 1) {
c *= a;
c %= m;
}
b >>= 1;
a *= a;
a %= m;
}
return c;
}
现在通过了!我们可以再做一个失败的测试:
int main()
{
return test(0, 0, 2, 1)
+ test(1005, 16, 100, 25)
+ test(100000005, 16, 1000000000, 587890625)
;
}
modpow(100000005, 16, 1000000000) returned 919214529; expected 587890625
现在我们需要使用更大的类型来计算乘法:
unsigned modpow(unsigned a, unsigned b, unsigned m) {
unsigned long long c = 1;
while (b) {
if (b & 1) {
c = (unsigned long long)c * a % m;
}
b >>= 1;
a = (unsigned long long)a * a % m;
}
return c;
}
一旦我们对 modpow()
函数有足够的信心,我们就可以调试算法的其余部分了。
请注意,如果您的整数大小与我的不同,您需要在测试中使用不同的值来复制结果。我选择了以 005
结尾的大数字,因为我们知道最后两位数字是不变的 25
,无论其幂如何。您可能会发现 dc -e '???|p'
有助于生成测试用例(在其标准输入中提供三个参数,它将打印预期值)。