TopCoder SRM 607 Div1 Hard PulleyTautLine

問題

座標平面上にn個の滑車と2本の釘がある。
i番目の滑車は((i+1) * d, 0)の位置にあり、釘は(0, 0)および((n + 1) * d, 0)にある。
滑車の半径はrで、どの二つの滑車も接触しないことが保証されている。

原点の釘からもう一方の釘へ、ロープをたるまないように通す。
ロープは太さを無視でき、滑車に接触させたり何度でも巻きつけたりしてもよい。

ロープの通し方のうち、k番目に短いもの（ただし同じ長さでも違うルートを通る通し方は異なるものとして数える）の長さを求めよ。

制約条件

n≦50
r＜5億
2*r+1≦d≦10億
k≦10^18

方針

上下に関して対称なので、最初の滑車にロープは上側から接触するとしてよい。
ロープの動きを考えると、

隣の滑車へまっすぐ移る
隣の滑車へ斜めに移る
半回転巻きついて方向を変える

のどれかしかなく、それぞれの個数を決めればロープの長さも一意に定まることがわかる。

なので、まっすぐの移動の回数をi回、斜めの移動をj回、半回転巻きつきの回数をk回使って最後の滑車へ行く場合の数をそれぞれ効率的に求められれば、
短い順にこの動かし方を試していき、それぞれの場合の数をkから引いていって、0になったときの動かし方が求める答え。

であるが、巻きつきの回数は非常に大きくなりうるのでこれではダメ。
100万回とか1000万回とか巻きつきうるので。

じゃあどうしたらいいかというと、
巻きつきは後から巻くことにして、半回転は連続しないようにしてやる。
すると半回転は横移動の直後にしか行えないので、横移動の回数以下になり、
横移動の回数も、80回を超えれば、横移動だけで80C40通りくらいあり、10^18を十分に超えるので、
横移動、半回転を併せて80回ずつ以下くらいDPして十分間に合うことがわかる。

で、巻く回数を取り除いてしまうと、
短い順に巻き方を試すということが不可能になるので、いい方法を考えなければならない。
これは、先にロープの長さを決めて、その長さ以下の巻き方の個数を数えるという二分探索をすればいい。
そうすれば、i, j, kの回数を決めたら、巻ける回数も決まって、combinationで簡単に巻く回数を入れた場合の数も求まる。

巻ける回数の上限がわかったら、
→が横移動の操作で、
○を巻く操作とすれば
→→→→○○○○○と仕切り一個を入れて、それを並べ替える個数が場合の数。
仕切りは一番右の→よりも後に置くと考える。

→→→→○○○｜○○こんな感じで、仕切りよりも後の○は使わないと考えると、
(横移動＋巻く回数＋１）C（横移動）倍すれば巻く回数が入る。

更に、二分探索をするようにしたので、
DPをするときに横移動は一緒くたにしてしまって、二分探索のときだけ斜めの回数を全部試してコンビネーションをかければよい。
これでDPの計算量も減って余裕で通る。

ソースコード

//spaghetti sourceの幾何ライブラリを使用
//tanCPは一点を通る円の接線を求める関数（省略）
inline void add(ll &a, ll b){
	a += b;
	if(a + b > 1e18) a = (ll)1e18;
}
inline void mul(ll &a, ll b){
	if(1.0 * a * b > 1e18) a = (ll)1e18;
	else a *= b;
}
inline ll C(ll n, ll k){
	k = min(k, n - k);
	if(k > 50) return (ll)1e18;
	
	double d = 1;
	ll res = 1;
	rep(i, k){
		d *= n - i, d /= i + 1;
		res *= n - i, res /= i + 1;
	}
	if(d > 1e18) return (ll)1e18;
	if(abs(res - d) / d > 0.1) return (ll)d;
	return res;
}

ll dp[80][80][60][2]; //横移動, 回転, 位置, 直前に回ったか

class PulleyTautLine {
	public:
	double getLength(int d, int r, int n, long long k) {
		
		vector<P> ps1 = tanCP(P(d, 0), r, P(0, 0));
		vector<P> ps2 = tanCP(P(d / 2.0, 0), r, P(0, 0));
		double A = abs(ps1[1]) + (arg(ps1[1] - P(d, 0)) - PI / 2) * r;
		double B = abs(ps2[1]) + (arg(ps2[1] - P(d / 2.0, 0)) - PI / 2) * r;
		double D = 2 * PI * r;
		B *= 2;
		
		k = (k + 1) / 2;
		if(n == 1) return 2 * A + D * (k - 1);
		
		memset(dp, 0, sizeof(dp));
		dp[0][0][0][0] = 1;
		
		rep(i, 80) rep(j, 80) rep(pos, n) rep(f, 2){
			int cur = dp[i][j][pos][f];
			if(!cur) continue;
			//回転は連続しないようにする
			if(j + 1 < 80 && !f) add(dp[i][j + 1][pos][1], cur);
			
			int npos = j % 2 ? pos - 1 : pos + 1;
			if(i + 1 < 80 && 0 <= npos && npos < n) add(dp[i + 1][j][npos][0], cur);
		}
		
		double lo = 0, hi = 1e20;
		rep(it, 100){
			double mid = (lo + hi) / 2;
			ll sum = 0;
			
			rep(i, 80) rep(j, 80){
				ll way = dp[i][j][n - 1][0];
				if(!way) continue;
				
				rep(l, i + 1){
					double len = 1.0 * l * d + (i - l) * B + D / 2 * j;
					if(len > mid) continue;
					
					ll t = (ll)(min((mid - len) / D, 1e18));
					ll w = C(i, l);
					mul(w, C(i + t + 1, t));
					mul(w, way);
					
					add(sum, w);
					if(sum >= k) goto SKIP;
				}
			}
			SKIP:
			if(sum >= k) hi = mid;
			else lo = mid;
		}
		return hi + 2 * A;
	}
};