2019ICPC西安邀请赛 - B. Product - 数论
打印的时候麻烦把：https://blog.csdn.net/skywalkert/article/details/50500009这个打印下来。
求\(\prod\limits_{i=1}^{n} \prod\limits_{j=1}^{n} \prod\limits_{k=1}^{n} m^{gcd(i,j)[k|gcd(i,j)]} mod p\)
欧拉定理：
\(m^{ \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} \sum\limits_{k=1}^{n} {gcd(i,j)[k|gcd(i,j)]} mod (p-1) }mod p\)
算出指数直接套快速幂就可以了。考虑指数怎么算。为了方便把取模省略。
$\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} \sum\limits_{k=1}^{n} {gcd(i,j)[k|gcd(i,j)]} $
这种带个整除的，一般都是交换求和到前面，然后后面算他的倍数。
$=\sum\limits_{k=1}^{n} \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} {gcd(i,j)[k|gcd(i,j)]} $
带方括号的，枚举g，去括号。把g放在前面。
$=\sum\limits_{g=1}{n}g\sum\limits_{k=1} \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} {[gcd(i,j)==g][k|g]} $
k与i,j无关了，提到前面。
$=\sum\limits_{g=1}{n}g\sum\limits_{k=1} [k|g] \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} {[gcd(i,j)==g]} $
那么这个k就是d(g)了。
$=\sum\limits_{g=1}^{n}g*d(g) \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} {[gcd(i,j)==g]} $
后面那个就把g提出来，好像还是很套路的。在什么鬼GCD统计里见过下式：
$\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} {[gcd(i,j)g]} \(
\) = \sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor} \sum\limits_{j=1}^{\lfloor\frac{n}{g}\rfloor} {[gcd(i,j)1]} $
那么这个两个数的互质对，规定i比j大于等于，那么就是欧拉函数，特别的，\(\varphi(1)=1\)。(i,j)和(j,i)互换重复计算(1,1)。
$ = \sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor} \sum\limits_{j=1}^{\lfloor\frac{n}{g}\rfloor} {[gcd(i,j)==1]} \(
\) = 2 ( \sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor} \varphi(i) ) - 1 $
这个东西用杜教筛板子都不用思考任何东西，直接弄出来（说起来我的杜教筛好像多个log）。
记上面的算式为 \(S(\lfloor\frac{n}{g}\rfloor)\) ，那么原来的简化为：
$=\sum\limits_{g=1}^{n}g*d(g) S(\lfloor\frac{n}{g}\rfloor) $
换个喜欢的字母！！！
$=\sum\limits_{i=1}^{n}i*d(i) S(\lfloor\frac{n}{i}\rfloor) $
这些都跟g有关，太自闭了，大半时间卡在这里。后面的函数明显可以分块到 \(O(\sqrt{n})\)，对于一段连续的i求出结果。
要是可以计算前面的 \(\sum\limits_{i=1}^{n}i*d(i)\) 的前缀和，那么可以再花\(O(1)\)求出一段连续的i的求和结果和后面相乘。
最后的问题就是怎么快速计算 \(\sum\limits_{i=1}^{n}i*d(i)\) 这个狗东西。
首先可以线性预处理到5e6左右，花费一大堆内存，因为d(i)我是有线性筛的（虽然不会他的杜教筛）。
（菜鸡的）第一次重要的突破！
每个数的贡献和他的因子个数成线性。
然后我突发奇想，枚举每个因子d，枚举到 \(\sqrt{n}\) ，每个因子会使得他的每个倍数都恰好多贡献1倍，那么就是 \(d+2d+3d+...+\lfloor\frac{n}{d}\rfloor * d\)。
等差数列求和嘛！
那么就变成：
\(\sum\limits_{i=1}^{n}i*d(i)\)
\(=\sum\limits_{d=1}^{\lfloor\sqrt{n}\rfloor} \frac{d*(1+\lfloor\frac{n}{d}\rfloor)*(\lfloor\frac{n}{d}\rfloor)}{2}\)
已经变成 \(O(\sqrt{n})\)求的了，当时出不了第三个样例，怎么想也是哦，能出就奇怪了。
（菜鸡的）第二次重要的突破！
233，看了半天才发现，上面不是也可以分块求的吗？对\(\lfloor\frac{n}{d}\rfloor\)分块啊。对一段连续的d，\(\lfloor\frac{n}{d}\rfloor\)都是同一个值，想了半天才领会！
那么就可以 \(O(n^{\frac{1}{4}})\) 求出\(\sum\limits_{i=1}^{n}i*d(i)\) 了，把这个狗东西记为 \(T(i)\)
$=\sum\limits_{i=1}^{n}T(i) S(\lfloor\frac{n}{i}\rfloor) $
杜教筛预处理 \(O(n^{\frac{2}{3}})\) ，单次询问S(x)差不多 \(O(1)\) ，虽然我的杜教筛比dq的多了一点东西。
单次询问T(i)差不多 \(O(n^{\frac{1}{4}})\) 。意思是对单个i，处理出答案需要 \(O(n^{\frac{1}{4}})\) 。
求和时，第一次分块对i分块用了 \(O(\sqrt{n})\)，每段连续的i共用一个T(i) S(\lfloor\frac{n}{i}\rfloor) ，所以就是 \(O(n^{\frac{3}{4}})\) 。n是1e9，刚刚好够。
其他注意点：
1.应该用欧拉定理，这里卡了很久！幂取模不是直接取的，不要搞假算法！
2.涉及除法的运算，恰好求和公式可以约掉那个2，不要对那个数字取模！
总结：
\(\sum\limits_{i=1}^{n}i*d(i)=\sum\limits_{d=1}^{\lfloor\sqrt{n}\rfloor} \frac{d*(1+\lfloor\frac{n}{d}\rfloor)*(\lfloor\frac{n}{d}\rfloor)}{2}\) 可以 \(O(n^{\frac{1}{4}})\) 计算。
同理可得 \(\sum\limits_{i=1}^{n}d(i)=\sum\limits_{d=1}^{\lfloor\sqrt{n}\rfloor} \frac{(1+\lfloor\frac{n}{d}\rfloor)*(\lfloor\frac{n}{d}\rfloor)}{2}\) 也是一样的。
说白了是我不会计算\(d(i)\)前缀和的方法，过于依赖线性筛了，去关注一下单个xx函数的求法吧。
当时写的代码，270分钟1A，好像挤到金银分隔线。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN=5e6;
ll n,m,p,mod;
ll inv2;
ll qpow(ll x,ll n)
{
ll res=1%p;
while(n)
{
if(n&1)
res=(res*x)%p;
x=(x*x)%p;
n>>=1;
}
return res%p;
}
ll Qarray[MAXN+5];
//Q(n)=\sum_i=1^n i*d(i)
inline ll Q(ll n)
{
if(n<=MAXN)
return Qarray[n];
ll res=0;
for(ll l=1,r;l<=n; l=r+1){
r=n/(n/l);
ll t=n/l;
ll tmp1=((t+1)*t/2)%mod;
ll tmp2=((l+r)*(r-l+1)/2)%mod;
res=(res+(tmp1*tmp2%mod))%mod;
}
return res%mod;
}
unordered_map<int,ll> Sphi;
ll sump[MAXN+5];
int pri[MAXN+5],pritop;
bool notpri[MAXN+5];
ll d[MAXN+5],a[MAXN+5];
void sieve(int n=MAXN)
{
pritop=0;
notpri[1]=1;
sump[1]=1;
d[1]=1;
for(int i=2; i<=n; i++)
{
if(!notpri[i])
{
pri[++pritop]=i;
d[i]=2;
a[i]=1;
sump[i]=(i-1)%mod;
}
for(int j=1; j<=pritop&&i*pri[j]<=n; j++)
{
notpri[i*pri[j]]=1;
if(i%pri[j])
{
d[i*pri[j]]=d[i]*d[pri[j]];
a[i*pri[j]]=1;
sump[i*pri[j]]=sump[i]*sump[pri[j]]%mod;
}
else
{
d[i*pri[j]]=d[i]/(a[i]+1)*(a[i]+2);
a[i*pri[j]]=a[i]+1;
sump[i*pri[j]]=sump[i]*pri[j]%mod;
break;
}
}
}
Qarray[0]=0;
Qarray[1]=1;
for(int i=2; i<=n; i++)
{
sump[i]=(sump[i-1]+sump[i])%mod;
Qarray[i]=(Qarray[i-1]+1ll*i*d[i])%mod;
}
return;
}
inline ll GetSphi(int n)
{
if(n<=MAXN)
return sump[n];
if(Sphi.count(n))
return Sphi[n];
ll ret=n;
if(ret%2==0)
ret=(ret/2)*(1ll+n);
else
ret*=(1ll+n)/2;
for(ll l=2,r;l<=n; l=r+1)
{
ll t=n/l;
r=n/(t);
ret-=(r-l+1)*GetSphi(t);
}
ret%=mod;
return Sphi[n]=ret;
}
int main()
{
#ifdef local
freopen("Yinku.in","r",stdin);
#endif // local
scanf("%lld%lld%lld",&n,&m,&p);
mod=p-1;
sieve();
ll ans=0;
for(ll l=1,r; l<=n; l=r+1)
{
ll t=n/l;
r=n/t;
ll tmp=(Q(r)-Q(l-1)+mod)%mod;
ll tmp2=GetSphi(t)%mod;
ans=(ans+tmp*tmp2%mod)%mod;
}
//printf("ans=%lld\n",ans);
ans=(ans*2ll)%mod;
ans=(ans-Q(n)+mod)%mod;
//printf("ans=%lld\n",ans);
ll Ans=qpow(m,ans%mod);
printf("%lld\n",Ans);
return 0;
}
