Makoto and a Blackboard

$p_i (v) \coloneqq$「$i$回目で数字が$v$となる確率」とすると，漸化式$p_{i+1}(v) \leftarrow \sum_{w \text{ s.t. } v \in D(w)} p_{i} (w) / |D(w)| \ (D(v) \coloneqq v$の約数$)$が立つ．これに従うと求める値は$\sum _{w \in D(n)} w p _ k (w)$と書けるが，これを求めるには$n$の約数を再帰的に因数分解する必要がありそうだし，それを除いても$O(|D(n)|k)$は掛かりそう．約数の個数$|D(n)|$の評価は難しいようだが少なくとも間に合わない．ここで詰まった．

ここで「$v$の約数$d$を選ぶ」というのは，$v= \prod p_i ^ {k _ i}$と素因数分解したとき，$d = \prod p_i ^ {k'_i} \text{ s.t. } k' _ i \le k _ i$を選ぶ，すなわち各素因数$p _ i$の個数$k' _ i$を決定することである． 例えば$v=72=2^3 \times 3^2$のとき，約数は$\{1,2,3,4,6,8,9,12,18,24,36,72\}$であるが，これは$2^3$の約数$\{1,2,4,8\}$と$3^2$の約数$\{1,3,9\}$の積である．このように約数を選ぶことは各素因数の組からそれぞれ好きなものを選ぶことと等しい． この決定は互いに独立であるから，「積の期待値＝期待値の積」となり，$E[v]=E[\prod p _ i ^ {k _ i}] = \prod E[p _ i ^ {k _ i}]$という風に各素因数に対して求めた結果の積を取ればよい^[1]．

以下では素因数$p$について求める．いま初期値$n$に含まれる$p$が$l$個のとき，$dp[i][j] := i$回目で数字が$p^j$となる確率とおく．すると素朴には$dp[i+1][j] \gets \sum_{j \le w \le l} dp[i][w]/(j+1)$と更新できる．この計算量は$O(kl^2)$であり，$l$は最大でも$O(\log_2 n)$なので$O(k \log ^2 n)$．

実際はこれだと間に合わない（手元やAtCoderだと1sだがCodeforces上だと3s．C++なら間に合うかも？）ため，DPを高速化する必要がある．これは$j$を降順に試すとき，$dp[i+1][j] \gets dp[i+1][j+1]+dp[i][j]/(j+1)$であることに着目すると$w$が不要になることを利用する^[2]．これにより計算量は$O(k \log n)$となる．

これを各素因数について行うと計算量は$O(L k \log n)$ $(L$ $\coloneqq$素因数の個数$)$と書ける．素数を小さい順に14個掛け合わせると約$1.3 \times 10^{16} \gt n$であるから$L \lt 14$であり間に合いそうな気がしてくる．以上から最初の素因数分解とあわせて全体の計算量は$O(\sqrt n + k \log n)$．

初めての解答方式だったが，これは逆元を求めればよい．最初は「分数$P/Q$で保持して最後に$P \times inv(Q)$しなければならない」と思っていたのだが，実際は分数$p/q$が出てくるたびに$p \times inv(q)$として扱って大丈夫．逆元の実装はeditorialにもあるFermatの小定理版だとかなり低速だったので，拡張ユークリッド法によるものを用いた（初実装．参考記事）．因みに他の人の提出を見ていると，modつき$_nC_r$を計算するときのように配列に前計算しているケースも多かった．

ACコード（TLEギリギリ）: 解答例

追記（19/03/07）：素因数分解でリストを作っている箇所を大きさを決め打った配列に直したら560ms速くなった．あるいは大きさを推定せずともHashtblを用いてもよさそう（手元で5回測定した平均だとテストケース"10^15 10^4"で前者のほうが10ms早い程度）．今後はこのような処理ではHashtblを用いたほうがよいのかもしれない？

約数を素因数毎に分解して考えるというのは盲点だった．一方で約数の個数を求める式も似た考えで導出できることから，そんなに突飛な考え方でもなかったのかもしれない．

---

1. 記法は適当． ↩︎

2. 加えて速度にはそこまで影響しなそうだが$j$を降順にすることで$i$の番地も不要になるため，解答ではそのようにしている． ↩︎