Могу предложить мою C++ программу "матричное решето"
#include <cstdlib>
#include <iostream>
#include <math.h>
#include <ctime>
using namespace std;
main( )
{
/* 3S MATRIX SIEVE*/
/*FINDING PRIMES IN THE RANGE (N1;N2)*/
/*Authors: B.Sklyar, brs.sklr@mail.ru, D.Sklyar, I.Sklyar*/
/* N1>23; N2<2 000 000 000 000 000 000; N2-N1<500 000*/
unsigned long long int N1=/*enter N1*/; unsigned long long int N2 =/*enter N2*/;
unsigned long long int pr1=floor(N1/6); unsigned long long int pr2=ceil( N2/6);
int r=84000; int R2[r]; int rm=pr2-pr1;unsigned long long int S2[r];int r3, r4;
int q=84000; int R1[q] ; int qm=rm; unsigned long long int S1[q] ; int q2, q1;
for (q=1;q<qm;q++)
R1[q] =1;
for (r=1;r<rm;r++)
R2[r] =1;
unsigned long long int i, j, P1, P2, P3, P4, B, K;
unsigned long long int i2= sqrt( pr2/6)+2;
long long int j1, j2;
int l1=0;int l2=0;
for ( i=1;i<i2;i++)
{ j2=(pr2+i+1)/( 6*i+1)+1;j1=(pr1+i+1)/( 6*i+1);
B=5+5*( i-1); K=7+6*( i-1);
if ( i>j1) j1=i;
for(j=j1; j<j2;j++)
{ P1=B+K*( j-1);
if(( P1>pr1)&&( P1<pr2))
{ q1=P1-pr1; R1[ q1] =0;
} }
j2=(pr2-i+1)/( 6*i-1)+1;j1=(pr1-i+1)/( 6*i-1);
if (j1<1) j1=1;
B=5+7*( i-1); K=5+6*( i-1);
if ( i>j1-1) j1=i+1;
for(j=j1; j<j2;j++)
{ P2=B+K*( j-1);
if(( P2>pr1)&&( P2<pr2))
{ q2=P2-pr1; R1[ q2] =0;
} }
j2=(pr2+i+1)/( 6*i-1)+1;j1=(pr1+i+1)/( 6*i-1);
B=3+5*( i-1); K=5+6*( i-1);
if ( i>j1) j1=i;
for(j=j1; j<j2;j++)
{ P3=B+K*( j-1);
if(( P3>pr1)&&( P3<pr2))
{ r3=P3-pr1; R2[ r3]=0;
} }
j2=(pr2-i+1)/( 6*i+1)+1;j1=(pr1-i+1)/( 6*i+1);
B=7+7*(i-1); K=7+6*(i-1);
if ( i>j1) j1=i;
for(j=j1; j<j2;j++)
{ P4=B+K*( j-1);
if(( P4>pr1)&&( P4<pr2))
{ r4=P4-pr1; R2[ r4] =0;
} } }
cout<<"\ni2="<<i2<<";pr2="<<pr2<<" \n";
cout<<"\nP=pr1+q; pr1="<<pr1<<" \n";
for ( q=1;q<qm;q++) { S1[q] =R1[q]*((pr1+q)*6+5); if (S1[q]%5==0) continue; l1=l1+1;
cout<<"q="<<q<<"; Prime in S1[P]=6*P+5="<<S1[ q]<<" \n";}
cout<<"\nP=pr1+r; pr1="<<pr1<<" \n";
for ( r=1;r<rm;r++) { S2[r] =R2[r]*((pr1+r)*6+7);if (S2[r]%5==0) continue;l2=l2+1;
cout<<"r="<<r<<"; Prime in S2[P]=6*P+7="<<S2[ r]<<" \n";}
cout<<"l1="<<l1<<" \n";
cout<<"l2="<<l2<<" \n";
cout<<"number of Primes in the range (N1;N2) =l1+l2\n";
cout<<" run time (ms)=";
cout<<clock();
system("PAUSE");
return EXIT_SUCCESS;
}
Программа находит простые числа в заданном диапазоне (N1;N2), N2<2*10^18, N2-N1<500000 и выводит их в виде двух последовательностей S1=6*P+5 и S2=6*P+7