/*
    Count the 4-dimensional polyhypercubes, treating mirror-images as
    different (A255487); but also count the symmetric ones, to get A068870.

    Each poly has a canonical form, such that two polys are equivalent
    if their canonical forms are identical. Compute it by examining all
    192 rotations of a form:
    - move each rotated poly close to the origin (coordinates are
      non-negative and minimal),
    - sort the sequence of hypercube coordinates, and
    - take the least of the 192 (using an arbitrary but consistent
      ordering).
*/

#include <algorithm>
#include <cstring>
#include <iostream>
#include <set>
using namespace std;

typedef unsigned long int coord; // four 8-bit numbers
static const int
    Shift1 = 8,
    Shift2 = 2 * Shift1,
    Shift3 = 3 * Shift1;
static const coord
    Wcoord = 1,
    Xcoord = Wcoord << Shift1,
    Ycoord = Xcoord << Shift1,
    Zcoord = Ycoord << Shift1,
    Endcoord = Zcoord << Shift1,
    Wmask = Xcoord - Wcoord,
    Xmask = Ycoord - Xcoord,
    Ymask = Zcoord - Ycoord,
    Zmask = Endcoord - Zcoord;

static const int MaxSize = 11;
static int PrevSize, Size;
    /* The program extends polys with PrevSize cells to make those with
       Size cells. PrevSize+1 == Size always. */

struct poly {
    coord coords[MaxSize];
    coord min, max;

    bool operator == (const poly& rhs) const {
	return memcmp(&coords[0], &rhs.coords[0], Size * sizeof (coord))
	       == 0;
    }
    bool operator < (const poly& rhs) const {
	return memcmp(&coords[0], &rhs.coords[0], Size * sizeof (coord))
	       < 0;
    }
    // Use comparators only on size-Size polys.
    /* set<poly> and canonicalize() use operator<
       It's an arbitrary but consistent ordering. */

    void move(int sx, coord displ);
    void updateExtrema(coord newco);
    poly canonicalize() const;
    void flipW();
    bool isSymmetric() const;
  private:
    void rotWtoX(poly& canon);
    void rotWtoY(poly& canon);
    void rotWtoZ(poly& canon);
    void rotXtoY(poly& canon);
    void rotYtoX(poly& canon);
    void rotXtoZ(poly& canon);
    void do2D(poly& canon); // do all 4 2-d rotations, WX plane
    void do3D(poly& canon); // do all 24 3-d rotations, WXY space
    void do4D(poly& canon); // do all 192 4-d rotations
};

static void printPoly(const poly& po) {
    for (int cx = 0; cx < Size; cx += 1) {
	coord co = po.coords[cx];
	cout << " ["
	     <<  (co & Wmask)            << ' '
	     << ((co & Xmask) >> Shift1) << ' '
	     << ((co & Ymask) >> Shift2) << ' '
	     << ((co & Zmask) >> Shift3) << ']';
    }
    cout << " {"
	     <<  (po.min & Wmask)            << ' '
	     << ((po.min & Xmask) >> Shift1) << ' '
	     << ((po.min & Ymask) >> Shift2) << ' '
	     << ((po.min & Zmask) >> Shift3) << '}';
    cout << " {"
	     <<  (po.max & Wmask)            << ' '
	     << ((po.max & Xmask) >> Shift1) << ' '
	     << ((po.max & Ymask) >> Shift2) << ' '
	     << ((po.max & Zmask) >> Shift3) << '}';
    cout << endl;
}

void poly::move(int sz, coord displ) {
    for (int cx = 0; cx < sz; cx += 1) {
	coords[cx] += displ;
    }
    min += displ;
    max += displ;
}

void poly::updateExtrema(coord newco) {
    if ((newco & Wmask) < (min & Wmask)) {
	min = (min & ~Wmask) | (newco & Wmask);
    }
    if ((newco & Wmask) > (max & Wmask)) {
	max = (max & ~Wmask) | (newco & Wmask);
    }
    if ((newco & Xmask) < (min & Xmask)) {
	min = (min & ~Xmask) | (newco & Xmask);
    }
    if ((newco & Xmask) > (max & Xmask)) {
	max = (max & ~Xmask) | (newco & Xmask);
    }
    if ((newco & Ymask) < (min & Ymask)) {
	min = (min & ~Ymask) | (newco & Ymask);
    }
    if ((newco & Ymask) > (max & Ymask)) {
	max = (max & ~Ymask) | (newco & Ymask);
    }
    if ((newco & Zmask) < (min & Zmask)) {
	min = (min & ~Zmask) | (newco & Zmask);
    }
    if ((newco & Zmask) > (max & Zmask)) {
	max = (max & ~Zmask) | (newco & Zmask);
    }
}

void poly::rotWtoX(poly& canon) {
    coord wmax = max & Wmask;
    for (int cx = 0; cx < Size; cx += 1) {
	coord& co = coords[cx];
	co = (co & (~ Wmask & ~ Xmask))
	   | ((wmax - (co & Wmask)) << Shift1)
	   | ((co & Xmask) >> Shift1);
    }
    max = (max & (~ Wmask & ~ Xmask))
	| (wmax << Shift1)
	| ((max & Xmask) >> Shift1);
    sort(&coords[0], &coords[Size]);
    if (*this < canon) {
	canon = *this;
    }
}

void poly::rotWtoY(poly& canon) {
    coord wmax = max & Wmask;
    for (int cx = 0; cx < Size; cx += 1) {
	coord& co = coords[cx];
	co = (co & (~ Wmask & ~ Ymask))
	   | ((wmax - (co & Wmask)) << Shift2)
	   | ((co & Ymask) >> Shift2);
    }
    max = (max & (~ Wmask & ~ Ymask))
	| (wmax << Shift2)
	| ((max & Ymask) >> Shift2);
    sort(&coords[0], &coords[Size]);
    if (*this < canon) {
	canon = *this;
    }
}

void poly::rotWtoZ(poly& canon) {
    coord wmax = max & Wmask;
    for (int cx = 0; cx < Size; cx += 1) {
	coord& co = coords[cx];
	co = (co & (~ Wmask & ~ Zmask))
	   | ((wmax - (co & Wmask)) << Shift3)
	   | ((co & Zmask) >> Shift3);
    }
    max = (max & (~ Wmask & ~ Zmask))
	| (wmax << Shift3)
	| ((max & Zmask) >> Shift3);
    sort(&coords[0], &coords[Size]);
    if (*this < canon) {
	canon = *this;
    }
}

void poly::rotXtoY(poly& canon) {
    coord xmax = max & Xmask;
    for (int cx = 0; cx < Size; cx += 1) {
	coord& co = coords[cx];
	co = (co & (~ Xmask & ~ Ymask))
	   | ((xmax - (co & Xmask)) << Shift1)
	   | ((co & Ymask) >> Shift1);
    }
    max = (max & (~ Xmask & ~ Ymask))
	| (xmax << Shift1)
	| ((max & Ymask) >> Shift1);
    sort(&coords[0], &coords[Size]);
    if (*this < canon) {
	canon = *this;
    }
}

void poly::rotYtoX(poly& canon) {
    coord ymax = max & Ymask;
    for (int cx = 0; cx < Size; cx += 1) {
	coord& co = coords[cx];
	co = (co & (~ Xmask & ~ Ymask))
	   | ((ymax - (co & Ymask)) >> Shift1)
	   | ((co & Xmask) << Shift1);
    }
    max = (max & (~ Xmask & ~ Ymask))
	| (ymax >> Shift1)
	| ((max & Xmask) << Shift1);
    sort(&coords[0], &coords[Size]);
    if (*this < canon) {
	canon = *this;
    }
}

void poly::rotXtoZ(poly& canon) {
    coord xmax = max & Xmask;
    for (int cx = 0; cx < Size; cx += 1) {
	coord& co = coords[cx];
	co = (co & (~ Xmask & ~ Zmask))
	   | ((xmax - (co & Xmask)) << Shift2)
	   | ((co & Zmask) >> Shift2);
    }
    max = (max & (~ Xmask & ~ Zmask))
	| (xmax << Shift2)
	| ((max & Zmask) >> Shift2);
    sort(&coords[0], &coords[Size]);
    if (*this < canon) {
	canon = *this;
    }
}

void poly::do2D(poly& canon) {
    rotWtoX(canon);
    rotWtoX(canon);
    rotWtoX(canon);
}

void poly::do3D(poly& canon) {
    do2D(canon);
    rotWtoY(canon);
    do2D(canon);
    rotXtoY(canon);
    do2D(canon);
    rotXtoY(canon);
    do2D(canon);
    rotYtoX(canon);
    do2D(canon);
    rotWtoY(canon);
    do2D(canon);
}

void poly::do4D(poly& canon) {
    do3D(canon);
    rotWtoZ(canon);
    do3D(canon);
    rotWtoZ(canon);
    do3D(canon);
    rotXtoZ(canon);
    do3D(canon);
    rotWtoZ(canon);
    do3D(canon);
    rotWtoZ(canon);
    do3D(canon);
    rotWtoZ(canon);
    do3D(canon);
    rotXtoZ(canon);
    do3D(canon);
}

void poly::flipW() {
    coord wmax = max & Wmask;
    for (int cx = 0; cx < Size; cx += 1) {
	coord& co = coords[cx];
	co = (co & ~ Wmask) | (wmax - (co & Wmask));
    }
    sort(&coords[0], &coords[Size]);
}

poly poly::canonicalize() const {
    poly po = *this;
    po.move(Size, - po.min);
    sort(&po.coords[0], &po.coords[Size]);
    poly canon = po;
    po.do4D(canon);
    return canon;
}

bool poly::isSymmetric() const {
    poly fl = *this;
    fl.flipW();
    fl = fl.canonicalize();
    return fl == *this;
}

typedef set<poly> polySet;
typedef polySet::iterator polyIter;

static polySet PolysA, PolysB;
static polySet *PolyP1, *PolyP2;

static void initPolys() {
    PolyP1 = &PolysA;
    PolyP2 = &PolysB;
    poly po;
    po.coords[0] = 0;
    po.min = 0;
    po.max = 0;
    PolyP2->insert(po);
}

static int SymQuan;

static void expandPoly(poly po) { // pass by copy
    po.move(PrevSize, Wcoord + Xcoord + Ycoord + Zcoord);
	// so that all adjacent cells exist
    set<coord> cellset;
    for (int cx = 0; cx < PrevSize; cx += 1) {
	cellset.insert(po.coords[cx]);
    }
    poly expo = po;
    coord& exco = expo.coords[PrevSize];
    for (int cx = 0; cx < PrevSize; cx += 1) {
	coord co = po.coords[cx];
	for (coord delta = Wcoord; delta != Endcoord; delta <<= Shift1) {
	    exco = co + delta;
	    if (cellset.insert(exco).second) {
		expo.updateExtrema(exco);
		poly canon = expo.canonicalize();
		bool newpoly;

#pragma omp critical
		newpoly = PolyP2->insert(canon).second;

		if (newpoly && canon.isSymmetric()) {
#pragma omp critical
		    SymQuan += 1;
		}

		expo.min = po.min;
		expo.max = po.max;
	    }
	    exco = co - delta;
	    if (cellset.insert(exco).second) {
		expo.updateExtrema(exco);
		poly canon = expo.canonicalize();
		bool newpoly;

#pragma omp critical
		newpoly = PolyP2->insert(canon).second;

		if (newpoly && canon.isSymmetric()) {
#pragma omp critical
		    SymQuan += 1;
		}

		expo.min = po.min;
		expo.max = po.max;
	    }
	}
    }
}


static polyIter PolyIters[1500000];
static void expandPolys() {
    SymQuan = 0;
    int sz = PolyP1->size();
    if (sz > 1500000) {
	throw "too big";
    }
    polyIter* pp = &PolyIters[0];
    for (polyIter it = PolyP1->begin(); it != PolyP1->end(); ++it) {
	/* Bah! It looks like OMP can't handle this kind of loop, so
	   copy the iterators to a plain array. */
	*pp++ = it;
    }

#pragma omp parallel for
    for (int ix = 0; ix < sz; ix += 1) {
	expandPoly(*PolyIters[ix]);
    }

}

static void printPolys() {
    for (polyIter it = PolyP2->begin(); it != PolyP2->end(); ++it) {
	printPoly(*it);
    }
}

int main() {
    initPolys();
    for (Size = 2; Size <= MaxSize; Size += 1) {
	PrevSize = Size - 1;
	swap(PolyP1, PolyP2);
	PolyP2->erase(PolyP2->begin(), PolyP2->end());
	expandPolys();
	int A255487 = PolyP2->size();
	cout << "size=" << Size << " A255487=" << A255487
	     << " SymQuan=" << SymQuan
	     << " A068870=" << ((A255487+SymQuan)/2)
	     << endl;
	// printPolys();
    }
    return 0;
}