求一个最小的素数p，使得 p+2293是2的幂

mathe

第一个代码对于小素数p,求$2^x=n(mod p)$的解

1. #include <stdio.h>
   
2. #define PC 1000000
   
3. #define N 2293
   
4. char isprime[PC];
   
5. 
   
6. void init()
   
7. {
   
8.     long long i;
   
9.     isprime[0]=isprime[1]=0;
   
10.     for(i=2;i<PC;i++)isprime[i]=1;
    
11.     for(i=2;i*i<PC;i++){
    
12.         if(isprime[i]){
    
13.             int j;
    
14.             for(j=i*i;j<PC;j+=i)isprime[j]=0;
    
15.         }
    
16.     }
    
17. }
    
18. 
    
19. void solve(int p)
    
20. {
    
21.     int i,r,t=N%p,f=0;
    
22.     if(t==0)return;
    
23.     i=0,r=1;
    
24.     do{
    
25.        if(r==t){
    
26.            f=1;
    
27.            printf("%d,",i);
    
28.        }
    
29.        r*=2;if(r>=p)r-=p;i++;
    
30.        if(r==1){
    
31.            if(f==1) printf("%d\n",i);
    
32.            break;
    
33.        }
    
34.     }while(1);
    
35. }
    
36. 
    
37. int main()
    
38. {
    
39.     int i;
    
40.     init();
    
41.     for(i=3;i<PC;i++){
    
42.         if(isprime[i]){
    
43.             solve(i);
    
44.         }
    
45.     }
    
46.     return 0;
    
47. }
    

复制代码

mathe

第二个代码把上面代码输出简单处理一下，产生第一步筛选结果

1. #include <stdio.h>
   
2. #define IN "pp.out"
   
3. #define BUFC 100000000
   
4. #define MACC (8*BUFC)
   
5. char buf[BUFC];
   
6. 
   
7. #define BYTE(x)  ((x)>>3)
   
8. #define INDEX(x)  ((x)&7)
   
9. 
   
10. #define SET(x)   (buf[BYTE(x)]|=1<<INDEX(x))
    
11. #define GET(x)   (buf[BYTE(x)]&(1<<INDEX(x)))
    
12. 
    
13. int getd(int m){
    
14.     if(m%12==0)return 12;
    
15.     if(m%6==0)return 6;
    
16.     if(m%4==0)return 4;
    
17.     if(m%3==0)return 3;
    
18.     if(m%2==0)return 2;
    
19.     return 1;
    
20. }
    
21. 
    
22. int main()
    
23. {
    
24.     int d,m,h,s;
    
25.     FILE *in=fopen(IN,"r");
    
26.     if(in==NULL){
    
27.         fprintf(stderr, "Fail to open input file\n");
    
28.         return -1;
    
29.     }
    
30.     while(fscanf(in,"%d,%d",&d,&m)==2){
    
31.         int md=getd(m);
    
32.         if(md>1&&d%md!=1)continue;//skip data non compatible with 1(mod 12)
    
33.         m/=md;d%=m;
    
34.         if(m==1){
    
35.            fprintf(stderr,"Error, all data is filtered away\n");
    
36.            return -1;
    
37.         }
    
38.         //find s so that d+sm=1(mod 12)
    
39.         for(s=0;s<12;s++){
    
40.              if((d+s*m)%12==1)break;
    
41.         }
    
42.         if(s==12){
    
43.            fprintf(stderr,"Logic error\n");
    
44.            return -1;
    
45.         }
    
46.         for(h=(d+s*m-1)/12;h<MACC;h+=m){
    
47.             SET(h);
    
48.         }
    
49.     }
    
50.     fclose(in);
    
51.     for(h=0;h<MACC;h++){
    
52.        if(GET(h)==0){
    
53.             printf("%d\n",h*12+1);
    
54.        }
    
55.     }
    
56.     return 0;
    
57. }
    

复制代码

mathe

最后是伪素数判断

1. #include <gmpxx.h>
   
2. #include <stdio.h>
   
3. #include <time.h>
   
4. #define IN "ppc.out"
   
5. #define N  2293
   
6. 
   
7. int gc=0;
   
8. time_t sstime;
   
9. void test(int k)
   
10. {
    
11.     int i;
    
12.     mpz_class m=1;
    
13.     m <<=k;
    
14.     m -= N;//calculate m =2^k-N;
    
15.     mpz_class t=1;
    
16.     t<<=N+1;
    
17.     t%=m;
    
18.     mpz_class x=2;
    
19.     for(i=1;i<=k;i++){
    
20.         x*=x;
    
21.         mpz_class y=x>>k;
    
22.         mpz_class z=x-(y<<k);
    
23.         z+=N*y;x=z;
    
24.         y=x>>k;
    
25.         z=x-(y<<k);
    
26.         z+=N*y;x=z;
    
27.         while(x>=m)x-=m;
    
28.     }
    
29.     gc++;
    
30.     if(x==t){//test passed
    
31.           printf("%d\n",k);fflush(stdout);
    
32.          fprintf(stderr,"Pass for %d\n",k);
    
33.     }else{
    
34.        if(gc>=100){
    
35.          gc-=100;
    
36.          time_t ct=time(NULL);
    
37.          fprintf(stderr,"Searched for %d, time %d\n",k, (int)(ct-sstime));
    
38.        }
    
39.     }
    
40. }
    
41. 
    
42. int main(void)
    
43. {
    
44.     int k;
    
45.     FILE *f=fopen(IN,"r");
    
46.     if(f==NULL){
    
47.         fprintf(stderr, "Cannot open input file\n");
    
48.         return -1;
    
49.     }
    
50.     sstime = time(NULL);
    
51.     while(fscanf(f,"%d",&k)==1){//start to test for k
    
52.          test(k);
    
53.     }
    
54.     fclose(f);
    
55. }
    

复制代码

mathe

突然发现好像大家都没有注意到下面两个链接
https://oeis.org/A096502
以及
https://en.wikipedia.org/wiki/Riesel_number
里面的The dual Riesel number
当然没有google是不容易搜索到它们的