Generate primes < 10^9

primegen.cpp

#include <iostream>
#include <bitset>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>

using namespace std;

typedef long long ll;
const ll maxp = 1000000000;
const int blocksz = 2*3*5*7*11*13*17;
const int firstp = 7;

//const int wheelmod = 30;
//const int wheelsz = 8;
//const int wheel[] = { 6,4,2,4,2,4,6,2 };
//int rems[] = { 1,7,11,13,17,19,23,29 };
//int remtobit[32];
const int wheelmod = 210;
const int wheelsz = 48;
const int wheel[] =
  { 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2
  , 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2
  };
int rems[] =
  { 1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79
  , 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143, 149, 151
  , 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209
  };
int remtobit[256];

//const int bits = blocksz/wheelmod * wheelsz;
bitset<blocksz> presieve, sieve;
const int maxsievep = 31627;
bitset<maxsievep+1> sieve2;

int spsz=0;
int spval[4096];
int spsq[4096];
int spmult[4096];
int spfirstb[4096];
static unsigned char spwi[4096];
int spwheelprod[4096][wheelsz];

const int bufsz = 4*1024*1024;
unsigned char buf[bufsz];
int bufpos = 0;

FILE *outfile;
void flush() {
  fwrite(buf, bufpos, 1, outfile);
  bufpos = 0;
}

int cnt=0;
inline void output(int p) {
  //cnt++;
  //if (p >= maxp) return;
  //cout << p << "\n";
  *(int*)(buf+bufpos) = p;
  bufpos += 4;
  if (bufpos == bufsz) flush();
}

int main() {
  ios_base::sync_with_stdio(false);
  outfile = fopen("huge-prime-list", "wb");
  setbuf(outfile, 0);
  for (int i = 0; i < wheelmod; ++i)
    remtobit[i] = -1;
  for (int i = 0; i < wheelsz; ++i)
    remtobit[rems[i]] = i;

  for (int i = 2; i <= maxsievep; i++) {
    if (sieve2[i]) continue;
    spval[spsz] = i;
    spsq[spsz] = i*i;
    spmult[spsz] = (i*i) % blocksz;
    spwi[spsz] = remtobit[i%wheelmod];
    spfirstb[spsz] = (i*i) / blocksz;
    for (int j = 0; j < wheelsz; ++j)
      spwheelprod[spsz][j] = wheel[j]*i;
    ++spsz;
    for (int j = i*i; j <= maxsievep; j += i)
      sieve2.set(j);
  }
  for (int i = 0; i < firstp; ++i)
    output(spval[i]);

  for (int i = 0; i < firstp; ++i) {
    int p = spval[i];
    for (int j = p; j < blocksz; j += p)
      presieve.set(j);
  }

  int blocks = (maxp+blocksz-1)/blocksz;
  sieve = presieve;
  sieve.set(1);
  for (int b = 0; b < blocks; ++b) {
    int lo = b*blocksz;
    for (int i = firstp; i < spsz; ++i) {
      if (spfirstb[i] > b) break;
      const int *wp = spwheelprod[i];
      int j = spmult[i];
      int idx = spwi[i];
      int k = wheelsz-idx;
      while (j < blocksz) {
        while(k--) {
          if (j >= blocksz) break;
          sieve.set(j);
          j += wp[idx];
          ++idx;
        }
        if (j < blocksz) idx=0,k=wheelsz;
        else break;
      }
      spmult[i] = j % blocksz;
      spwi[i] = idx;
    }
#define FOO(x) if (!sieve[i]) output(lo + i); i += x;
    int i = 1;
    while (i < blocksz) {
      FOO(10); FOO(2); FOO(4); FOO(2); FOO(4); FOO(6); FOO(2); FOO(6);
      FOO(4); FOO(2); FOO(4); FOO(6); FOO(6); FOO(2); FOO(6); FOO(4); FOO(2);
      FOO(6); FOO(4); FOO(6); FOO(8); FOO(4); FOO(2); FOO(4); FOO(2); FOO(4);
      FOO(8); FOO(6); FOO(4); FOO(6); FOO(2); FOO(4); FOO(6); FOO(2); FOO(6);
      FOO(6); FOO(4); FOO(2); FOO(4); FOO(6); FOO(2); FOO(6); FOO(4); FOO(2);
      FOO(4); FOO(2); FOO(10); FOO(2);
    }
    sieve = presieve;
  }
  //cout << cnt << endl;
  flush();
}