#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 <boost/integer.hpp>

#include "mmalloc.hpp"
#include "timer.hpp"

using namespace std;
using boost::optional;

#ifndef BIT_WIDTH
#define BIT_WIDTH 16
#endif

typedef typename boost::uint_t<BIT_WIDTH + 1>::least State;
const State logState = BIT_WIDTH;
const State maxState = (State) 1 << logState;

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;
typedef array<uint8_t, maxState> pbitset;

void iterState(function<void(int)> f, optional<string> msg = optional<string>(), bool parallel = false) {
  PerfPrinter perfPrinter(msg);
  int numThreads=1;
  if (parallel) {
    numThreads = min((State) thread::hardware_concurrency(), maxState);
  }
  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, optional<string> msg = optional<string>(), bool parallel = false) {
  PerfPrinter perfPrinter(msg);
  auto msg2 = [=,&msg] (int i) {
    return msg ? (*msg + string(" ") + to_string(times-i) + string("/") + to_string(times)) : msg; };
  while (times--) {
    iterState(f, msg2(times), parallel);
  }
}

void init(Trans &t) {
  iterState([&](int s) {
      t[s] = update(s);
    }, (string) "single step transition table", true);
}

void findCycle(Trans &t, Trans &c, pbitset &reachable) {
  // compute reachability
  iterState([&](int s) {
      reachable[t[s]] = 1;
    }, (string) "reachability", true);
  // forward to t=maxState; now every state is in a cycle
  iterTrans(logState, [&](int s) {
      t[s] = t[t[s]];
    }, (string) "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];
    }, (string) "init cycles", true);
  iterTrans(logState, [&](int s) {
      c[s] = min<State>(c[s], c[t[s]]);
      t[s] = t[t[s]];
    }, (string) "cycles");
}

void cycleStat(Trans &c, pbitset &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++;
    }, (string) "basin & eden");

  // how long is the cycle, what is the actual minimal state
  { PerfPrinter perfPrinter((string) "cycle length");
  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)
  map<State, set<State>> ccyc;
  { PerfPrinter perfPrinter((string) "find duplicate cycles");
  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;
  };
  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);
  }
}

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>();
    pbitset *r = mmalloc<pbitset>();
    init(*t);
    findCycle(*t, *c, *r);
    cycleStat(*c, *r);
  }
  return 0;
}