FORMULES POUR NOMBRES PREMIERS

#include <iostream>
#include <math.h>
#include <time.h>
int nnp, *np;      // nombre de nombres premiers et liste des nnp premiers nombres premiers
unsigned char *pr; // suite d'octets utilisés comme suite de bits pour marquer tous les nombres
int npr, npm;      // nombre d'octets utiles dans la suite pr[] et simplement npm=np[nnp] 
int s(int i) {  // le suivant de i dans la suite de tous les entiers : 0 1 -1 2 -2 3 -3 4 -4 ...
    int x;
    x=-i;
    if(x<0) return x;
    return x+1;
}
void prems(){                 // calcul des nnp premiers nombres premiers dans np[i],i=1,n
    int i,j,k,m,p,t;          // et de chaque i-ième bit de pr[] selon i premier ou composé
    unsigned char zero='\x00';
    int *mp;
    np=(int *)malloc((nnp+1)*sizeof(int));           // np[i] utile pour i de 1 à nnp
    mp=(int *)malloc((nnp+1)*sizeof(int));           // mp[i] utile pour i de 3 à nnp
    np[0]=1;  np[1]=2;  np[2]=3;  np[3]=5;           // np[0], mp[0], mp[1] et mp[2]    
    mp[0]=1;  mp[1]=4;  mp[2]=9;  mp[3]=5*5;         // ne serviront jamais à rien !!
    p=4;                       // c'est le rang du prochain nombre premier recherché
    m=nnp+1;                   // on s'arrêtera quand np[p] sera trouvé pour p=nnp
    t=2;                       // les multiples de 2 et 3 seront évités : on fera
    for(i=7; p<m; i=i+t) {     // i = 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, ... 
        t=6-t;                 // maintenant np[1] et np[2] ne seront plus utilisés           
        k=3;                   // comme annoncé, on commence à 3 ( j = np[3] = 5 )
        for(;;) { 
            j=np[k];                // à tout moment de la recherche on tient à jour      
            if( i==mp[k] ) break;   // dans mp[k] le plus petit multiple de np[k]      
            if( j*j>=i )   break;   // qu'il faudra éviter ensuite quand un prochain
            if(i>mp[k]) mp[k]=mp[k]+2*j;  // nombre candidat examiné i aura sa valeur
            k=k+1;                  // les nombres pairs ne sont pas examinés => +2*j 
        }                                     
        if(i!=mp[k]) {                         
            mp[p]=i*i;   // si i est premier on le garde
            np[p++]=i;   // et on initialise son plus petit multiple utile
        }
    }
    free(mp);     // pour libérer la mémoire du tableau auxiliaire mp[]
    npm=np[nnp];  // c'est le dernier nombre premier calculé et rangé
    npr=npm/8+1;  // nombre d'octets utiles dans pr[] ( npm n'est pas multiple de 8 ! )
    pr=(unsigned char *)malloc(npr*sizeof(unsigned char));
    for(i=0; i<npr;  i++) pr[i]=zero;   
    for(i=1; i<=nnp; i++) {
        j=np[i]/8;                      // tous les n-ièmes bits de pr[] sont mis à 1
        k=np[i]-8*j;                    // quand n est premier et les autres sont à 0
        pr[j] = pr[j] | 1<<(k-1);  
    }
    return;
}
int rang(int n){           // pour np[1] <= n <= np[nnp]
    int i,i1,i2;           // on calcule le rang i tel que
    i1=1;                  // np[i] soit égal à n ou bien
    i2=nnp;                // le plus voisin supérieur à n avec 
    while( i1 <= i2 ) {    // ici une méthode en dichotomie habituelle
        i=(i1+i2)/2;       
        if(n==np[i]) return i;     // quand n est trouvé
        if(n<np[i])  i2=i-1;
        else         i1=i+1;
    }
    return i1;  // quand n n'est pas trouvé ( maintenant : i1 > i2 ! )
}
int main(int argc, char *argv[]) {
    FILE * pFile;
    int i,j,k,l,m,n,r,r1,im,jm,km,lm,mm,nm;
    int i1,i2,j1,j2,k1,k2,l1,l2,m1,m2;
    int x,px,qx,rx,vu;
    if(argc!=12) {
        printf("\nIl faut faire :\n\nformules nnp i1 i2 j1 j2 k1 k2 l1 l2 m1 m2\n\n");
        system("pause"); 
        return 1;
    }
    nnp=atoi(argv[1]);
    i1=atoi(argv[2]);  i2=atoi(argv[3]); 
    j1=atoi(argv[4]);  j2=atoi(argv[5]); 
    k1=atoi(argv[6]);  k2=atoi(argv[7]); 
    l1=atoi(argv[8]);  l2=atoi(argv[9]); 
    m1=atoi(argv[10]); m2=atoi(argv[11]); 
    time_t t1=time(NULL);
    pFile=fopen("formules.txt","w");
    prems();
    printf("\nCalcul de %i nombres premiers termin\x82 !\n\n",nnp);
    printf(" %i < i < %i   %i < j < %i   ",i1,i2,j1,j2);
    printf("%i < k < %i   %i < l < %i   %i < m < %i \n\n",k1,k2,l1,l2,m1,m2);
    fprintf(pFile,"\nLe %i-ième nombre premier vaut : %i\n",nnp,npm);
    fprintf(pFile,"\nPour :\n");
    fprintf(pFile,"\n %i < i < %i ",i1,i2);
    fprintf(pFile,"\n %i < j < %i ",j1,j2);
    fprintf(pFile,"\n %i < k < %i ",k1,k2);
    fprintf(pFile,"\n %i < l < %i ",l1,l2);
    fprintf(pFile,"\n %i < m < %i ",m1,m2);
    fprintf(pFile,"\n\nOn trouve :");
    fprintf(pFile,"\n\n  i  j  k  l  m    n\n");    
    fclose(pFile);
    if(m1<2)     m1=2;        // m est nécessairement un nombre premier
    if(m2>npm)   m2=npm;      // il doit être l'un des :  np[i],i=1,nnp
    r1=rang(m1);
    nm=0;
    vu=0;
    for (i=i1; i<=i2; i=s(i)) { 
        for (j=j1; j<=j2; j=s(j)) {
            for (k=k1; k<=k2; k=s(k)) {
                for (l=l1; l<=l2; l=s(l)) {
                    if(i==0 && j==0 && k==0 && l==0) goto fin;        // p(x) = cte est refusée
                    r=r1;
                    for (m=np[r]; m<=m2; m=np[++r]) {      
                        n=0;         // on va compter le nombre de nombres premiers consécutifs
                        for(x=0; x<=100; x++) {   // pour : p(x) = i.x4 + j.x3 + k.x2 + l.x + m
                            px=(((i*x+j)*x+k)*x+l)*x+m;                       // de x = 0 à 100
                            if(px<2)   break;
                            if(px>npm) break;                         // px n'est pas trouvable
                            qx=px/8;
                            rx=px-8*qx;
                            if( ( pr[qx]>>(rx-1) & 1 ) == 1 ) n=n+1;    // quand px est premier
                            else break;                                 // px n'est pas premier
                        } 
                        if(n>nm) {               // si la formule est meilleure on l'enregistre
                            pFile=fopen("formules.txt","a");
                            fprintf(pFile,"\n  %i  %i  %i  %i  %i  ( %i ) ",i,j,k,l,m,n);
                            fclose(pFile);
                            im=i;
                            jm=j;
                            km=k;
                            lm=l;
                            mm=m;
                            nm=n;
                        }
                        vu=vu+1;
                        if(vu==50000000) {      // avancement actuel de la recherche
                        printf("\ron est \x85  i : %i  j : %i  k : %i  ",i,j,k);
                        printf("l : %i  m : %i  et  nm : %i     ",l,m,nm);
                        vu=0;
                        }
                    }
                    fin: n=0;
                }
            }
        }
    }
    time_t t2=time(NULL);
    pFile=fopen("formules.txt","a");
    fprintf(pFile,"\n\nExécution en %i secondes\n",(int)difftime(t2,t1));
    fclose(pFile);
    pFile=fopen("verif.bat","w");  // pour aider à vérifier la meilleure formule
    fprintf(pFile,"@echo off\nverif.exe %i %i %i %i %i %i",nnp,im,jm,km,lm,mm);
    fclose(pFile);
    free(np);
    free(pr);
    return 0;
}