「ペア・ドロップ」問題

「ペア・ドロップ」問題

締め切りが 2017/12/15 10:00 AM なので，その 1 分後に投稿されるように予約

n を自然数とします。1 から n までの自然数が 1 つずつ書かれた n 枚のカードが 2 組あります。
 
これら 2n 枚のカードをよく混ぜ、A と B の 2 人に n 枚ずつ配ります。
A と B は、それぞれ自分の持ち札の中に番号が一致するカードがあればその 2 枚を捨てます。
 
捨てられるカードを全て捨てた後の A の持っているカードの枚数が k である確率を F(n, k) とします。
例えば F(1, 0)＝0, F(1, 1)＝1, F(2, 0)＝1/3, F(2, 1)＝0, F(2, 2)＝2/3, F(4, 2)＝24/35 です。
 
標準入力から、自然数 n と k （1 ≦ k ≦ n ≦ 100）が半角空白区切りで与えられます。
標準出力に 10^6×F(n, k) の整数部分の値を出力するプログラムを書いてください。
 
例えば、(n, k)＝(2, 2) であれば 666666 と出力してください。
(n, k)＝(4, 2) であれば 685714 と出力してください。
なお、全てのテストケースで、 10^6×F(n, k) の値とその最寄りの整数は 10^(-3) 以上離れていることが保証されています。
 
===================================

末尾に掲載するプログラムで，n, k の小さい部分を求める。

n について和を求めたものは，オンライン整数列大辞典の A000984 Central binomial coefficients a(n) = (2n)! / (n!)^2 であるこ。
更に，A005430（0, 2, 12, 60, 280, 1260, ...）の a(n-1) が Sum_{k=0..floor(n/2)} k*C(n,k)*C(n-k,k)*2^(n-2*k) であることより，以下の短いプログラムを得る。

 f = function(n, k) {
  k = n %/% 2 - k %/% 2
  cat(floor(1e6 * choose(n, k) * choose(n-k, k) * 2^(n-2*k) / choose(2*n, n)))
}

# arg = scan(file("stdin", "r"))
# f(arg[1], arg[2])

f(2, 2) # 666666
f(4, 2) # 685714
f(5, 1) # 238095
f(13, 9) # 211187
f(30, 20) # 67130
f(41, 17) # 126481
f(80, 24) # 226
f(99, 47) # 137249

＝＝ g(n, k) の分布＝＝＝＝

> g = function(n) {
+   a = combn(rep(1:n, 2), n)
+   a = apply(a, 2, function(x) {
+             b = table(x)
+             b = b[b == 1]
+             length(b)
+   })
+   table(a)
+ }
> g(1)
a
1
2
> g(2)
a
0 2
2 4
> g(3)
a
 1  3
12  8
> g(4)
a
 0  2  4
 6 48 16
> g(5)
a
  1   3   5
 60 160  32
> g(8)
a
   0    2    4    6    8
  70 2240 6720 3584  256