ICPC2019 西安邀请赛 B.Product
题目链接
https://nanti.jisuanke.com/t/39269
题目大意
计算 $\prod\limits_{i=1}^n\prod\limits_{j=1}^n\prod\limits_{k=1}^n{m ^ {gcd(i, j)[k | gcd(i, j)]}}$ ，答案对质数 $p$ 取模。
题目分析
由于 $p$ 是质数，所以由费马小定理可知 $m ^ {p – 1} = 1$ ，因此只需要计算出质数部分对 $p – 1$ 取模的答案即可。因此只需要计算下式对 $p – 1$ 取模的答案即可：
$$\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{k=1}^n{gcd(i, j)[k | gcd(i, j)]}$$
$k$ 的部分很明显是求 $gcd(i, j)$ 的约数个数，所以可以化简为
$$\sum\limits_{i=1}^n\sum\limits_{j=1}^n{gcd(i, j)d(gcd(i, j))}$$
那么可以考虑枚举 $gcd(i, j)$ 的值，
$$\sum\limits_{x=1}^n{xd(x)\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{x}\rfloor}{[gcd(i, j) = 1]}}$$
显然后半部分可以化简为 $\varphi$ 函数，即
$$\sum\limits_{x=1}^n{xd(x)((\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}{2\varphi(i)}) – 1)}$$
可以发现可以使用分块求和计算。后半部分可以使用杜教筛来解决。再来考虑 $\sum\limits_{x=1}^n{xd(x)}$，考虑枚举因数，则显然有
$$\sum\limits_{d=1}^n\sum\limits_{d | x}x$$
可以化为
$$\sum\limits_{d=1}^n{d\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}x}$$
也可以使用分块求和解决。
注意这一部分计算的复杂度还是不小的，需要在预处理时就预处理出来一部分，筛的时候记录下来已经筛过的。
最后在用快速幂计算答案即可。
代码
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <string>
#include <algorithm>
#include <map>
#include <unordered_map>
using namespace std;
const int maxn = 5e6 + 5;
const int m = 5e6;
int n, k, p;
int prime[maxn / 10];
bool isnp[maxn];
int phi[maxn];
int d[maxn];
int pc = 0;
int mod;
unordered_map <int, int> phim;
unordered_map <int, int> dm;
int init(){
int i, j;
phi[0] = 0;
phi[1] = 1;
d[1] = 1;
for(i=2;i<=m;i++){
if(!isnp[i]){
pc++;
prime[pc] = i;
phi[i] = i - 1;
d[i] = 2;
}
for(j=1;j<=pc and i * prime[j]<=m;j++){
isnp[i * prime[j]] = true;
if(i % prime[j] == 0){
phi[i * prime[j]] = phi[i] * prime[j];
int cnt = 1, y = i;
while(y % prime[j] == 0){
y /= prime[j];
cnt++;
}
d[i * prime[j]] = d[i] / cnt * (cnt + 1);
break;
}else{
phi[i * prime[j]] = phi[i] * phi[prime[j]];
d[i * prime[j]] = d[i] * 2;
}
}
}
for(i=2;i<=m;i++){
phi[i] = (phi[i] + phi[i - 1]) % mod;
d[i] = 1ll * d[i] * i % mod;
d[i] = (d[i] + d[i - 1]) % mod;
}
return 0;
}
long long qpow(long long b, long long x){
long long ret = 1;
while(x){
if(x & 1){
ret = ret * b % mod;
}
b = b * b % mod;
x >>= 1;
}
return ret;
}
int calc(int x){
int ret = (1ll * x * (x + 1) / 2) % mod;
int i, j;
if(x < m){
return phi[x];
}
if(phim.count(x)){
return phim[x];
}
for(i=2;i<=x;i=j+1){
j = x / (x / i);
ret = ((ret - 1ll * (j - i + 1) * calc(x / i)) % mod + mod) % mod;
}
return phim[x] = ret;
}
int sumd(int t){
int i, j;
int ret = 0;
if(t < m){
return d[t];
}
if(dm.count(t)){
return dm[t];
}
for(i=1;i<=t;i+=j){
j = t / (t / i) - i + 1;
ret += 1ll * (i + i + j - 1) * j / 2 % mod * ((1ll * (1 + t / i) * (t / i) / 2) % mod) % mod;
ret %= mod;
}
return dm[t] = ret;
}
int main(){
int i, j;
int ans = 0;
cin >> n >> k >> p;
mod = p - 1;
init();
for(i=1;i<=n;i+=j){
j = n / (n / i) - i + 1;
ans += 1ll * (sumd(i + j - 1) - sumd(i - 1)) % mod * ((2 * calc(n / i) - 1) % mod) % mod;
ans = (ans % mod + mod) % mod;
}
mod = p;
cout << qpow(k, ans) << endl;
return 0;
}