【POJ 2154】Color

http://poj.org/problem?id=2154
还是先套上Burnside引理：$$\begin{aligned}
ans & =\sum_{i=1}^n n^{(i,n)-1} \
& = \sum_{d=1}^n [d|n]\sum_{i=1}^n [d|i]\left[\left(\frac id,\frac nd\right)=1\right]n^{d-1} \
& = \sum_{d=1}^n [d|n]n^{{d-1}\sum_{i=1}}\left[\left(i,\frac nd\right)=1\right] \
& = \sum_{d=1}^n [d|n]n^{d-1}\varphi\left(\frac nd\right)
\end{aligned}$$
因为n非常大，为了方便，求$\varphi$时进行质因子分解，求$\varphi(d)$单次$O\left(\frac{\sqrt d}{\log d}\right)$，时间复杂度$O\left(T\frac{n^{\frac 34}}{\log n}\right)$。
我是这么写的，时间复杂度非常不科学。还有一种更快的$O(1)$求$\varphi$的做法，就是先欧拉筛$\sqrt n$以内的$\varphi$并且预处理n的大于$\sqrt n$的质因子$p_{last}$（如果有的话），求$\varphi(d)$时如果$d>\sqrt n$，那么返回$\varphi\left(\frac d{p_{last}}\right)*\left(p_{last} - 1\right)$，否则直接返回$\varphi(d)$。时间复杂度$O\left(T\sqrt n\right)$。
这么简单的题我竟然说了这么多。。。

#include<cmath>
#include<cstdio>
#include<bitset>
#include<cstring>
#include<algorithm>
using namespace std;

const int N = 1000000003;

int n, p, sq, num = 0, prime[100003];
bitset <100003> notp;

void Euler_shai() {
	for (int i = 2; i <= sq; ++i) {
		if (!notp[i]) prime[++num] = i;
		for (int j = 1; j <= num && prime[j] * i <= sq; ++j) {
			notp[prime[j] * i] = 1;
			if (i % prime[j] == 0)
				break;
		}
	}
}

int ipow(int a, int b) {
	int ret = 1, w = a % p;
	while (b) {
		if (b & 1) (ret *= w) %= p;
		(w *= w) %= p;
		b >>= 1;
	}
	return ret;
}

int PHI(int nu) {
	int ret = 1;
	for (int i = 1; i <= num && prime[i] <= nu; ++i)
		if (nu % prime[i] == 0) {
			(ret *= ((prime[i] - 1) % p)) %= p;
			nu /= prime[i];
			while (nu % prime[i] == 0) {
				(ret *= (prime[i] % p)) %= p;
				nu /= prime[i];
			}
		}
	if (nu != 1) (ret *= ((nu - 1) % p)) %= p;
	return ret;
}

int work() {
	sq = ceil(sqrt(n)); int ret = 0;
	for (int i = 1; i < sq; ++i)
		if (n % i == 0) {
			(ret += ipow(n, i - 1) * PHI(n / i) % p) %= p;
			(ret += ipow(n, n / i - 1) * PHI(i) % p) %= p;
		}
	if (sq * sq == n) (ret += ipow(n, sq - 1) * PHI(sq) % p) %= p;
	return ret;
}

int main() {
	sq = sqrt(N);
	Euler_shai();
	
	int T; scanf("%d", &T);
	while (T--) {
		scanf("%d%d", &n, &p);
		printf("%d\n", work());
	}
	return 0;
}