/* Author : Anand Krishnamoorthi HC_OEIS.cpp: Program to count number of Hamiltonian Cycles in P_m x P_n rectangular grids. */ /* Compilation: g++ -O2 -Wall A213813.cc -lstdc++ */ #if defined(__GCC__) #define FORCE_INLINE __attribute__((always_inline)) #elif defined(_MSC_VER) #define FORCE_INLINE __forceinline #define _SECURE_SCL 0 #define _HAS_ITERATOR_DEBUGGING 0 #else #define FORCE_INLINE inline #endif #include <vector> #include <iostream> #include <string> #include <map> #include <algorithm> #include <cassert> #include <cstdio> typedef unsigned long long ullong; struct Vertex { Vertex* path_end; unsigned char edges; unsigned char row; unsigned char col; }; struct Edge { Vertex* v1; Vertex* v2; Vertex* v1_path_end; Vertex* v2_path_end; bool created; FORCE_INLINE Edge(Vertex* v1, Vertex* v2) { this->v1 = v1; this->v2 = v2; created = false; assert(v1->edges != 2); // Make sure vertices are not fully connected.. if (v2->edges != 2) { // Check for cycles. This is a key semantic and performance constraint. if (v1 != v2->path_end) { v1_path_end = v1->path_end; v2_path_end = v2->path_end; /* sanity checks*/ assert(v1_path_end->path_end == v1); assert(v2_path_end->path_end == v2); created = true; // Update ends of the merged path to point to each other v1_path_end->path_end = v2_path_end; v2_path_end->path_end = v1_path_end; v1->edges++; v2->edges++; } } } FORCE_INLINE ~Edge() { if (created) { // Restore path end-points. v1_path_end->path_end = v1; v2_path_end->path_end = v2; v1->edges--; v2->edges--; } } FORCE_INLINE operator bool() const { return created; } }; struct LargeInt { static const ullong base = 1000ULL * 1000ULL * 1000ULL; std::vector<ullong> polynomial; FORCE_INLINE LargeInt(ullong v=0) { assert(v < base); polynomial.reserve(6); polynomial.clear(); polynomial.push_back(v); } FORCE_INLINE LargeInt& operator=(ullong v) { assert(v < base); polynomial.clear(); polynomial.reserve(6); polynomial.push_back(v); return *this; } bool operator==(const LargeInt& other) const { return polynomial == other.polynomial; } bool operator!=(const LargeInt& other) const { return polynomial != other.polynomial; } FORCE_INLINE LargeInt& operator+=(const LargeInt& other) { if (this == &other) { LargeInt tmp = other; *this += other; return *this; } if (polynomial.size() < other.polynomial.size()) polynomial.resize(other.polynomial.size()); ullong carry = 0; for (size_t i=0; i < polynomial.size(); ++i) { ullong term = polynomial[i] + carry; if (i < other.polynomial.size()) term += other.polynomial[i]; if (term >= base) { polynomial[i] = term - base; carry = 1; } else { polynomial[i] = term; carry= 0; } } if (carry) polynomial.push_back(carry); return *this; } void print(std::ostream& os) const { std::string num; num.reserve(1000); int digits_per_term = 0; ullong b = base; while (b/10) { ++digits_per_term; b = b / 10; } for(int p=0; p < int(polynomial.size())-1; ++p) { ullong term = polynomial[p]; for(int d=0; d < digits_per_term; ++d) { num.push_back((term % 10) + '0'); term /= 10; } } ullong term = polynomial.back(); while (term) { num.push_back((term % 10) + '0'); term /= 10; } if (num.empty()) num.push_back('0'); std::reverse(num.begin(), num.end()); os << num; } }; std::ostream& operator<<(std::ostream& os, const LargeInt& large_int) { large_int.print(os); return os; } FORCE_INLINE bool operator<(const Vertex& v1, const Vertex& v2) { return v1.path_end < v2.path_end; } typedef std::map<std::vector<Vertex>, LargeInt> Cache; typedef std::vector< std::vector<Vertex> > Grid; void solve(Grid& grid, const LargeInt& cur_count, int col, Cache& next_row_cache) { const int C = grid[0].size() - 1; int row = 0; if (col == C) { // Current row's been processed. next_row_cache[grid[row+1]] += cur_count; return; } Vertex* vertex = &grid[row][col]; Vertex* vertex_right = &grid[row][col+1]; Vertex* vertex_down = &grid[row+1][col]; if (vertex->edges == 2) return solve(grid, cur_count, col+1, next_row_cache); if (vertex->edges == 0) { // Vertex *must* be connected on both right and down. if (const Edge& edge_right = Edge(vertex, vertex_right)) if(const Edge& edge_down = Edge(vertex, vertex_down)) solve(grid, cur_count, col+1, next_row_cache); return; } // room->connections == 1. // Connect first to the right and explore. if (const Edge& edge_right = Edge(vertex, vertex_right)) solve(grid, cur_count, col+1, next_row_cache); // Then connect down and explore. if (const Edge& edge_down = Edge(vertex, vertex_down)) solve(grid, cur_count, col+1, next_row_cache); } FORCE_INLINE void init_row(int r, std::vector<Vertex>& row) { for(int c=0; c < row.size(); ++c) { Vertex* vertex = &row[c]; vertex->path_end = vertex; // This assignment simplifies connection. vertex->edges = 0; vertex->row = r; vertex->col = c; } row.back().edges = 2; } void solve(int R, int C) { std::vector< std::vector<Vertex> > grid; grid.resize(2, std::vector<Vertex>(C+1)); init_row(0, grid[0]); Cache cache; cache[grid[0]] = LargeInt(1); for (int row=0; row < R; ++row) { Cache next_row_cache; init_row(row+1, grid[1]); for(Cache::iterator itr = cache.begin(); itr != cache.end(); ++itr) { const LargeInt& count = itr->second; grid[0] = itr->first; solve(grid, count, 0, next_row_cache); } // Find the number of hamiltonian cycles till this row. // If the next now has just a single active path that starts and // ends at neighboring vertices, then that count is the number of // hamiltonian paths. // We choose vertex 0 and vertex 1. Any two vertices would work just as well. std::swap(grid[1][0].path_end, grid[1][1].path_end); LargeInt solutions = next_row_cache[grid[1]]; if (solutions != 0) std::cout << (row+1) << " " << solutions << std::endl; // move on to next row std::swap(cache, next_row_cache); std::swap(grid[0], grid[1]); } } int main(int argc, char** argv) { if (argc != 3) { std::cout << "usage: HC_OEIS cols rows" << std::endl; return 0; } int R = 0; int C = 0; sscanf(argv[1], "%d", &C); sscanf(argv[2], "%d", &R); if (C > 18) { std::cout <<"Solving for " << C <<" cols. May run out of memory\n"; } solve(R, C); }