【POJ 2154】Color

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

#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
", work());
	}
	return 0;
}