/* 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; }