打印的时候麻烦把：。

求\(\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}^{n} \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}^{n} [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
    
     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
     
       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;}