https://arc064.contest.atcoder.jp/tasks/arc064_d
f(l)=周期がちょうどlの条件を満たす数列の数とすると、
f(l)=K^[(l+1)/2]-Σf(l/p)となって、求める値はΣ(lは奇数)f(l)*l+Σ(lは偶数)f(l)*(l/2)です。
f(l)の式は高速メビウス変換の式に良く似ています。
http://omochan.hatenablog.com/entry/2017/09/26/235331
ただ集合が0と1以外の値も取るので注意しましょう。
dp[k+1][a0a1....ak....am]=dp[k][a0a1....ak....am]-dp[k][a0a1....(ak-1)....am]のように遷移すれば良いです。
かなり高速でN<10^18でも通りそうです。多倍長必要ですが。
でもN=10^9程度だったら(lとしてありうるものの個数)<=2000らしいので、わざわざdpにして高速化しなくても2000*2000で通ります。
ll N, K, M; ll ans; map<int, int> P; vector<pi> vp; map<vi, ll> dp[11]; ll diviser(const vi& vec) { ll n = 1; rep(i, 0, M) { n *= mod_pow(vp[i].fst, vec[i]); } return n; } ll loop(vi vec, int at) { if(dp[at].find(vec) != dp[at].end()) return dp[at][vec]; if(at == 0) { ll l = diviser(vec); return dp[at][vec] = mod_pow(K, (l + 1) / 2); } else { int p = vec[at - 1]; ll res = 0; if(p - 1 >= 0) { vec[at - 1] = p - 1; ADD(res, mod - loop(vec, at - 1)); } vec[at - 1] = p; ADD(res, loop(vec, at - 1)); return dp[at][vec] = res; } } void loop2(vi vec, int at) { if(at == M) { ll l = diviser(vec); if(l % 2 == 0) ADD(ans, (l / 2) * loop(vec, M) % mod); else ADD(ans, l * loop(vec, M) % mod); return; } vec.resize(at + 1); rep(i, 0, vp[at].sec + 1) { vec[at] = i; loop2(vec, at + 1); } } void solve() { cin >> N >> K; ans = 0; ll n = N; for(int i = 2; i * i <= n; i++) { while(n % i == 0) { n /= i; P[i]++; } } if(n != 1) P[n]++; vp = vector<pi>(all(P)); M = sz(vp); loop2({}, 0); cout << ans << "\n"; }
まず場合の数パートが難しい。もっとサラッと考えられないかなぁ。