// Algorithms from Theorem 2 // including the prefix sums AND the memory optimization // Using the GMPlib big integers from the CGAL library // Runtime O(n^4), Memory O(n^3) times the bigint-overhead (roughly another O(nlogn)). // Just compile and run it to get a quick demo or look at the main function for some examples: // Compilation steps with CGAL: // Run: cgal_create_cmake_script // To the file CMakeLists.txt add this line: set(CMAKE_CXX_FLAGS \${CMAKE_CXX_FLAGS} " -std=c++11") // Run: cmake . // Run: make // Execute: ./motzkin_leftright_bignum #include #include #include #include #include #include #include typedef CGAL::Gmpz bignum; using namespace std; using in = long long int; using PII = pair; using VI = vector; using VB = vector; #define sz(v) (v).size() #define forn(i,n) for(in i=0; i<(n); ++i) #define forv(i,v) forn(i,sz(v)) #define DEB(x) // Replace by "DEB(x) x" to get debug output using namespace std; bignum tob(in x) { return (bignum)(double)x; } // M(n,A): Unweighted Version vector> LFMemo; vector> PSLFMemo; // Prefix sums // D(n,d): Weighted Version vector> WLFMemo; vector> PSWLFMemo; // Prefix sums int globalw; const in MAXW = 110; void init() { LFMemo.resize(MAXW); PSLFMemo.resize(MAXW); WLFMemo.resize(MAXW); PSWLFMemo.resize(MAXW); } bignum LFcount(in A, in w, in l); bignum SumLFcount(in A, in w, in la, in lb) { // sum for l in [la,lb] if(la > lb) return 0; if(la < 0) return 0; if(la==lb) return LFcount(A,w,la); if(PSLFMemo[w].find({A,w,lb}) == PSLFMemo[w].end()) { PSLFMemo[w][{A,w,lb}] = SumLFcount(A,w,0,lb-1) + LFcount(A,w,lb); } return PSLFMemo[w][{A,w,lb}] - SumLFcount(A,w,0,la-1); } bignum LFcount(in A, in w, in l) { if(w>3) { LFMemo[w-3].clear(); PSLFMemo[w-3].clear(); } if(A==0 && w == 0 && l==0 ) { return 1; } if(A < 0 || w <= 0 || l < 0) return 0; if(A > w*w/2) return 0; if(LFMemo[w].find({A,w,l}) != LFMemo[w].end()) { return LFMemo[w][{A,w,l}]; } if(w + 5 < globalw) cout << "A=" << A << " w=" << w << " l=" << l << endl; // if(w>5) { // // cout << "del for w " << w << endl; // PSLFMemo[w-5].clear(); // LFMemo[w-5].clear(); // } bignum res = 0; res += SumLFcount(A-l, w-1, l, w/2); // for(in ll = l; ll<=w/2; ll++) { // res += LFcount(A-l, w-1, ll); // } if(l>0) { res += SumLFcount(A-(2*l-1), w-2, l-1, w/2); // for(in ll = l-1; ll<=w/2; ll++) { // res += LFcount(A-(2*l-1), w-2, ll); // } } LFMemo[w][{A,w,l}] = res; // cout << " computed " << "LFcount(" << A << "," << w << "," << l << ") = " << res << endl; return res; } bignum LFcount(in A, in w) { bignum res = 0; for(in l = 0; l <= w; l++) { bignum x = LFcount(A,w,l); // cout << "add " << "LFcount(" << A << "," << w << "," << l << ") =... " << x << endl; res += x; } return res; } bignum LFcount(in w) { bignum res = 0; for(in A = 0; A <= w*w; A++) { res += LFcount(A,w); } return res; } // Weighted Version bignum WLFcount(in A, in w, in l); bignum SumWLFcount(in A, in w, in la, in lb) { // sum for l in [la,lb] if(la > lb) return 0; if(la < 0) return 0; if(la==lb) return WLFcount(A,w,la); if(PSWLFMemo[w].find({A,w,lb}) == PSWLFMemo[w].end()) { PSWLFMemo[w][{A,w,lb}] = SumWLFcount(A,w,0,lb-1) + WLFcount(A,w,lb); } return PSWLFMemo[w][{A,w,lb}] - SumWLFcount(A,w,0,la-1); } bignum WLFcount(in A, in w, in l) { if(w>3) { LFMemo[w-3].clear(); PSLFMemo[w-3].clear(); } // cout << "A " << A << " w " << w << " l " << l << endl; if(A==0 && w == 0 && l==0 ) { // cout << "base" << endl; return 1; } if(A < 0 || w <= 0 || l < 0) return 0; if(A > w*w/2) return 0; if(WLFMemo[w].find({A,w,l}) != WLFMemo[w].end()) { return WLFMemo[w][{A,w,l}]; } bignum res = 0; res += tob(2*l+1)*SumWLFcount(A-l, w-1, l, w/2); if(l>0) { res += tob(l)*tob(l)*SumWLFcount(A-(2*l-1), w-2, l-1, w/2); } WLFMemo[w][{A,w,l}] = res; return res; } bignum WLFcount(in A, in w) { bignum res = 0; for(in l = 0; l <= w; l++) { bignum x = WLFcount(A,w,l); res += x; } return res; } bignum WLFcount(in w) { bignum res = 0; for(in A = 0; A <= w*w; A++) { res += WLFcount(A,w); } return res; } int main() { init(); in A,w,h; in WStop = 100; cout << "\nMotzkin Path Counting Demo\n\n"; cout << "Motzkin Numbers from 1 to 10\n"; cout << "see: https://en.wikipedia.org/wiki/Motzkin_number\n"; cout << "see: http://oeis.org/A001006\n"; for(in w = 1; w<=WStop; w++) { bignum res = LFcount(w); cout << "M(" << w << ") = " << res << endl; } cout << "\n\nM(n,A): Print Motzkin Path Counts by width and area using the last fall DP\n"; cout << "see: http://oeis.org/A129181\n"; for(w = 1; w<=WStop; w++) { for(A = 0; A <= w*w; A++) { bignum res = LFcount(A,w); // assert(res == weighted_count(A,w)); if(res>0) cout << res << "\t"; } cout << endl; } cout << "\n\nPrint Factorial Numbers w! = Sum of Weighted Motzkin Path Counts by width for a fixed w\n"; cout << "see: http://oeis.org/A000142\n"; for(w = 1; w<=100; w++) { globalw = w; bignum res = WLFcount(w); cout << w << "! = " << res << endl; } cout << "\n\nD(n,A): Print Weighted Motzkin Path Counts by width and area using the last fall DP\n"; cout << "see: https://oeis.org/A062869\n"; for(w = 1; w<=WStop; w++) { for(A = 0; A <= w*w; A++) { bignum res = WLFcount(A,w); // assert(res == weighted_count(A,w)); if(res>0) cout << res << "\t"; } cout << endl; } }