Finds the wheel for 2 x 3 x 5 x ... x pi http://en.wikipedia.org/wiki/Wheel_factorization Complexity O(n*lg(n))

sievepattern.cc

#include <iostream>
#include <cstring>
#include <vector>
#define MAXN 50000000

using namespace std;

bool isp[MAXN];

vector<int> sieve(int q = 3) {
  vector<int> J;
  int P[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
  memset(isp, true, sizeof isp);
  for (int i = 0; i < q; ++i)
    for (int j = P[i] * P[i]; j < MAXN; j+= P[i])
      isp[j] = false;
  int jump = 1;
  for (int i = P[q] + 1; i < MAXN; ++i) {
    if (isp[i]) {
      J.push_back(jump);
      jump = 1;
    } else {
      ++jump;
    }
  }
  return J;
}

void display(const vector<int> &A, int limit = 1 << 30) {
  limit = min((int)A.size(), limit);
  if (limit > 0) {
    cout << A[0];
    for (int i = 1; i < limit; ++i)
      cout << " " << A[i];
  }
  cout << endl;
}

int pi[MAXN+1];

void kmpPreprocess(const vector<int> &P) {
  int i = 0, j = -1, m = P.size();
  pi[0] = -1;
  while(i < m) {
    while(j>=0 and P[i] != P[j])
      j = pi[j];
    i ++, j++;
    pi[i] = j;
  }
}

pair<int, int> searchPeriod(const vector<int> &J) {
  for (int i = 1; 2*i <= J.size(); ++i) {
    bool valid = true;
    for (int j = 2 * i; valid and j <= J.size(); j += i)
      if (pi[j] != j - i)
        valid = false;
    if (valid) {
      int per = 0;
      for (int j = 0; j < i; ++j)
        per+= J[j];
      return make_pair(i, per);
    }
  }
  return make_pair(-1, -1);
}

int main() {
  vector<int> J = sieve(3); // de 1 a 8
  kmpPreprocess(J);
  pair<int, int> P = searchPeriod(J);
  cout << P.first << " -> " << P.second << endl;
  display(J, P.first);
}