费尔马二平方素数

[medie2005]

题目：

除了2这个特别的素数外，所有的素数都可以分成两类：第一类是被4除余
1的素数，如5，13，17，29，37，41；第二类是被4除余3的素数，如3，7，
11，19，23，31。第一类素数都能表示成两个整数的平方和(第二类不能)，
例如：5=1*1+2*2、13=2*2+3*3、17=1*1+4*4、29=2*2+5*5...这就是著名的
费尔马“二平方”定理。有趣的是：上述等式右侧的数有的又恰恰是两个素
数的平方，如13、29，我们就把这样的素数叫作费尔马“二平方”素数，即
是如果一个素数能够表示成两个素数的平方和的形式，例如：F=X*X+Y*Y
(1)，其中F、X、Y都是素数，它就是费尔马“二平方”素数。

请求出并打印40亿以内的所有的费尔马“二平方”素数。

我已经写了一个，要6秒左右，麻烦高手们想一个算法使运行时间小于1秒。

[kingbo2006]

能把你程序贴出来看看吗

[mathe]

也就是说p^2+4=q
可以先找出所有小于63245=sqrt(4*10^9)的素数。
对于每个这样的奇素数p，计算p^2+4,然后我们需要判断p^2+4是否是素数。
如果这个数字小于63245,那么我们可以利用现成的表格判断，不然，可以用所有小于p的型如
4K＋1的素数去除，如果都不能整除，那么p^2+4就是素数，输出

[mathe]

如果发现这样计算速度太慢，我们可以使用素数探测方法来加速。
比如对于q=p^2+4,q比较大时，
我们首先计算2^(q-1) (mod q),如果这个数字不是1，那么q肯定不是素数。
如果2^(q-1)(mod q)=1,那么我们无法判断，然后再使用前面较慢的方法。

[mathe]

计算a^b (mod q)的代码如下：
power_mod(int a, int b, int q){
int mul=a%q, result=1;
while(b){
if(b&1){
result*=mul;result%=q;
}
mul*=mul; mul%=q;
b>>=1;
}
return result;
}

[mathe]

需要注意的是上面乘法可能会溢出，所以要使用64位的整数，在Windows里面用__int64,在Linux里面用long long.

[kingbo2006]

质数除了2以外都是奇数，这样你可以缩小一半的搜索量啊

[mmmcd]

用上述想法实现的程序

找到368个

#include <iostream.h>
#include <fstream.h>
const int N=63246;
bool nprime[N];
int prime[N],pl;
ofstream fout("out.txt");
void init()
{
int i,j;
for(i=2;i<N;i++){
if(nprime[i])
continue;
if(i%4==1)
prime[pl++]=i;
for(j=i+i;j<N;j+=i){
nprime[j]=true;
}
}
}

void output(__int64 p){
if(p>=10)
output(p/10);
cout<<(int)(p%10);
fout<<(int)(p%10);
}

int main()
{
int i,j,d;
__int64 p;
init();
for(i=0;i<pl;i++){
p=d=prime[i];
p=p*p+4;
if(d<=251){
if(!nprime[p]){
cout<<(int)p<<'='<<"2*2+"<<d<<'*'<<d<<endl;
fout<<(int)p<<'='<<"2*2+"<<d<<'*'<<d<<endl;
}
}else{
for(j=0;j<i;j++){
if(p%prime[j]==0)
break;
}
if(j>=i){
output(p);
cout<<'='<<"2*2+"<<d<<'*'<<d<<endl;
fout<<'='<<"2*2+"<<d<<'*'<<d<<endl;
}
}
}
return 0;
}

[mmmcd]

你的方法：for(i=53;i<4000000000;i+=120)//对于<4000000000的每个整数i，判断sqrt(i-4)是否素数；这个循环效率太低了

建议这样：对于<63246的每个素数i，判断i*i+4是否素数；

我们判断素数的方法一样，
但素数个数远远小于整数的个数，时间就相差在这里

[mmmcd]

(4000000000－53)/120=33333332.9

而小于63245的素数只有 6337 个

[kingbo2006]

对于MMMCD的法，我也赞同。判断10万以内的的素数所费的时间要比10万以上的素数要费时。不信的话你可以做一下试验。大约要多费0.001~0.01S的时间，这样在比较的判断中就会产生很大时间。

[mmmcd]

改了一下，一开始不应把模4余3的素数去掉。
找到718个。（运行楼上的程序似乎也是找718个哦）

#include <iostream.h>
#include <fstream.h>
#include <windows.h>
const int N=63246;
bool nprime[N];
int prime[N],pl;
ofstream fout("out.txt");
void init()
{
int i,j;
for(i=2;i<N;i++){
if(nprime[i])
continue;
//if(i%4==1)
prime[pl++]=i;
for(j=i+i;j<N;j+=i){
nprime[j]=true;
}
}
}

void output(__int64 p){
if(p>=10)
output(p/10);
//cout<<(int)(p%10);
fout<<(int)(p%10);
}

int main()
{
int i,j,d,c=0;
int time1,time2;
__int64 p;
time1=GetTickCount();
init();
time2=GetTickCount();
cout<<"筛选素数法，用时:"<<(time2-time1)/1000.0<<"s"<<endl;
for(i=0;i<pl;i++){
p=d=prime[i];
p=p*p+4;
if(d<=251){
if(!nprime[p]){
//cout<<(int)p<<'='<<"2*2+"<<d<<'*'<<d<<endl;
fout<<(int)p<<'='<<"2*2+"<<d<<'*'<<d<<endl;
c++;
}
}else{
for(j=0;j<i;j++){
if(p%prime[j]==0)
break;
}
if(j>=i){
output(p);
//cout<<'='<<"2*2+"<<d<<'*'<<d<<endl;
fout<<'='<<"2*2+"<<d<<'*'<<d<<endl;
c++;
}
}
}
time2=GetTickCount();
cout<<"总用时:"<<(time2-time1)/1000.0<<"s"<<endl;
cout<<"找到"<<c<<"个"<<endl;
return 0;
}


我的机子比较高级一点，输出如下

筛选素数法，用时:0s
总用时:0.203s
找到718个

[kingbo2006]

祝贺我们从不同的角度，解决了这样一个问题。以后多多合作！谢谢！

[medie2005]

高手们，我的程序错在哪里？麻烦大家看看啊，解决后立即揭贴。

如下：

//计算4000000000以内的费尔马素数，用时6-7秒
#include<iostream>
#include<cmath>
#include<windows.h>
using namespace std;
int primelist[ 6300 ];
bool isprime(int n)
{
unsigned int i=0;
unsigned int sn=(int)sqrt(n);
for(;primelist[i]<sn+1;++i)
if(n%primelist[ i ]==0)
return false;
return true;
}
void bornprimelist()
{
primelist[ 0 ]=3;
unsigned int i,count=0;
int step=4;
for(i=5;i<62500;i+=step)
{
if(isprime(i))
count++,primelist[ count ]=i;
step=6-step;
}
}
int main()
{
unsigned int i;
cout<<"13==2*2+3*3"<<endl;
cout<<"29==2*2+5*5"<<endl;
int count=2;
int time=GetTickCount();
bornprimelist();
for(i=53;i<4000000000;i+=120)
{
double dt=sqrt(i-4);
if(dt==(int)dt)
if(isprime((int)dt))
if(isprime(i))
{
cout<<i<<"==2*2+"<<dt<<'*'<<dt<<endl;
count++;
}
}
cout<<"total :"<<count<<endl;
cout<<"Use time :"<<GetTickCount()-time<<"ms"<<endl;
return 1;
}

[kingbo2006]

我的程序设计如下：
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <fstream.h>
#include <windows.h>
using namespace std;
const int IsSuXu(const int x)//判断x是不是素数；
{
unsigned __int64 result;
unsigned __int64 temp=0,k=3;
int mod=((x%2)!=0);
unsigned __int64 _End= (int)sqrt((double)x);
while(k<=_End&&(x%k)!=0&&mod)
{
if((k%3)==0&&(k!=3))
{
k+=2;
temp++;
}
else
{
k+=2;
temp++;
}
}
if(k>_End) result=1;
else result=0;
return result;
}
int main(int argc, char *argv[])
{
double time1,time2;
unsigned __int64 i=3,count1=1,temp=0,count2=0;
int SuXu[6337];
ofstream fout("out.txt");
time1=GetTickCount();
SuXu[0]=2;
while(i<=63245)
{
if((i%3)==0&&(i!=3))
{
i+=2;
temp++;
}
else
{
if(IsSuXu(i)){SuXu[count1]=i;count1++;}
i+=2;
temp++;
}
}
time2=GetTickCount();
cout<<(time2-time1)/1000.0<<" s"<<endl;
time1=GetTickCount();
for(i=1;i<6337;i++)
{
temp=SuXu[i];
int j=0;
while(j<=i&&(temp*temp+4)%SuXu[j]!=0)
j++;
if(j>i)
{
//cout<<temp<<"*"<<temp<<"+4="<<temp*temp+4<<endl;
fout<<temp<<"*"<<temp<<"+4="<<temp*temp+4<<endl;
count2++;
}
}

//cout<<Two_Power(100)<<endl;;
cout<<IsSuXu(63241)<<endl;
time2=GetTickCount();
cout<<"count1="<<count1<<endl;
cout<<"count2="<<count2<<endl;
cout<<(time2-time1)/1000.0<<" s"<<endl;
system("PAUSE");
return 0;
}
程序的结果如下：
费尔马二平方素数共计716个：
用时0.3s左右。

[medie2005]

谢谢！！