2019 西安邀请赛 B. Product（莫比乌斯反演 + 杜教筛 + 欧拉降幂）

如图：求 ${\prod }_{i=1}^n{\prod }_{j=1}^n{\prod }_{k=1}^n{m}^{gcd(i,j)[k|gcd(i,j)]}$

把求和符号放到幂次就变成求和符号，幂次式子变成：$$\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^{n}gcd(i,j)[k|gcd(i,j)]$$开始莫比乌斯反演：$$\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^{n}gcd(i,j)[k|gcd(i,j)]$$$$=\sum _{k=1}^{n}k\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{k}\rfloor }gcd(i,j)$$$$=\sum _{k=1}^{n}k\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{p|gcd(i,j)}\varphi (p)$$$$=\sum _{k=1}^{n}k\sum _{p=1}^{\lfloor \frac{n}{k}\rfloor }\varphi (p)\sum _{i=1}^{\lfloor \frac{n}{k\ast p}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{k\ast p}\rfloor }=\sum _{k=1}^{n}k\sum _{p=1}^{\lfloor \frac{n}{k}\rfloor }\varphi (p)\ast {\lfloor \frac{n}{k\ast p}\rfloor }^2$$$$=\sum _{T=1}^n\sum _{k|T}k\ast \varphi (\frac{T}{k})\ast {\lfloor \frac{n}{T}\rfloor }^2$$
令$S(T)={\sum }_{k|T}k\ast \varphi (\frac{T}{k})$，需要用杜教筛预处理S函数的前缀和。
$S=\varphi \ast Id$，$S\ast I=\varphi \ast Id\ast I=\mu \ast Id\ast I\ast Id=Id\ast Id=n\ast d(n)$
$n\ast d(n)$ 不能 O(1) 计算前缀和。
转化式子：${\sum }_{i=1}^{n}i\ast d(i)={\sum }_{i=1}^{n}i{\sum }_{j|i}1={\sum }_{j=1}^{n}j{\sum }_{i=1}^{\lfloor \frac{n}{j}\rfloor }i$
转化成这样可以暴力在$\sqrt{n}$时间内求和，加上记忆化，会使得杜教筛复杂度达到大概 $O(n^{\frac{3}{4}})$左右

分块求出幂次和之后求快速幂即可，幂次取模用欧拉定理，模数 mod = $\varphi (p)$

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 2e6 + 10;
bool ispri[maxn];
int pri[maxn],cnt[maxn],tot;
ll sum[maxn],tp[maxn];
ll n,m,mod,inv2,inv6,sqr;
ll ans[maxn];
map<ll,ll> mp,np;
ll fpow(ll a,ll b) {
	ll r = 1;
	while(b) {
		if(b & 1) r = r * a % mod;
		b >>= 1;
		a = a * a % mod;
	}
	return r;
}
void sieve(int n) {
	ispri[0] = ispri[1] = true;
	pri[0] = 0;sum[0] = 0;cnt[0] = 0;
	sum[1] = 1;
	for(int i = 2; i <= n; i++) {
		if(!ispri[i]) {
			pri[++pri[0]] = i;
			tp[i] = i;
			cnt[i] = 1;
			sum[i] = 2 * i - 1;
		}
		for(int j = 1; j <= pri[0] && i * pri[j] <= n; j++) {
			ispri[i * pri[j]] = true;
			if(i % pri[j] == 0) {
				tp[i * pri[j]] = tp[i] * pri[j];
				cnt[i * pri[j]] = cnt[i] + 1;
				ll c = cnt[i * pri[j]];
				ll tmp = ((c + 1) * tp[i * pri[j]] % mod - c * tp[i] % mod + mod) % mod;
				sum[i * pri[j]] = sum[i / tp[i]] * tmp % mod;
				break;
			}
			tp[i * pri[j]] = pri[j];
			cnt[i * pri[j]] = 1;
			sum[i * pri[j]] = sum[i] * sum[pri[j]] % mod;
		}
	}
	for(int i = 1; i <= n; i++)
		sum[i] = (sum[i] + sum[i - 1]) % mod;
}
ll cal(ll x) {
	return x * (x + 1) / 2 % mod;
}
ll cal2(ll x) {
	return x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;
}
ll getans(ll x) {
	ll tmp = 0;
	for(ll i = 1,j; i <= x; i = j + 1) {
		j = x / (x / i);
		tmp += (cal(j) - cal(i - 1) + mod) % mod * cal(x / i) % mod;
		tmp %= mod;
	}
	return tmp;
}
ll getsum(ll x) {
	if(x <= 2000000) return sum[x];
	if(mp[x]) return mp[x];
	ll res = getans(x);
	for(ll i = 2,j; i <= x; i = j + 1) {
		j = x / (x / i);
		res = (res + mod - (j - i + 1) * getsum(x / i) % mod) % mod;
	}
	return mp[x] = res;
}
int main() {
	scanf("%lld%lld%lld",&n,&m,&mod);
	mod--;
	sieve(2000000);
	inv2 = fpow(2,mod - 2);inv6 = fpow(6,mod - 2);
	ll tot = 0;
	for(ll i = 1,j; i <= n; i = j + 1) {
		j = n / (n / i);
		tot += (getsum(j) + mod - getsum(i - 1)) % mod * (n / i) % mod * (n / i) % mod;
		tot %= mod; 
	}
	mod++;
	printf("%lld\n",fpow(m,tot));
	return 0;
}