【发布时间】:2015-11-16 08:02:53
【问题描述】:
我目前正在解决 Project Euler 问题 216。 首先,我在 Python 中实现了 Miller-Rabin 素数测试:
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
看起来不错,因为 PE 中 N=10000 匹配的示例。但是python比C++慢,所以我把这段代码翻译成C++。在Problem中,N=5e7,所以我们应该使用long long 64-bit int,对于指数,我们应该使用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;
}
这是确定性代码,因为如果 n
但令人惊讶的是,它打印 2203,而 python 打印 2202。我测试了小数的素数(
我还试图确定是否有任何整数溢出证据,如果isPrime 方法中的 val 得到
为什么会出现这个错误?我在 Windows 上使用了 g++ 4.6.3。
经过一些调试,我发现 C++ 说 2*1939*1939-1 是素数,但实际上不是。
【问题讨论】:
-
尝试确定哪个数字在 C++ 中是素数,但在 Python 中不是素数。然后问 WolframAlpha,它是否真的是素数。也许这会给你一个线索。
-
@user38034 C++ 说 2*1939*1939-1 是素数,但实际上不是。
-
一眼看去,
(!d&1)好像很奇怪,你是说(!(d&1))吗? -
@M.M 我认为你是对的。现在可以了!
标签: c++