找回密码
 欢迎注册
楼主: hujunhua

[提问] 求一个最小的素数p,使得 p+2293是2的幂

[复制链接]
发表于 2016-11-4 19:13:54 | 显示全部楼层
第一个代码对于小素数p,求$2^x=n(mod p)$的解
  1. #include <stdio.h>
  2. #define PC 1000000
  3. #define N 2293
  4. char isprime[PC];

  5. void init()
  6. {
  7.     long long i;
  8.     isprime[0]=isprime[1]=0;
  9.     for(i=2;i<PC;i++)isprime[i]=1;
  10.     for(i=2;i*i<PC;i++){
  11.         if(isprime[i]){
  12.             int j;
  13.             for(j=i*i;j<PC;j+=i)isprime[j]=0;
  14.         }
  15.     }
  16. }

  17. void solve(int p)
  18. {
  19.     int i,r,t=N%p,f=0;
  20.     if(t==0)return;
  21.     i=0,r=1;
  22.     do{
  23.        if(r==t){
  24.            f=1;
  25.            printf("%d,",i);
  26.        }
  27.        r*=2;if(r>=p)r-=p;i++;
  28.        if(r==1){
  29.            if(f==1) printf("%d\n",i);
  30.            break;
  31.        }
  32.     }while(1);
  33. }

  34. int main()
  35. {
  36.     int i;
  37.     init();
  38.     for(i=3;i<PC;i++){
  39.         if(isprime[i]){
  40.             solve(i);
  41.         }
  42.     }
  43.     return 0;
  44. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-11-4 19:15:01 | 显示全部楼层
第二个代码把上面代码输出简单处理一下,产生第一步筛选结果
  1. #include <stdio.h>
  2. #define IN "pp.out"
  3. #define BUFC 100000000
  4. #define MACC (8*BUFC)
  5. char buf[BUFC];

  6. #define BYTE(x)  ((x)>>3)
  7. #define INDEX(x)  ((x)&7)

  8. #define SET(x)   (buf[BYTE(x)]|=1<<INDEX(x))
  9. #define GET(x)   (buf[BYTE(x)]&(1<<INDEX(x)))

  10. int getd(int m){
  11.     if(m%12==0)return 12;
  12.     if(m%6==0)return 6;
  13.     if(m%4==0)return 4;
  14.     if(m%3==0)return 3;
  15.     if(m%2==0)return 2;
  16.     return 1;
  17. }

  18. int main()
  19. {
  20.     int d,m,h,s;
  21.     FILE *in=fopen(IN,"r");
  22.     if(in==NULL){
  23.         fprintf(stderr, "Fail to open input file\n");
  24.         return -1;
  25.     }
  26.     while(fscanf(in,"%d,%d",&d,&m)==2){
  27.         int md=getd(m);
  28.         if(md>1&&d%md!=1)continue;//skip data non compatible with 1(mod 12)
  29.         m/=md;d%=m;
  30.         if(m==1){
  31.            fprintf(stderr,"Error, all data is filtered away\n");
  32.            return -1;
  33.         }
  34.         //find s so that d+sm=1(mod 12)
  35.         for(s=0;s<12;s++){
  36.              if((d+s*m)%12==1)break;
  37.         }
  38.         if(s==12){
  39.            fprintf(stderr,"Logic error\n");
  40.            return -1;
  41.         }
  42.         for(h=(d+s*m-1)/12;h<MACC;h+=m){
  43.             SET(h);
  44.         }
  45.     }
  46.     fclose(in);
  47.     for(h=0;h<MACC;h++){
  48.        if(GET(h)==0){
  49.             printf("%d\n",h*12+1);
  50.        }
  51.     }
  52.     return 0;
  53. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-11-4 19:16:25 | 显示全部楼层
最后是伪素数判断
  1. #include <gmpxx.h>
  2. #include <stdio.h>
  3. #include <time.h>
  4. #define IN "ppc.out"
  5. #define N  2293

  6. int gc=0;
  7. time_t sstime;
  8. void test(int k)
  9. {
  10.     int i;
  11.     mpz_class m=1;
  12.     m <<=k;
  13.     m -= N;//calculate m =2^k-N;
  14.     mpz_class t=1;
  15.     t<<=N+1;
  16.     t%=m;
  17.     mpz_class x=2;
  18.     for(i=1;i<=k;i++){
  19.         x*=x;
  20.         mpz_class y=x>>k;
  21.         mpz_class z=x-(y<<k);
  22.         z+=N*y;x=z;
  23.         y=x>>k;
  24.         z=x-(y<<k);
  25.         z+=N*y;x=z;
  26.         while(x>=m)x-=m;
  27.     }
  28.     gc++;
  29.     if(x==t){//test passed
  30.           printf("%d\n",k);fflush(stdout);
  31.          fprintf(stderr,"Pass for %d\n",k);
  32.     }else{
  33.        if(gc>=100){
  34.          gc-=100;
  35.          time_t ct=time(NULL);
  36.          fprintf(stderr,"Searched for %d, time %d\n",k, (int)(ct-sstime));
  37.        }
  38.     }
  39. }

  40. int main(void)
  41. {
  42.     int k;
  43.     FILE *f=fopen(IN,"r");
  44.     if(f==NULL){
  45.         fprintf(stderr, "Cannot open input file\n");
  46.         return -1;
  47.     }
  48.     sstime = time(NULL);
  49.     while(fscanf(f,"%d",&k)==1){//start to test for k
  50.          test(k);
  51.     }
  52.     fclose(f);
  53. }
复制代码

点评

另外需要注意函数test中假设了$2^k>N$,所以对于太小的k是不能使用这个函数的  发表于 2016-11-4 21:39
其实前面两个小程序相对都很快,最主要需要对最后一个程序优化  发表于 2016-11-4 19:28
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-11-6 08:13:49 | 显示全部楼层
突然发现好像大家都没有注意到下面两个链接
https://oeis.org/A096502
以及
https://en.wikipedia.org/wiki/Riesel_number
里面的The dual Riesel number
当然没有google是不容易搜索到它们的

评分

参与人数 1威望 +6 金币 +6 贡献 +6 经验 +6 鲜花 +6 收起 理由
northwolves + 6 + 6 + 6 + 6 + 6 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 03:20 , Processed in 0.057240 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表