// Motzkin Path Enumerator
// It can count, enumerate, sample and print all Motzkin paths of a given area, width and height.
// Just compile and run it* to get a quick demo or look at the main function for some examples.
// g++ motzkin.cpp -o m -std=c++11 -O2; ./m

#include <iostream>
#include <vector>
#include <map>
#include <algorithm>
#include <cassert>

using namespace std;
using in = long long int;
using PII = pair<in,in>;
using VI = vector<in>;

#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;

in nCr(in n, in k) { // n choose k
	DEB(cout << "n=" << n << " choose k=" << k;)
	in res = 1;
	for(in i=1; i<=k; i++) {
		res *= n-i+1;
		res /= i;
	}
	DEB(cout << " res=" << res << endl;)
	return res;
}

in exp(in a, in b) { // comput a ^ b
	in r = 1;
	forn(i,b) r*= a;
	return r;
}

VI bitmask(in n, in k, in x) {
	DEB(cout << "bitmask n=" << n << " k=" << k << " x=" << x << ": ";)
	VI b(n,0);
	// Place the bits from right to left
	in l = 0;
	for(in i=0; i<k; i++) {
		in sum = 0;
		for(in j=n-1; j>=l; j--) {
			in c = nCr(n-j-1, k-i-1);
			if(sum+c > x) {
				b[j] = 1;
				l = j;
				x -= sum;
				break;
			}
			sum+=c;
		}
	}
	DEB(forv(i,b) cout << b[i];
	cout << endl;)
	return b;
}

struct MotzkinPath {
	VI deltas;
	VI profile;

	void read() {
		string s;
		getline(cin,s);
		deltas = VI(sz(s));
		forv(i,s) {
			deltas[i] = (s[i]=='U') ? 1 : (s[i]=='D' ? -1 : 0);
			DEB(cout << deltas[i];)
		}
		DEB(cout << endl;)
	}
	void set_profile() {
		profile = VI(sz(deltas)+1);
		profile[0] = 0;
		forv(i,deltas) {
			profile[i+1] = deltas[i] + profile[i];
			DEB(cout << profile[i];)
		}
		DEB(cout << endl;)
	}

	in height() {
		in maxh=0;
		forv(i,deltas) maxh = max(maxh,profile[i]);
		return maxh;
	}

	void print() {
		forn(i,3+sz(deltas)) cout << "~"; cout << endl;
		in h = height();
		for(in i=h; i>=0; i--) {
			cout << i << " ";
			forv(j,deltas) {
				in jh = min(profile[j],profile[j+1]);
				if(jh==i) {
					if(deltas[j]==1) cout << "/";
					if(deltas[j]==0) cout << "_";
					if(deltas[j]==-1) cout << "\\";
				} else {
					if(jh>i) {
						cout << "X";
					} else {
						cout << " ";
					}
				}
			}
			cout << endl;
		}
		// forn(i,3+sz(deltas)) cout << "~"; cout << endl;
	}

	// Add some new flats among the peaks
	/*             __   _
	   /\_/\/\ -> /  \_/ \/\
	*/
	void add_flats_on_top(VI bitmask) {
		DEB(cout << "add flats on top: ";
		forv(i,bitmask) cout << bitmask[i];
		cout << endl;)
		set_profile();
		VI new_deltas;
		in h = height();
		in j = 0;
		forv(i,profile) {
			if(profile[i]==h) {
				while(j<sz(bitmask) && bitmask[j]==0) {
					new_deltas.push_back(0);
					j++;
				}
				j++;
			}
			if(i<sz(deltas)) new_deltas.push_back(deltas[i]);
		}
		deltas = new_deltas;
		set_profile();
	}
	// Replace some of the top flats by peaks
	/*  __   _        _/\   /\
	   /  \_/ \/\ -> /   \_/  \/\
	*/
	void add_peaks_on_top(VI bitmask) {
		set_profile();
		VI new_deltas;
		in h = height();
		in j = 0;
		forv(i,deltas) {
			if(profile[i]==h && deltas[i]==0) {
				if(bitmask[j]==1) { // Replace by a peak
					new_deltas.push_back(1);
					new_deltas.push_back(-1);
				} else { // Keep the flat
					new_deltas.push_back(0);
				}
				j++;
			} else {
				new_deltas.push_back(deltas[i]);
			}
		}
		deltas = new_deltas;
		set_profile();
	}
};

map<VI,in> Memo;

in count(in A, in w, in h, in p) {
	DEB(cout << "call A=" << A << " w=" << w << " h=" << h << " p=" << p << endl;)
	if(A==0 && h == 0 && p==0 && w == 0) {
		DEB(cout << "return base case" << endl;)
		return 1;
	}
	if(A < 0 || w <= 0 || h <= 0) return 0;
	if(Memo.find({A,w,h,p}) != Memo.end()) {
		return Memo[{A,w,h,p}];
	}
	DEB(cout << "count A=" << A << " w=" << w << " h=" << h << " p=" << p << endl;)

	in outer_sum = 0;
	for(in f=0; f<=w; f++) {
		DEB(cout << " h=" << h << " flats: " << f << endl;)
		in outer_bin = nCr(p+f, f);
		in inner_sum = 0;
		for(in pp=(h==1 ? 0 : 1); pp<=(h==1 ? 0 : w); pp++) {
			DEB(cout << " h=" << h << " inner peaks: " << pp << endl;)
			in inner_bin = pp==0 ? 1 : nCr(p+f+pp-1, pp-1);
			in inner_dp = count(A - p*(2*h-1)-f*(h-1),
								w - 2*p - f,
								h-1,
								pp);
			inner_sum += inner_bin * inner_dp;
		}
		DEB(cout << "outer_bin=" << outer_bin << " * inner_sum=" << inner_sum << endl;)
		outer_sum += outer_bin * inner_sum;
	}
	Memo[{A,w,h,p}] = outer_sum;
	return outer_sum;
}

// Sum out number of peaks on the top
in count(in A, in w, in h) {
	return count(A,w,h+1,0);
}

in count(in A, in w) {
	in result = 0;
	for(in h=0; h<=w; h++) {
		result += count(A,w,h);
	}
	return result;
}

in count(in w) {
	in result = 0;
	for(in A=0; A<=w*w; A++) {
		result += count(A,w);
	}
	return result;	
}


// Get the x.th Motzkin Path of area A, width w, height h and top peak count p
MotzkinPath sample(in A, in w, in h, in p, in x) {
	DEB(cout << "sample A=" << A << " w=" << w << " h=" << h << " p=" << p << " x=" << x << endl;)
	if(A==0 && h == 0 && p==0 && w == 0) {
		return MotzkinPath(); // Empty path
	}
	// Sample the number of flats
	in sum = 0;
	for(in f=0; f<=w; f++) {
		in outer_bin = nCr(p+f, f);
		in inner_sum = 0;
		for(in pp=(h==1 ? 0 : 1); pp<=(h==1 ? 0 : w); pp++) {
			DEB(cout << " h=" << h << " inner peaks: " << pp << endl;)
			in inner_bin = pp==0 ? 1 : nCr(p+f+pp-1, pp-1);
			in inner_dp = count(A - p*(2*h-1) - f*(h-1),
								w - 2*p - f,
								h - 1,
								pp);
			in c = outer_bin * inner_bin * inner_dp;
			if(sum+c > x) {
				// We found f and pp
				DEB(cout << "found f=" << f << " and pp=" << pp << endl;)
				in inner_x = (x-sum) % (inner_dp);
				DEB(cout << "inner_x = " << inner_x << " = " << x-sum << "%%" << inner_dp << endl;)
				in outer_x = (x-sum) / (inner_dp);
				in pf_x = outer_x % inner_bin; // Iterator for the peak+flat positions
				in f_x = outer_x / inner_bin; // Iterator for the flats
				MotzkinPath path = sample(A - p*(2*h-1) - f*(h-1),
											w - 2*p - f,
											h - 1,
											pp,
											inner_x);
				VI flats_pattern = bitmask(p+f+pp-1, pp-1,pf_x);
				if(pp==0) flats_pattern = bitmask(p+f,0,0);
				path.add_flats_on_top(flats_pattern);
				VI peaks_pattern = bitmask(p+f, p, f_x);
				path.add_peaks_on_top(peaks_pattern);
				return path;
			}
			sum += c;
		}
	}
	cerr << "Sample path " << x << " out of " << sum << " is out of range!" << endl;
	return MotzkinPath();
}

// Get the x.th Motzkin Path of area A, width w and height h
MotzkinPath sample(in A, in w, in h, in x) { 
	return sample(A,w,h+1,0,x);
}

// Get the x.th Motzkin Path of area A and width w
MotzkinPath sample(in A, in w, in x) { 
	// decide on the height
	in sum = 0;
	for(in h=0; h<=w; h++) {
		in c = count(A,w,h);
		if(sum+c > x) {
			return sample(A,w,h+1,0,x-sum);
		}
		sum += c;
	}
	cerr << "Sample path " << x << " out of " << sum << " is out of range!" << endl;
	return MotzkinPath();
}

// Get the x.th Motzkin Path of width w
MotzkinPath sample(in w, in x) { 
	// decide on the area
	in sum = 0;
	for(in A=0; A<=w*w; A++) {
		in c = count(A,w);
		if(sum+c > x) {
			return sample(A,w,x-sum);
		}
		sum += c;
	}
	cerr << "Sample path " << x << " out of " << sum << " is out of range!" << endl;
	return MotzkinPath();
}

MotzkinPath random_sample(in A, in w, in h, in p) { 
	return sample(A,w,h,p,rand()%count(A,w,h,p));
}

MotzkinPath random_sample(in A, in w, in h) { 
	return sample(A,w,h,rand()%count(A,w,h));
}

MotzkinPath random_sample(in A, in w) { 
	return sample(A,w,rand()%count(A,w));
}

MotzkinPath random_sample(in w) { 
	return sample(w,rand()%count(w));
}

map<VI,in> WeightedMemo;

in weighted_count(in A, in w, in h, in p) {
	DEB(cout << "call A=" << A << " w=" << w << " h=" << h << " p=" << p << endl;)
	if(A==0 && h == 0 && p==0 && w == 0) {
		DEB(cout << "return base case" << endl;)
		return 1;
	}
	if(A < 0 || w <= 0 || h <= 0) return 0;
	if(WeightedMemo.find({A,w,h,p}) != WeightedMemo.end()) {
		return WeightedMemo[{A,w,h,p}];
	}
	DEB(cout << "weighted count A=" << A << " w=" << w << " h=" << h << " p=" << p << endl;)

	in outer_sum = 0;
	for(in f=0; f<=w; f++) {
		DEB(cout << " h=" << h << " flats: " << f << endl;)
		in outer_bin = nCr(p+f, f);
		in inner_sum = 0;
		in inner_exp = exp(2*h-1,f);
		for(in pp=(h==1 ? 0 : 1); pp<=(h==1 ? 0 : w); pp++) {
			DEB(cout << " h=" << h << " inner peaks: " << pp << endl;)
			in inner_bin = pp==0 ? 1 : nCr(p+f+pp-1, pp-1);
			in inner_dp = weighted_count(A - p*(2*h-1)-f*(h-1),
								w - 2*p - f,
								h-1,
								pp);
			inner_sum += inner_bin * inner_dp;
		}
		DEB(cout << "outer_bin=" << outer_bin << " * inner_sum=" << inner_sum << endl;)
		outer_sum += inner_exp * outer_bin * inner_sum;
	}
	in outer_exp = exp(h,2*p);
	// cout << "outer exp h=" << h << " p=" << p << " -> " << outer_exp << endl;
	in result = outer_exp * outer_sum;
	WeightedMemo[{A,w,h,p}] = result;
	return result;
}

// Sum out number of peaks on the top
in weighted_count(in A, in w, in h) {
	return weighted_count(A,w,h+1,0);
}

in weighted_count(in A, in w) {
	in result = 0;
	for(in h=0; h<=w; h++) {
		result += weighted_count(A,w,h);
	}
	return result;
}

in weighted_count(in w) {
	in result = 0;
	for(in A=0; A<=w*w; A++) {
		result += weighted_count(A,w);
	}
	return result;	
}

int main() {
	in A,w,h;

	cout << "\nMotzkin Path Enumerator Demo\n\n";

	cout << "Motzkin Numbers from 1 to 10\n";
	cout << "see: https://en.wikipedia.org/wiki/Motzkin_number\n";
	for(in w = 1; w<=10; w++) {
		in res = count(w);
		cout << "M(" << w << ") = " << res << endl;
	}

	w = 4;
	cout << "\n\nList all Motzkin Paths of width=" << w << endl;
	forn(x,count(w)) {
		MotzkinPath M = sample(w,x);
		M.print();
	}

	A = 9;
	w = 7;
	h = 2;
	cout << "\n\nList all Motzkin Paths of area=" << A << " width=" << w << " height=" << h << endl;
	forn(x,count(A,w,h)) {
		MotzkinPath M = sample(A,w,h,x);
		M.print();
	}

	A = 14;
	w = 9;
	h = 3;
	srand(time(0));
	cout << "\n\nA random Motzkin Path of area=" << A << " width=" << w << " height=" << h << endl;
	MotzkinPath M = random_sample(A,w,h);
	M.print();

	cout << "\n\nWeighted Motzkin Path Counts = n! from 1 to 10\n";
	for(in w = 1; w<=10; w++) {
		in res = weighted_count(w);
		cout << w << "! = " << res << endl;
	}

	cout << "\n\nWeighted Motzkin Path Counts = n! from 1 to 10\n";
	w = 4;
	for(A = 0; A <= 16; A++) {
		in res = weighted_count(A, w);
		cout << "weighted_count(A=" << A << ", w=" << w << ") = " << res << endl;
	}

	cout << "\n\nPrint Weighted Motzkin Path Counts as in Table 1 of the paper\n";
	for(w = 1; w<=10; w++) {
		for(A = 0; A <= w*w; A++) {
			in res = weighted_count(A,w);
			if(res>0) cout << res << "\t";
		}
		cout << endl;
	}

}