動的計画法

まずはこの問題を解くための DP の式を立ててみましょう。
問題文の通りに DP の式を立てると時間 \(\mathrm{O}(NM)\) ・空間\(\mathrm{O}(M)\) 程度の計算量で答えを求められる式が立ちますが、もちろん遅いです。こういう場合は

• \(f(x) :=\) 「あと \(x\) 倍までは掛けても操作を停止しない」という状態の時、サイコロを何回振れるかという期待値

とおいて式を立てましょう。すると、

\[ \begin{aligned} &f(0) = 0 \\ &f(x) = 1 + \frac{1}{N} \sum_{i=1}^N f(\lfloor x/i \rfloor) \\ \iff &f(x) = \frac{N}{N-1} \left(1 + \frac{1}{N}\sum_{i=2}^N f(\lfloor x/i \rfloor)\right) \end{aligned} \]

という式を得ることができて、答えは \(f(M)\) になります。
さて、上式を利用すればメモ化再帰により答えを求めることができますが、このアルゴリズムは(適切に計算すると)どれくらいの計算量になるでしょうか？

証明のためにいくつかの観察を列挙していきます。

観察 1

\(d\) を正整数とする。 \(\lfloor M / d\rfloor\) の値としてあり得るものは \(\mathrm{O}(\sqrt{M})\) 通りである。

これは ABC230-E でも登場した事実なので詳しい説明は省略しますが、\(L = \lfloor \sqrt{M} \rfloor, K = \lfloor M/(L+1) \rfloor\) とおいたときに、\(x = \lfloor M/d \rfloor\) を満たす正整数 \(d\) が存在する \(x\) は以下の \(2\) 種類で尽くされます。

• \(x=1,2,\dots,L\)
• \(x=\lfloor M/1\rfloor,\lfloor M/2\rfloor, \dots,\lfloor M/K\rfloor\)

この事実を利用すると、ABC230-E と同様に商ごとにまとめて計算するテクニックが使えます。すなわち、

• \(Q(x)\) : (\(q = \lfloor x / d\rfloor\) を満たす正整数 \(d\) が存在する \(n\) 以下の正整数 \(q\) からなる集合)
• \(g(x, q) : \) (サイコロの出目 \(d\) が \(\lfloor x/d\rfloor = q\) を満たす確率)

として

\[f(x) = \frac{N}{N-1} \left(1 + \sum_{q \in Q(x) \setminus \lbrace 1 \rbrace} g(x, q) f(q)\right)\]

と変形できるので、\(f(x)\) を \(\mathrm{O}(\sqrt{x})\) 個の \(f(q)\) の和で表すことができます。

観察 2

正整数 \(i,j,M\) の間に以下の事実が成り立つ。

\[\lfloor \lfloor M/i \rfloor / j \rfloor = \lfloor M / (ij) \rfloor\]

証明としては、\(M\) は \(0 \leq s \lt i, 0 \leq r \lt j\) を満たす \(s,r\) 、および整数 \(q\) を用いて \(M = (qj + r)i + s\) と一意に表せるので、この式を両辺の \(M\) に代入すれば示すことができます。もちろん、\(f(\lfloor \lfloor \lfloor M / i \rfloor / j \rfloor / k \rfloor) = f(\lfloor M /(ijk) \rfloor)\) のように 切り捨て除算が \(3\) 回以上の場合も同様に示すことができます。

観察 3

上記のメモ化再帰のアルゴリズム内部において、 \(f(x)\) を計算する必要がある \(x\) としてあり得るのは \(\lfloor M / d \rfloor = x\) を満たす正整数 \(d\) が存在する整数だけである。

これは観察 2 から従います。再帰の内部では \(f(M)\) から \(f(\lfloor M / i \rfloor)\) を呼び出して、そこから \(f(\lfloor \lfloor M / i \rfloor / j \rfloor)\) を呼び出して… という様子が見られますが、観察 \(2\) から引数はすべて \(\lfloor M / d \rfloor\) (\(d\) は登場する除数の積) と等しい値なのがわかります。

観察 1, 2, 3 を合わせると次のような事実を導けます。

観察 4

\(f(M)\) は \(\mathrm{O}(M^{3/4})\) で計算できる。

今までの観察から

• \(f(M)\) の計算に必要な \(f(x)\) の値は全部で \(\mathrm{O}(\sqrt{M})\) 状態である。
• \(f(x)\) は \(\mathrm{O}(\sqrt{x})\) 回の演算で計算することができる。

ということがわかっています。上記の事実から計算量は高々 \(\mathrm{O}(M)\) なのがわかりますが、実はさらに小さいです。
\(M\) の商としてあり得る \(x\) は\(L = \lfloor\sqrt{M}\rfloor\) とおくと

• \(x=1,2,\dots,L\) と
• \(x=\lfloor M/1\rfloor,\lfloor M/2\rfloor, \dots,\lfloor M/L\rfloor\)

の和集合になることが観察 1 で示されています。それぞれの部分に分けて計算量を求めると、

• \(x=1,2,\dots,L\) の部分： \(\mathrm{O}(\sqrt{L}) \times L = \mathrm{O}(M^{3/4})\)
• \(x=\lfloor M/1\rfloor,\lfloor M/2\rfloor, \dots,\lfloor M/L\rfloor\)　の部分：\(\mathrm{O}(\sum_{i=1}^L \sqrt{M/i}) = \mathrm{O}(\int_{t=1}^L \sqrt{M/t}) = \mathrm{O}(M^{3/4})\)

となりいずれも \(\mathrm{O}(M^{3/4})\) になることが確認できました。

以上の観察により、この問題は \(f(M)\) をメモ化に計算する時間計算量 \(\mathrm{O}(M^{3/4})\) ・空間計算量 \(\mathrm{O}(M^{1/2})\) のアルゴリズムで解くことができて、適切に実装すると制約下で十分高速に動作するので AC することができます。

• \(\mathrm{O}(M)\) 解法との切り分けの都合上 TL が厳しめになっていて、std::map 等を使用した解法は落ちると思います。ご了承ください。

C++ による実装例を以下に貼ります。

#include <cmath>
#include <iostream>
using namespace std;
#include "atcoder/modint.hpp"
using mint = atcoder::modint1000000007;

mint low[50000], high[50000];

int main() {
  int N, M;
  cin >> N >> M;
  int L = llround(sqrt(M)), K = M / (L + 1);
  for (int i = 1; i <= L + K; i++) {
    int x = i <= L ? i : M / (L + K + 1 - i);
    mint sum = 0;
    for (int l = 2, r, le = min(N, x); l <= le; l = r) {
      int q = x / l;
      r = min(N, x / q) + 1;
      mint cur = q <= L ? low[q] : high[M / q];
      sum += cur * (r - l);
    }
    (i <= L ? low[i] : high[L + K + 1 - i]) = (sum + N) / (N - 1);
  }
  cout << (M <= L ? low[M] : high[1]).val() << "\n";
}

• \(f(x)\) の値を std::unordered_map で管理すると少し楽ができます。(定数倍によっては TLE するのでご注意ください。)

#include <iostream>
#include <unordered_map>
using namespace std;
#include "atcoder/modint.hpp"
using mint = atcoder::modint1000000007;

mint solve(int N, int x) {
  static unordered_map<int, mint> memo;
  if (memo.count(x)) return memo[x];
  mint res = 0;
  for (int R = x, L, q; R; R = L) {
    q = x / R, L = x / (q + 1);
    if (q != x) res += solve(N, q) * (min(N, R) - min(N, L));
  }
  return memo[x] = (N + res) / (N - 1);
}

int main() {
  int N, M;
  cin >> N >> M;
  cout << solve(N, M).val() << "\n";
}

発展：乗法的関数の prefix sum

正整数 \(a, b\) が互いに素であるときに \(f(a)f(b) = f(ab)\) が成り立つような関数を 乗法的関数 と呼びます。
実はこの問題は乗法的関数の prefix sum \(\sum_{i=1}^N f(i)\) を高速に計算するときに出てくるアルゴリズムを参考にして出題しています。
詳しくは以下の文献を参考にしてください。

• maspyのHP, Dirichlet 積と、数論関数の累積和
• Codeforces, [Tutorial] Math note — Dirichlet convolution

上記のブログに載っているテクニックを利用すると、元の問題は \(\mathrm{\tilde{O}}(M^{2/3})\) で解くこともできます。

また、以下の問題は同様のテクニックを利用して解くことができます。

• Xmas Contest 2019 D : Sum of (-1)^f(n)
• SPOJ : DIVCNT2
• Library Checker : Sum of Totient Function