#include <cassert> #include <iomanip> #include <iostream> #include <limits> #include <list> #include <boost/multiprecision/cpp_int.hpp> using namespace std; typedef boost::multiprecision::uint1024_t value_type; typedef list<value_type> State_Description; static State_Description operator+(State_Description lhs, const value_type &rhs) { lhs.push_back(rhs); return lhs; } static State_Description & operator+=(State_Description &lhs, const value_type &rhs) { lhs.push_back(rhs); return lhs; } static ostream & operator<<(ostream &os, const State_Description &sd) { os << '<'; if(!sd.empty()) { auto st = sd.begin(); os << *st; while(++st != sd.end()) os << ", " << *st; } os << '>'; return os; } typedef list<pair<State_Description, value_type>> SD_Values; static SD_Values operator+(SD_Values lhs, const pair<State_Description, value_type> &rhs) { lhs.push_back(rhs); return lhs; } static SD_Values & operator+=(SD_Values &lhs, const pair<State_Description, value_type> &rhs) { lhs.push_back(rhs); return lhs; } static value_type sum(const SD_Values &lhs) { value_type sum = 0; for(auto lt = lhs.begin(); lt != lhs.end(); ++lt) sum += lt->second; return sum; } static value_type operator*(const SD_Values &lhs, const SD_Values &rhs) { value_type sum = 0; assert(lhs.size() == rhs.size()); auto lt = lhs.begin(); auto rt = rhs.begin(); while(lt != lhs.end() && rt != rhs.end()) { assert(lt->first == rt->first); sum += lt++->second * rt++->second; } return sum; } static ostream & operator<<(ostream &os, const SD_Values &sdv) { os << '{'; if(!sdv.empty()) { auto st = sdv.begin(); os << st->first << ':' << st->second; while(++st != sdv.end()) os << ", " << st->first << ':' << st->second; } os << '}'; return os; } static value_type factorial(const int64_t n, const value_type accumulator = 1) { if(n < 2) return accumulator; return factorial(n - 1, accumulator * n); } static value_type combination(const int64_t n, const int64_t r) { return factorial(n) / (factorial(n - r) * factorial(r)); } static value_type permutation(const int64_t n, const int64_t r) { return factorial(n) / factorial(n - r); } static void bw_gen_states(SD_Values &sdv, State_Description sd, const value_type &multiplier, const int64_t num_blocks, int64_t largest_stack_size, int64_t num_of_largest) { sd += largest_stack_size; if(largest_stack_size == 1) { sdv += make_pair(sd, multiplier); return; } const value_type stack_0 = multiplier * permutation(num_blocks, largest_stack_size) / num_of_largest++; const int64_t remaining_blocks = num_blocks - largest_stack_size; if(!remaining_blocks) { sdv += make_pair(sd, stack_0); return; } if(remaining_blocks < largest_stack_size) { largest_stack_size = remaining_blocks; num_of_largest = 1; } do { bw_gen_states(sdv, sd, stack_0, remaining_blocks, largest_stack_size--, num_of_largest); num_of_largest = 1; } while(largest_stack_size != 0); } static void bw_gen_states(SD_Values &sdv, const int64_t num_blocks) { for(int64_t largest_stack_size = num_blocks; largest_stack_size != 0; --largest_stack_size) bw_gen_states(sdv, State_Description(), 1, num_blocks, largest_stack_size, 1); } static value_type bw_num_states(const int64_t num_blocks) { SD_Values sdv; bw_gen_states(sdv, num_blocks); return sum(sdv); } static void bw_gen_sap(SD_Values &sdv, State_Description sd, const int64_t num_blocks, int64_t largest_stack_size, const int64_t num_preceding_stacks) { sd += largest_stack_size; if(largest_stack_size == 1) { const int64_t num_stacks = num_preceding_stacks + num_blocks - 1; sdv += make_pair(sd, num_stacks * num_stacks - num_blocks); return; } const int64_t remaining_blocks = num_blocks - largest_stack_size; if(!remaining_blocks) { sdv += make_pair(sd, num_preceding_stacks * num_preceding_stacks); return; } if(remaining_blocks < largest_stack_size) largest_stack_size = remaining_blocks; do { bw_gen_sap(sdv, sd, remaining_blocks, largest_stack_size--, num_preceding_stacks + 1); } while(largest_stack_size != 0); } static void bw_gen_sap(SD_Values &sdv, const int64_t num_blocks) { for(int64_t largest_stack_size = num_blocks; largest_stack_size != 0; --largest_stack_size) bw_gen_sap(sdv, State_Description(), num_blocks, largest_stack_size, 1); } static value_type bw_num_sap(const int64_t num_blocks) { SD_Values sdv; SD_Values sda; bw_gen_states(sdv, num_blocks); bw_gen_sap(sda, num_blocks); // cerr << sdv << endl; // cerr << sda << endl; return sdv * sda; } static tuple<int64_t, value_type, value_type, value_type> bw_num_states_sap(const int64_t num_blocks) { SD_Values sdv; SD_Values sda; bw_gen_states(sdv, num_blocks); bw_gen_sap(sda, num_blocks); // cerr << sdv << endl; // cerr << sda << endl; return make_tuple(int64_t(sdv.size()), sum(sdv), sum(sda), sdv * sda); } class comma_numpunct : public numpunct<char> { protected: virtual char do_thousands_sep() const { return ','; } virtual string do_grouping() const { return "\03"; } }; int main(int argc, char **argv) { // this creates a new locale based on the current application default // (which is either the one given on startup, but can be overriden with // locale::global) - then extends it with an extra facet that // controls numeric output. locale comma_locale(locale(), new comma_numpunct()); // tell cout to use our new locale. cerr.imbue(comma_locale); cout.imbue(comma_locale); cerr.precision(numeric_limits<value_type>::max_digits10); cout.precision(numeric_limits<value_type>::max_digits10); for(int i = 1; i < argc; ++i) { const int64_t num_blocks = atoi(argv[i]); auto num_states_sap = bw_num_states_sap(num_blocks); cout << num_blocks << ' ' << get<3>(num_states_sap) << std::endl; } }