リュカ・レーマー・テスト（メルセンヌ素数）を求めるC言語プログラム

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define N 10000
#define A 2
#define B 650
#define TXT "luc11b.txt"

/*
技術評論社　奥村晴彦　著
C言語による最新アルゴリズム事典　２３３０円＋税
P160〜161より引用・改変

A=始める数字
B=計算をここまでの数字までやる
TXT=実行結果を書き込むファイル名

AからBまでの素数のみリュカ・レーマー・テストをする
メルセンヌ素数は、2,3,5,7,13,17,19,31,61,89,107,127,
521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,
19937,21701,23209,44497,86243,110503,132049,216091


リュカ・レーマー・テストのアルゴリズムは

M(P)=2^P-1;　として、

x(1)=4;
x(i+1)=(x(i)^2-2) mod M(P)
x(P-1)=0　なら素数


*/

int prime(int p);
double gettime();
double st,pt;

main(){
   int p;
   int a,b,c,m,s,i;
   FILE *fp;

   pt=gettime();
   fp=fopen(TXT,"a");
   if(fp==NULL){
      printf("fail\n");
      exit(1);
   }
   fprintf(fp,"from A=%d to B=%d\n",A,B);
   fprintf(fp,"file :%s\n",TXT);
   fprintf(fp,"CALC START:\n");
   fclose(fp);

   for(m=A; m<=B; m++){
      s=m;
      i=2;
      while(i<=s){
         if(s<i*i){
            goto F1;
         }
         if(s%i==0){
            s/=i;
            i=2;
         }
         else{
            if(i==2) i++;
            else     i+=2;
         }
      }
F1:i=i;
      if(s==m){
         if(prime(m)){
             printf("p=%d は素数 t=%4.2f pt=%4.2f\n",m,gettime()-st,gettime()-pt);
             fp=fopen(TXT,"a");
             if(fp==NULL){
                printf("fail\n");
                exit(1);
             }
             fprintf(fp,"p=%d は素数 t=%4.2f pt=%4.2f\n",m,gettime()-st,gettime()-pt);
             fclose(fp);
           
         }
         else         printf("p=%d は素数ではない t=%4.2f pt=%4.2f\n",m,gettime()-st,gettime()-pt);
      }
   }
   fp=fopen(TXT,"a");
   if(fp==NULL){
      printf("fail\n");
      exit(1);
   }
   fprintf(fp,"CALC END:   t=%4.2f pt=%4.2f\n",gettime()-st,gettime()-pt);
   fclose(fp);
}


double gettime(){
   struct timeval tv;
   gettimeofday(&tv,NULL);
   return tv.tv_sec+(double)tv.tv_usec*1e-6;
}

int prime(int p){
   char a[N+1],x[N];
   int h,i,j,k,s;

   st=gettime();
   for(i=0;i<p;i++) a[i]=0;
   a[2]=1;
   for(k=2;k<p;k++){
      printf("%d/%d t=%4.2f pt=%4.2f\n",k,p,gettime()-st,gettime()-pt);
      for(i=0;i<p;i++){
         x[i]=a[i];
         a[i]=1;
      }
      a[1]=0;
      for(i=0;i<p;i++){
         if(x[i]){
            s=0;
            h=i;
            for(j=0;j<p;j++){
               s=(s>>1)+a[h]+x[j];
               a[h]=s&1;
               h=(h+1)%p;
            }
            if(s>1){

               while(a[h]){
                   a[h]=0;
                   h=(h+1)%p;
               }
               a[h]=1;
            }
         }
      }
   }
   a[p]=1-a[0];
   i=1;
   while(a[i]==a[0]) i++;
   return (i==p);
}