EULER AURAIT 303 ANS

/***********************************************************************************************
*  formules : recherche de p(x)=i*x4+j*x3+k*x2+l*x+m qui fournit n premiers pour x de 0 à n-1  *
*               Utilisation :  ...>formules.exe nnp i1 i2 j1 j2 k1 k2 l1 l2 m1 m2              *
***********************************************************************************************/
#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 ...
    if(i>0) return -i;
    return 1-i;
}
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]) {                         
            np[p]=i;     // si i est premier on le garde
            mp[p]=i*i;   // on initialise son plus petit multiple utile pour la suite
            p=p+1;
        }
    }
    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++) {             // tous les p-ièmes bits de pr[] sont mis à 1
        j=(np[i]-1)/8;                  // quand p est premier et les autres sont à 0
        k=np[i]-8*j;                    // on range les indicateurs 1 à 8 dans pr[0],
        pr[j] = pr[j] | 1<<(k-1);       // 9 à 16 dans pr[1], 17 à 24 dans pr[2], ...
    }                               // à noter : aucun np[i] ne vaut 8 ou un multiple
    return;                         // de 8 => ici, k est donc toujours entre 1 et 7 
}
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,ijk;
    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);
    prems();
    pFile=fopen("formules.txt","w");
    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 p(x)=i*x4+j*x3+k*x2+l*x+m qui fournit n nbs prems pour x de 0 à n-1 :");
    fprintf(pFile,"\n\n        i        j        k        l        m\n");    
    fclose(pFile);
    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);
    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)) {
                ijk=i+j+k;                          // si i+j+k+l est impair p(1) est pair sauf
                for (l=l1; l<=l2; l=s(l)) {         // pour m=2 et pour aller plus vite on peut
                    if((ijk+l)%2!=0) continue;      // éliminer cette formule, ce qui n'est pas
                    r=r1;                           // justifié pour p(1)=2 et m=2 : c'est rare
                    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;                //    quand px n'est pas trouvable
                            qx=(px-1)/8;                     // les bits 1 à 8 sont dans pr[0], 
                            rx=px-8*qx;                      // les bits 9 à 16 dans pr[1], ...
                            if( ( pr[qx]>>(rx-1) & 1 ) == 1 ) n=n+1;    // quand px est premier
                            else break;                                 // px n'est pas premier
                        } 
                        for(x=-1; x>-100; x--) {                           // et pour x négatif
                            px=(((i*x+j)*x+k)*x+l)*x+m;                      
                            if(px<2)   break;
                            if(px>npm) break;            
                            qx=(px-1)/8;                 
                            rx=px-8*qx;                  
                            if( ( pr[qx]>>(rx-1) & 1 ) == 1 ) n=n+1; 
                            else break;                              
                        } 
                        if(n>nm) {        // pour archiver choisir if(n>nm) ou, par exemple, if(n>37)
                            x=x+1;        // c'est la valeur la plus négative où p(x) est premier
                            mm=m+l*x+k*x*x+j*x*x*x+i*x*x*x*x;
                            lm=l+2*k*x+3*j*x*x+4*i*x*x*x;   // faire le décalage d'origine de x pour
                            km=k+3*j*x+6*i*x*x;             // la formule p(x) "définitive" qui fournit
                            jm=j+4*i*x;                     // n nombres premiers avec x de 0 à n-1
                            im=i;
                            x=n-1;                                 // pour examiner les formules symétriques
                            if(mm!=((((im*x+jm)*x+km)*x+lm)*x+mm)) goto f3;               // p(0) = p(n-1) ?
                            x=n-2;
                            if((im+jm+km+lm+mm)!=((((im*x+jm)*x+km)*x+lm)*x+mm)) goto f3; // p(1) = p(n-2) ?
                            if((n+1)/2<=nm) goto f4;        // la formule symétrique est-elle insuffisante ?
                            x=n/2;                          // pour enregistrer la formule symétrique avec
                            n=(n+1)/2;                      // uniquement la partie appropriée
                            mm=mm+lm*x+km*x*x+jm*x*x*x+im*x*x*x*x;  
                            lm=lm+2*km*x+3*jm*x*x+4*im*x*x*x;
                            km=km+3*jm*x+6*im*x*x;
                            jm=jm+4*im*x;
                            f3: pFile=fopen("formules.txt","a");
                            fprintf(pFile,"\n%9i%9i%9i%9i%9i    ( n = %2i ) ",im,jm,km,lm,mm,n);
                            fclose(pFile);
                            if(n>nm) nm=n;
                        }
                        f4: 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;
                        }
                    }
                }
            }
        }
    }
    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;
}