#include <stdio.h>
#define PC 1000000
#define N 2293
char isprime;
void init()
{
long long i;
isprime=isprime=0;
for(i=2;i<PC;i++)isprime=1;
for(i=2;i*i<PC;i++){
if(isprime){
int j;
for(j=i*i;j<PC;j+=i)isprime=0;
}
}
}
void solve(int p)
{
int i,r,t=N%p,f=0;
if(t==0)return;
i=0,r=1;
do{
if(r==t){
f=1;
printf("%d,",i);
}
r*=2;if(r>=p)r-=p;i++;
if(r==1){
if(f==1) printf("%d\n",i);
break;
}
}while(1);
}
int main()
{
int i;
init();
for(i=3;i<PC;i++){
if(isprime){
solve(i);
}
}
return 0;
}
第二个代码把上面代码输出简单处理一下，产生第一步筛选结果
#include <stdio.h>
#define IN "pp.out"
#define BUFC 100000000
#define MACC (8*BUFC)
char buf;
#define BYTE(x)((x)>>3)
#define INDEX(x)((x)&7)
#define SET(x) (buf|=1<<INDEX(x))
#define GET(x) (buf&(1<<INDEX(x)))
int getd(int m){
if(m%12==0)return 12;
if(m%6==0)return 6;
if(m%4==0)return 4;
if(m%3==0)return 3;
if(m%2==0)return 2;
return 1;
}
int main()
{
int d,m,h,s;
FILE *in=fopen(IN,"r");
if(in==NULL){
fprintf(stderr, "Fail to open input file\n");
return -1;
}
while(fscanf(in,"%d,%d",&d,&m)==2){
int md=getd(m);
if(md>1&&d%md!=1)continue;//skip data non compatible with 1(mod 12)
m/=md;d%=m;
if(m==1){
fprintf(stderr,"Error, all data is filtered away\n");
return -1;
}
//find s so that d+sm=1(mod 12)
for(s=0;s<12;s++){
if((d+s*m)%12==1)break;
}
if(s==12){
fprintf(stderr,"Logic error\n");
return -1;
}
for(h=(d+s*m-1)/12;h<MACC;h+=m){
SET(h);
}
}
fclose(in);
for(h=0;h<MACC;h++){
if(GET(h)==0){
printf("%d\n",h*12+1);
}
}
return 0;
}
最后是伪素数判断#include <gmpxx.h>
#include <stdio.h>
#include <time.h>
#define IN "ppc.out"
#define N2293
int gc=0;
time_t sstime;
void test(int k)
{
int i;
mpz_class m=1;
m <<=k;
m -= N;//calculate m =2^k-N;
mpz_class t=1;
t<<=N+1;
t%=m;
mpz_class x=2;
for(i=1;i<=k;i++){
x*=x;
mpz_class y=x>>k;
mpz_class z=x-(y<<k);
z+=N*y;x=z;
y=x>>k;
z=x-(y<<k);
z+=N*y;x=z;
while(x>=m)x-=m;
}
gc++;
if(x==t){//test passed
printf("%d\n",k);fflush(stdout);
fprintf(stderr,"Pass for %d\n",k);
}else{
if(gc>=100){
gc-=100;
time_t ct=time(NULL);
fprintf(stderr,"Searched for %d, time %d\n",k, (int)(ct-sstime));
}
}
}
int main(void)
{
int k;
FILE *f=fopen(IN,"r");
if(f==NULL){
fprintf(stderr, "Cannot open input file\n");
return -1;
}
sstime = time(NULL);
while(fscanf(f,"%d",&k)==1){//start to test for k
test(k);
}
fclose(f);
}
突然发现好像大家都没有注意到下面两个链接
https://oeis.org/A096502
以及
https://en.wikipedia.org/wiki/Riesel_number
里面的The dual Riesel number
当然没有google是不容易搜索到它们的
页:
1
[2]