SDOI2018 反迴文串（莫比烏斯反演+Pollard-Rho）

題目連結

題目大意

求所有的串，滿足其所有迴圈同構串中至少有一個串是迴文串。
$n\le 10^{18}$

題解

第一步我就想偏了orz……我以為要分析這樣串的性質……
考慮所有迴文串，共有$k^{\left\lceil\frac{n}{2}\right\rceil}$個，我們考慮把每個迴文串第一個字元挪到後面，不斷進行這樣的操作，直到形成了新的迴文串為止。假設這樣操作了$t$次，那麼這個迴文串的貢獻就是$t$。
比如baab，一次變成aabb，第二次變成abba，是個迴文串，因此這個迴文串的貢獻是2.
考慮原串的最小整週期為$k$，顯然一個週期也是個迴文串。若$k$為奇數，那麼貢獻就是$k$，否則貢獻是$\frac k2$。
令$f(i)$表示週期為$i$的因數的所有迴文串個數，$g(i)$表示最小整週期恰好為$i$的所有迴文串個數。則有：
$f(i)=\sum_{d|i}g(d)=2^{\left\lceil\frac{i}{2}\right\rceil}$
再令$h(i)$表示最小整週期為$i$時的貢獻，即$i$為奇數是$h(i)=i$，否則$h(i)=\frac i2$。於是先使用莫比烏斯反演計算出$g$，就可以統計答案了。
$ans=\sum_{i|n}g(i)h(i)=\sum_{i|n}h(i)\sum_{d|i}\mu(d)f(\frac id)\\ =\sum_{i|n}f(i)\sum_{d|\frac ni}\mu(d)h(id)$
但是這樣計算還是會超時，考慮後面那個sigma的性質。$h(id)\neq d\cdot h(i)$當且僅當$i$為奇數且$d$為偶數。如果$d$是偶數，那麼$\frac ni$也必然是偶數。此時考慮任意的奇數$x$，必然有$\mu(x)h(ix)+\mu(2x)h(2ix)=0$，而如果後面sigma中的$d$含有4這個因數的話$\mu(d)=0$。因此：若$\frac ni$為偶數，$i$為奇數時後面的sigma取值為0.於是剩下來所有情況都是$h(id)=d\cdot h(i)$了。
$ans=\sum_{i|n,i~is~even~or \frac ni is ~odd}f(i)h(i)\sum_{d|\frac ni}\mu(d)d$
顯然後面那個sigma也是個積性函式，於是就可以得到：
$\sum_{d|i}\mu(d)d=\prod_{p|\frac ni,p~is~prime}(1-p)$
於是這題再用Pollard-Rho分解一下$n$就做完了。

#include <bits/stdc++.h>
namespace IOStream {
	const int MAXR = 1 << 23;
	char _READ_[MAXR], _PRINT_[MAXR];
	int _READ_POS_, _PRINT_POS_, _READ_LEN_;
	inline char readc() {
	#ifndef ONLINE_JUDGE
		return getchar();
	#endif
		if (!_READ_POS_) _READ_LEN_ = fread(_READ_, 1, MAXR, stdin);
		char c = _READ_[_READ_POS_++];
		if (_READ_POS_ == MAXR) _READ_POS_ = 0;
		if (_READ_POS_ > _READ_LEN_) return 0;
		return c;
	}
	template<typename T> inline void read(T &x) {
		x = 0; register int flag = 1, c;
		while (((c = readc()) < '0' || c > '9') && c != '-');
		if (c == '-') flag = -1; else x = c - '0';
		while ((c = readc()) >= '0' && c <= '9') x = x * 10 + c - '0';
		x *= flag;
	}
	template<typename T1, typename ...T2> inline void read(T1 &a, T2 &...x) {
		read(a), read(x...);
	}
	inline int reads(char *s) {
		register int len = 0, c;
		while (isspace(c = readc()) || !c);
		s[len++] = c;
		while (!isspace(c = readc()) && c) s[len++] = c;
		s[len] = 0;
		return len;
	}
	inline void ioflush() {
		fwrite(_PRINT_, 1, _PRINT_POS_, stdout), _PRINT_POS_ = 0;
		fflush(stdout);
	}
	inline void printc(char c) {
		_PRINT_[_PRINT_POS_++] = c;
		if (_PRINT_POS_ == MAXR) ioflush();
	}
	inline void prints(char *s) {
		for (int i = 0; s[i]; i++) printc(s[i]);
	}
	template<typename T> inline void print(T x, char c = '\n') {
		if (x < 0) printc('-'), x = -x;
		if (x) {
			static char sta[20];
			register int tp = 0;
			for (; x; x /= 10) sta[tp++] = x % 10 + '0';
			while (tp > 0) printc(sta[--tp]);
		} else printc('0');
		printc(c);
	}
	template<typename T1, typename ...T2> inline void print(T1 x, T2... y) {
		print(x, ' '), print(y...);
	}
}
using namespace IOStream;
using namespace std;
typedef long long ll;

const int B = 10, prime[B] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };
ll modmul(ll a, ll b, ll c) {
	return (ll)((__int128)a * b % c);
}
ll modpowi(ll a, ll b, int c) {
	ll res = 1;
	for (a %= c; b; b >>= 1) {
		if (b & 1) res = res * a % c;
		a = a * a % c;
	}
	return res;
}
ll modpow(ll a, ll b, ll c) {
	ll res = 1;
	for (a %= c; b; b >>= 1) {
		if (b & 1) res = modmul(res, a, c);
		a = modmul(a, a, c);
	}
	return res;
}
int miller_rabin(ll p) {
	for (int i = 0; i < B; i++) {
		if (prime[i] == p) return 1;
		else if (p % prime[i] == 0) return 0;
	}
	ll t = p - 1; int l = 0;
	for (; !(t & 1); t >>= 1) ++l;
	for (int i = 0; i < B; i++) {
		ll x = modpow(prime[i], t, p);
		for (int j = 0; j < l; j++) {
			ll y = modmul(x, x, p);
			if (y == 1 && x != 1 && x != p - 1) return 0;
			x = y;
		}
		if (x != 1) return 0;
	}
	return 1;
}
ll pollard_rho(ll n) {
	ll x, y = x = rand() % (n - 1) + 1, c = rand() % (n - 1) + 1, k = 1;
	for (int i = 1;; i++) {
		x = (modmul(x, x, n) + c) % n;
		if (x == y) return 1;
		ll g = __gcd(n, x < y ? y - x : x - y);
		if (g > 1 && g < n) return g;
		if (i == k) y = x, k <<= 1;
	}
}
map<ll, int> mp;
vector<pair<ll, int> > vec;
void fact(ll n) {
	if (miller_rabin(n)) { ++mp[n]; return; }
	ll p; while ((p = pollard_rho(n)) == 1);
	fact(p), fact(n / p);
}
ll ans, n, k; int T, p;
void dfs(int d, ll num, ll mul) {
	if (d < 0) {
		if ((num & 1) && !(n & 1)) return;
		ll h = (num & 1 ? num : num >> 1) % p;
		(ans += modpowi(k, (num + 1) >> 1, p) * mul % p * h) %= p;
		return;
	}
	ll t = 1, q = vec[d].first, e = vec[d].second;
	for (int i = 0; i <= e; i++, t *= q)
		dfs(d - 1, num * t, mul * (i < e ? 1 - q : 1) % p);
}
int main() {
	srand(time(0));
	for (scanf("%d", &T); T--;) {
		scanf("%lld%lld%d", &n, &k, &p);
		if (n == 1) { printf("%lld\n", k % p); continue; }
		mp.clear(), vec.clear(), ans = 0; fact(n);
		for (auto a : mp) vec.push_back(a);
		dfs(vec.size() - 1, 1, 1);
		printf("%lld\n", (ans + p) % p);
	}
	return 0;
}