好厉害的题！为了方便描述，以下认为题目中给的 $k$ 为 $m$。

查询有多少个串满足它本身或循环位移为回文串。首先一个串不同的循环位移数量由它的最短周期决定，其次一个回文串重复出现仍然是回文串，这都是很好理解的。基于此容易想到做一个简单的莫反，设 $f(x)$ 为存在周期长度为 $x$ 的回文串数量，$g(x)$ 的最短周期长度恰好为 $x$ 的数量。$f(x) = \sum\limits_{d|x}g(d)$ 即可得到 $g(x) = \sum\limits_{d|x}\mu(d)f(\dfrac{x}{d})$。现在只要求解比较容易的 $f(x)$ 即可。

直接的想法是 $f(x)$ 构造出回文串的一半，这部分是 $m^{\lfloor\frac{x+1}{2}\rfloor}$ 的。最终答案是 $\sum\limits_{x | n}xg(x)$，但是写一下发现并不对。因为可能存在循环同构的两个回文串，但是发现这种情况的充要条件为这个回文串长度为偶数而且这两个串形如 $A + \operatorname{rev}(A)$，$\operatorname{rev}(A) + A$。手玩小情况就可以理解，因为如果不是恰好这么循环位移那么就会“传递”两个字符相等，最后会得到全部字符相等，但是要求的是最小周期为 $x$。而最小周期为 $1$ 也不会产生这种情况。于是最终答案应为：

$$\sum\limits_{x|n}\dfrac{x}{2 - x\bmod 2}\sum\limits_{d | x}\mu(d)m^{\lfloor\frac{\frac{x}{d}+1}{2}\rfloor}$$

考虑如何加快计算这个式子。形式很猎奇，而且 $m^{\lfloor\frac{\frac{x}{d}+1}{2}\rfloor}$ 并不是积性函数，但是 $x \bmod 2$ 似乎可以通过分类讨论进行处理，于是我们尝试把 $\mu(d)$ 和 $\dfrac{x}{2 - x\bmod 2}$ 合并到一起，具体而言，为了方便设 $F(x) = m^{\lfloor\frac{\frac{x}{d}+1}{2}\rfloor}, G(x) = \dfrac{x}{2 - x\bmod 2}$，那么：

$$\sum\limits_{x|n}G(x)\sum\limits_{d|x}\mu(d)F(\dfrac{x}{d}) = \sum\limits_{d|n}F(d)\sum\limits_{k|\frac{n}{d}}G(kd)\mu(k)\\ = \sum\limits_{d|n}m^{\lfloor\frac{\frac{x}{d}+1}{2}\rfloor}d\sum\limits_{k|\frac{n}{d}}\dfrac{k}{2 - (kd\bmod 2)}\mu(k)$$

如果遮住后面 $kd \bmod 2$ 那么就是很漂亮的 $k\mu(k)$ 了。考虑对这些奇偶性讨论。

• $d$ 为偶数。
  • 后面这部分为 $\dfrac{1}{2}\sum\limits_{k|\frac{n}{d}}k\mu(k)$。
• $d$ 为奇数。

  • 如果 $\dfrac{n}{d}$ 为偶数。
    此时对于 $k$ 为偶数的情况，那么就是 $\sum\limits_{k|\frac{n}{2d}}\dfrac{2k(-\mu(k))}{2} = -\sum\limits_{k|\frac{n}{2d}}k\mu(k)$。对于 $k$ 为奇数的情况，就是 $\sum\limits_{k|\frac{n}{2d}}k\mu(k)$。我们惊喜的发现这两部分之和居然是 $0$！

  • 如果 $\dfrac{n}{d}$ 是奇数。 此时就是 $\sum\limits_{k|\frac{n}{d}}k\mu(k)$。

问题就变成了对于 $n$ 的所有因数 $d$ 求解 $\sum\limits_{k|d}k\mu(k)$。

因为和 $\mu(k)$ 相关，所以只需要保留 $d$ 的质因数存在集合即可。而在 $10^{18}$ 内最多的不同质因数只有 $18$ 个。于是可以想到用集合表示这个 $d$，那么这个求和就是一个高维前缀和。大质因数分解要用 Pollard-Rho。时间复杂度为 $O(d(n)\log n + C2^C)$，$d(n)$ 为 $n$ 约数数量。由于 $d(n)$ 最多是 $10^5$ 级别的，所以可以通过。

#include <bits/stdc++.h>
#define ll long long
#define i128 __int128
using namespace std;
const int N = 350;
const int prime[15] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
const int NN = 18;
ll n, m, Mod;
void upd(ll &x, ll y) {
	x = ((x + y >= Mod) ? (x + y - Mod) : (x + y));
}
ll modint(ll x, ll y) {
	return ((x + y >= Mod) ? (x + y - Mod) : (x + y));
}
ll arr[(1 << NN) + 10];
void fwt(ll *arr, int n) {
	for(int k = 1; k < (1 << n); k <<= 1)
		for(int i = 0; i < (1 << n); i += (k << 1))
			for(int j = 0; j < k; j++)
				upd(arr[i + j + k], arr[i + j]);
}
ll qpow(ll n, ll m) {
	ll res = 1;
	while(m) {
		if(m & 1ll) res = (i128)res * (i128)n % Mod;
		n = (i128)n * (i128)n % Mod;
		m >>= 1;
	}
	return res;
}
ll gcd(ll n, ll m) {
	if(!m) return n;
	else return gcd(m, n % m);
}
ll qmul(ll n, ll m, ll Mod) {
	ll res = 0;
	while(m) {
		if(m & 1ll) res = (res + n) % Mod;
		m >>= 1ll;
		n = (n + n) % Mod;
	}
	return res;
}
ll qpow(ll n, ll m, ll Mod) {
	ll res = 1;
	while(m) {
		if(m & 1ll) res = qmul(res, n, Mod);
		m >>= 1ll;
		n = qmul(n, n, Mod);
	}
	return res;
}

bool check(ll x) {
	if(x == 2) return 1;
	if(x % 2 == 0) return 0;
	if(x < 2) return 0;

	ll k = 0, cp = x - 1;
	while(cp % 2 == 0) k++, cp /= 2;

	bool flag = 1;
	for(int i = 1; i <= 12; i++) {
		ll d = qpow(prime[i], cp, x);
		if(d == 1) continue;
		if(x == prime[i]) return 1;
		if(x % prime[i] == 0) return 0;
		bool f = 0;
		for(int j = 1; j <= k; j++) {
			if(d == x - 1) {
				f = 1;
				break;
			}
			d = qmul(d, d, x);
		}
		if(!f) {
			flag = 0;
			break;
		}
	}
	return flag;
}

ll bigrand(ll n) {
	return 1ll * rand() * rand() % n + 1;
}
ll seed = 0;
ll go(ll x, ll n) {
	return ((i128)x * (i128)x + seed) % n;
}
ll floyd(ll n) {
	seed = bigrand(n);
	ll a, b; a = b = bigrand(n);
	b = go(b, n);
	while(a != b) {
		ll res = gcd(a - b + n, n);
		if(res != 1) return res;
		a = go(a, n);
		b = go(go(b, n), n);
	}
	return 1;
}
ll maxp = 0;
map <ll, int> pm;
void divide(ll n) {
	if(n == 1) return;
	if(check(n)) {
		pm[n]++;
		return ;
	}
	ll d = 0;
	while(1) {
		d = floyd(n);
		if(d != 1) break;
	}
	divide(d);
	divide(n / d);
}
struct numb {
	ll p; int c;
};

vector <numb> vec;
ll inv2 = 0;
ll f[(1 << NN) + 10], ans = 0;
bool hav2 = 0;
void dfs(int u, ll res, int state) {
	if(u >= vec.size()) {
		ll t = qpow(m, (res + 1) / 2) % Mod;
		if(res % 2 == 0)
			upd(ans, 1ll * t * f[state] % Mod * ((res / 2ll) % Mod) % Mod);
		else if(!hav2 || (hav2 && (state % 2 == 0)))
			upd(ans, 1ll * t * f[state] % Mod * (res % Mod) % Mod);
		return ;
	}
	ll bas = 1;
	for(int c = 0; c < vec[u].c; c++) {
		dfs(u + 1, res * bas, state | ((1 << (u - 1))));
		bas = 1ll * bas * vec[u].p;
	}
	dfs(u + 1, res * bas, state);
}
void init() {
	vec.clear(), pm.clear();

	cin >> n >> m >> Mod;
	inv2 = qpow(2, Mod - 2);
	divide(n);
	vec.push_back((numb){1, 0});
	for(map <ll, int>::iterator it = pm.begin(); it != pm.end(); it++)
		vec.push_back((numb){(*it).first, (*it).second});
	hav2 = (vec[1].p == 2);
	int lim = vec.size() - 1;
	for(int S = 0; S < (1 << lim); S++) {
		f[S] = 0;
		ll res = 1;
		for(int i = 1; i <= lim; i++)
			if((S >> (i - 1)) & 1) res = (i128)res * vec[i].p % Mod;
		int c = __builtin_popcount(S);
		if(c % 2 == 1) res = Mod - res;
		f[S] = res;
	}
	fwt(f, lim);
	ans = 0;
	dfs(1, 1, 0);
	cout << ans << '\n';
}
int main() {
	int T; cin >> T;
	while(T--) init();
}