公式解説を少し補足したいと思います。
あと、コードは書きません。
問題のリンク
step 1
1~mのサイコロを振ったとき、表がちょうど$N$回出る回数については、
$$\left[x^N\right]\prod_{i=1}^{m}\left(\frac{i^2-1}{i^2}+\frac{x}{i^2}\right)^A$$
$$=\left[x^N\right]\prod_{i=1}^{m}\left(1-\frac{1-x}{i^2}\right)^A$$
$$=\left[x^N\right]\prod_{i=1}^{m}\left(1-\frac{\sqrt{1-x}^2}{i^2}\right)^A$$
となっている。$m\to\infty$での極限を飛ばすと
$$=\left[x^N\right]\prod_{i=1}^{\infty}\left(1-\frac{\sqrt{1-x}^2}{i^2}\right)^A$$
となる。ここで知識があると以下のことを思い浮かべることができる。
$$\frac{\sin{\pi z}}{\pi z}=\prod_{n=1}^{\infty}\left(1-\frac{z^2}{n^2}\right)$$
すると、結局以下の通りになる。
$$\left[x^N\right]\left(\frac{\sin{(\pi\sqrt{1-x})}}{\pi\sqrt{1-x}}\right)^A$$
step2
テイラー展開を考えると、以下の通りになる。
$$N! \frac{d^N}{dx^N}\left.\left(\frac{\sin{(\pi\sqrt{1-x})}}{\pi\sqrt{1-x}}\right)^A\right|_{x=0}$$
$$=N! \frac{d^N}{dx^N}\left.\left(\frac{\sin{(\sqrt{-x\pi^2+\pi^2})}}{\sqrt{-x\pi^2+\pi^2}}\right)^A\right|_{x=0}$$
$$=N!{(-\pi^2)}^N\frac{d^N}{dx^N}\left.\left(\frac{\sin{(\sqrt{x})}}{\sqrt{x}}\right)^A\right|_{x=\pi^2}$$
step3
$f_k(x)=\frac{d^k}{dx^k}\sin^A{x}$とおく。すると以下のように書けるらしい。
$$\mathrm{ans}=\frac{N!{(-\pi^2)}^N}{2^N}\sum_{k=0}^{N}(-1)^k\cdot c_{A,k}\cdot f_{N-k}(x)\cdot (\sqrt{x})^{-(A+N+k)}$$
$c_{1,i}$を表に書くと以下のようになる。注意点として、$c_{1,i}$は$N$に依存する。
1 | |||||||
1 | 1 | ||||||
1 | 3 | 3 | |||||
1 | 6 | 15 | 15 | ||||
1 | 10 | 45 | 105 | 105 | |||
1 | 15 | 105 | 420 | 945 | 945 | ||
1 | 21 | 210 | 1260 | 4725 | 10395 | 10395 | |
1 | 28 | 378 | 3150 | 17325 | 62370 | 135135 | 135135 |
OEISに投げてみる。とりあえず対角線の部分(1,1,3,15,105,945,10395,135135)
対角線2列目も投げてみる。(1,6,45,420,4725,62370)
対角線3列目(1,10,105,1260,17325)
なんとなく以下の数式が見える。対角線$i$列目は以下のようになりそうである。
$$\frac{(2N+i)!}{2^Ni!N!}$$
ズレをなくすと、結局以下の通りになる。
$$c_{1,i}=\frac{(N + i)!}{(N - i)!i!2^i}$$
step4
解説によると、どうやら以下の漸化式が成り立つようである。
$$c_{A,i}=c_{A-1,i}+(N-i+1)c_{A,i-1}$$
これだと扱いにくいので、両辺を$(N-i)!$倍する。すると
$$(N-i)!c_{A,i}=(N-i)!c_{A-1,i}+(N-i+1)!c_{A,i-1}$$
$b_{A,i}=(N-i)!c_{A,i}$とすると、
$$b_{A,i}=b_{A-1,i}+b_{A,i-1}$$
つまり、
$$b_{A,i}=b_{A-1,i}+b_{A-1,i-1}+\cdots+b_{A-1,0}$$
となる。
すると、以下のようになる。
$$b_{A,i}=\left[x^i\right] \left(\sum_{k=0}^{\infty}b_{1,k}x^k\right)\left(\frac{1}{1-x}\right)^{A-1}$$
これは形式的べき級数のライブラリを使うことで$O(N(\log{N})^2)$または$O(N\log{N})$で解くことができる。
step 5
$f_i(x)$についても求める。
$$\left.\frac{d^n}{dx^n}\sin{x}\right|_{x=\pi}=\sin{\left(\frac{\pi}{2}n+\pi\right)}$$
である。
求め方としては、以下を使えば良いだろう。
$$f_K(\pi)=\sum_{k_1+\cdots+k_A=K}\frac{K!}{k_1!\cdots k_A!}\sin{(\pi k_1/2+\pi)}\cdots\sin{(\pi k_A/2+\pi)}$$
ここで、$k_1,\ldots,k_A$が全部奇数の場合のみ考えればよい。よって、
$k_1=2i_1+1,\ldots,k_A=2i_A+1$とすると、
$$f_K(\pi)=\sum_{i_1+\cdots+i_A=\frac{K-A}{2}}\frac{K!}{(2i_1+1)!\cdots (2i_A+1)!}{(-1)}^{i_1+1}\cdots {(-1)}^{i_A+1}$$
$$f_K(\pi)=\sum_{i_1+\cdots+i_A=\frac{K-A}{2}}\frac{K!}{(2i_1+1)!\cdots (2i_A+1)!}(-1)^{\frac{K+A}{2}}$$
$$f_K(\pi)=K!(-1)^{\frac{K+A}{2}}\left[x^{\frac{K-A}{2}}\right]\left(\sum_{i=0}^{\infty}\frac{1}{(2i+1)!}x^i\right)^A$$
これも多項式の累乗を使うと計算することができる。
また、$(K-A)$が奇数になる場合と負になる場合は、係数は0になる。
step 5'
(2021年1月7日追記)
これを書いたあとで気づいたことだが、実はもっと簡単に計算することができる。
$$f_K(\pi)=\left[x^K\right]K!\cdot \sin^A{(x-\pi)}=\left[x^K\right](-1)^A K!\cdot \sin^A{x}$$
ここで、$\sin^A{x}$のテイラー展開については、$\sin{x}$のテイラー展開について、形式的べき級数の積公式を使えば求めることができる。
この方法のメリットとしては、$K-A$の偶奇や符号を気にする必要がないことである。
実装上の注意点
pypyみたいな遅めの言語では、$O(N(\log{N})^2)$でもTLEになる場合があるので、形式的べき級数のexp,logを計算するライブラリを使うことで$O(N\log{N})$のプログラムを書くと良い。制限時間は5秒でギリギリ通るかもしれない。
感想
WTFの3141点枠で出そうだなと思った。