ARC064F: Rotated Palindromes

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";
}

まず場合の数パートが難しい。もっとサラッと考えられないかなぁ。