#include <assert.h>
#include <stdlib.h>
#include <string.h>

#include <array>
#include <bitset>
#include <functional>
#include <iostream>
#include <list>
#include <map>
#include <set>
#include <thread>

#include <sys/mman.h>

using namespace std;

#ifndef BIT_WIDTH
#define BIT_WIDTH 31
#endif

typedef uint32_t State;
const State logState = BIT_WIDTH;
const State maxState = 1 << logState;

// # of bits of largest memory fetch issued; machine-specific; used to
// garantue that sub-byte sized access of different threads never
// address the same word
const uint maxWordSize = 128;

bitset<8> rule(110);

State update(State s) {
  State r(0);
  bitset<logState> b(s);
  for (unsigned i=0; i<logState; i++)
    r |= rule[1 * b[(i + logState - 1) % logState]
	    + 2 * b[i]                
	    + 4 * b[(i + 1) % logState]] << i;
  return r;
}

typedef array<State, maxState> Trans;

void iterState(function<void(int)> f, bool parallel = false) {
  int numThreads=1;
  if (parallel) {
    numThreads = min(thread::hardware_concurrency(),
		     maxState / maxWordSize);
  }
  list<thread*> tasks;
  for (int t=0; t<numThreads; t++) {
    tasks.push_front(new thread([=]{
	  for (State s = maxState / numThreads * t; 
	       s < ((t == numThreads - 1) ? maxState : (maxState / numThreads * (t+1)));
	       s++)
	    f(s);
	}));
  }
  for (; !tasks.empty(); tasks.front()->join(), delete tasks.front(), tasks.pop_front());    
}

void iterTrans(int times, function<void(int)> f, char *msg = nullptr, bool parallel = false) {
  if (msg) cerr << msg << times << " left   \r";
  while (times--) {
    iterState(f, parallel);
    if (msg) cerr << msg << times << " left   \r";
  }
  if (msg) cerr << "                                  \r";
}

void init(Trans &t) {
  iterState([&](int s) {
      t[s] = update(s); }, true);
}

void findCycle(Trans &t, Trans &c, bitset<maxState> &reachable) {
  // compute reachability
  iterState([&](int s) {
      reachable[t[s]] = 1;
    }, true);
  // forward to t=maxState; now every state is in a cycle
  iterTrans(logState, [&](int s) {
      t[s] = t[t[s]]; }, (char*) "fwd time: ");
  // compute loop id (lowest occuring state no): go through the loop again and
  // record the lowest occuring state number
  iterState([&](int s) { 
      c[s] = t[s]; }, true);
  iterTrans(logState, [&](int s) {
      c[s] = min(c[s], c[t[s]]);
      t[s] = t[t[s]]; }, (char*) "cycles: ");
}

void cycleStat(Trans &c, bitset<maxState> &reachable) {
  struct Stat {
    State basin, len, eden, minState;
    Stat() : basin(0), len(0), eden(0), minState(maxState) {}
  };
  map<State, Stat> cyc;

  // How big is the basin of attraction?
  // How many garden of eden states does it contain?
  iterState([&](int s) {
      cyc[c[s]].basin++;
      if (!reachable[s])
	cyc[c[s]].eden++;
    });

  // how long is the cycle, what is the actual minimal state
  for (auto i : cyc) {
    Stat &s = cyc[i.first];
    State cur, start;
    cur = start = c[i.first]; 
    do {
      s.len++;
      s.minState = min(s.minState, cur);
      cur=update(cur);
    } while (cur != start);
  }

  // find duplicates cycles (only bitshifted)
  auto canonize = [](State s) {
    State cs = s;
    for (State i=0; i<logState; i++)
      cs = min(cs, (((s<<i) | (s>>(logState-i))) & (maxState - 1)));
    return cs;
  };
  map<State, set<State>> ccyc;
  for (auto i : cyc) {
    ccyc[canonize(i.second.minState)].insert(i.first);
  }

  // print it
  for (auto i : ccyc) {
    Stat &s = cyc[*(i.second.begin())];
    cout << bitset<logState>(i.first) << "\t"
	 << i.second.size() << "\t"
	 << s.len << "\t"
	 << s.basin << "\t"
	 << s.eden  << "\t"
	 << i.first << endl;
  }
}

void print(Trans &t) {
  for (auto s : t)
    cout << bitset<logState>(s) << endl;
}

void printTraj(State s, int count) {
  while (count--) {
    cout << bitset<logState>(s) << endl;
    s = update(s);
  }
}

template<typename T>
T* mmalloc() {
  void *p = MAP_FAILED;
  for (auto flag : {MAP_POPULATE | MAP_HUGETLB, MAP_POPULATE, 0})
    if (p == MAP_FAILED)
      p = mmap(NULL, sizeof(Trans), PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS | flag, -1, 0);
  assert(p != MAP_FAILED);
  return (T*) p;
}

int main(int argc, char **argv) {
  assert(argc >= 2);
  rule = atoi(argv[1]);
  if (!strcmp(argv[2], "traj")) {
    assert(argc == 5);
    printTraj(atoi(argv[3]), atoi(argv[4]));
  }
  if (!strcmp(argv[2], "cycle")) {
    Trans *t = mmalloc<Trans>(),
      *c = mmalloc<Trans>();
    bitset<maxState> *r = new bitset<maxState>(0);
    init(*t);
    findCycle(*t, *c, *r);
    cycleStat(*c, *r);
  }
  return 0;
}