题意

给你$n(n\leq 10^9)$n(n≤109)，$m(m\leq 2\times10^9)$m(m≤2×109)，$p(p\leq 2\times10^9 )$p(p≤2×109)，$p$p为质数求
$\prod_{i=1}^n\prod_{j=1}^n\prod_{k=1}^nm^{gcd(i,j)[k|gcd(i,j)]}$i=1∏n​j=1∏n​k=1∏n​mgcd(i,j)[k∣gcd(i,j)]
模$p$p的答案

思路

我们先只考虑指数部分和欧拉降幂那么有
$\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^ngcd(i,j)[k|gcd(i,j)]\ mod\ p-1$i=1∑n​j=1∑n​k=1∑n​gcd(i,j)[k∣gcd(i,j)] mod p−1
可以认为$\sum_{k=1}^ngcd(i,j)[k|gcd(i,j)]$k=1∑n​gcd(i,j)[k∣gcd(i,j)]
计算的是$gcd(i,j)*gcd(i,j)的约数个数$gcd(i,j)∗gcd(i,j)的约数个数，我们用$\sigma (d)$σ(d)表示$d$d的约数个数
那么原式就等于
$\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)\sigma(gcd(i,j))$i=1∑n​j=1∑n​gcd(i,j)σ(gcd(i,j))
令$d=gcd(i,j)$d=gcd(i,j)得
$\sum_{d=1}^nd\sigma(d)\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)==d]$d=1∑n​dσ(d)i=1∑n​j=1∑n​[gcd(i,j)==d]

$\sum_{d=1}^nd\sigma(d)\sum_{i=1}^{\left \lfloor\frac{n}{d}\right \rfloor}\sum_{j=1}^{\left \lfloor\frac{n}{d}\right \rfloor}[gcd(i,j)==1]$d=1∑n​dσ(d)i=1∑⌊dn​⌋​j=1∑⌊dn​⌋​[gcd(i,j)==1]
考虑后面的部分
$\sum_{i=1}^{\left \lfloor\frac{n}{d}\right \rfloor}\sum_{j=1}^{\left \lfloor\frac{n}{d}\right \rfloor}[gcd(i,j)==1]$i=1∑⌊dn​⌋​j=1∑⌊dn​⌋​[gcd(i,j)==1]
可以用莫比乌斯反演做但是这恰好$i$i和$j$j的上限都$n$n，可以认为是$1-n$1−n中互质的数对数，考虑每个$i$i的贡献其实就是$\varphi(i)$φ(i)欧拉函数的值，那么答案就是欧拉函数的前缀和，但是由于$i$i，$j$j的顺序是有影响的然后$(1,1)$(1,1)只计算一次所以有
$\sum_{i=1}^{\left \lfloor\frac{n}{d}\right \rfloor}\sum_{j=1}^{\left \lfloor\frac{n}{d}\right \rfloor}[gcd(i,j)==1]=2*\Phi(\left \lfloor\frac{n}{d}\right \rfloor)-1$i=1∑⌊dn​⌋​j=1∑⌊dn​⌋​[gcd(i,j)==1]=2∗Φ(⌊dn​⌋)−1
其中$\varphi$φ为欧拉函数的和函数
那么原式就转化为
$=\sum_{d=1}^nd\sigma(d)(2*\Phi(\left \lfloor\frac{n}{d}\right \rfloor)-1)$=d=1∑n​dσ(d)(2∗Φ(⌊dn​⌋)−1)
然后我们再考虑一下
$\sum_{i=1}^ni\sigma(i)$i=1∑n​iσ(i)
这个求和
可以把约数函数展开来
$=\sum_{i=1}^ni\sum_{d|i}1$=i=1∑n​id∣i∑​1
那么表示的就是对于每一个$i$i枚举他的因子有哪些，已知在$1-n$1−n中枚举$i$i的因子和枚举$i$i的倍数是等价的那么就有
$=\sum_{d=1}^n\sum_{d|i}i$=d=1∑n​d∣i∑​i
$i$i为$d$d的倍数，我们再把后面的展开来
$=\sum_{d=1}^n(d+2d+3d+\cdots+\left \lfloor\frac{n}{d}\right \rfloor d)$=d=1∑n​(d+2d+3d+⋯+⌊dn​⌋d)
后面就是一个等差数列
$=\sum_{d=1}^nd\frac{\left \lfloor\frac{n}{d}\right \rfloor(\left \lfloor\frac{n}{d}\right \rfloor+1)}{2}$=d=1∑n​d2⌊dn​⌋(⌊dn​⌋+1)​
令$Id(n)=\sum_{i=1}^ni\sigma(i)$Id(n)=i=1∑n​iσ(i)
原式想用分块优化的话那么有$last=\frac{n}{\frac{n}{i}}$last=in​n​有

$=\sum_{i=last+1}^n(Id(last)-Id(i-1))(2*\Phi(\left \lfloor\frac{n}{i}\right \rfloor)-1)$=i=last+1∑n​(Id(last)−Id(i−1))(2∗Φ(⌊in​⌋)−1)
$Id(n)$Id(n)可以用线性筛处理一部分(先求每个数的约数个数再求前缀和的时候乘以i)，大于的部分就用上面求的分块来处理了，欧拉函数求和用杜教筛即可
求完这个后最后的答案求一个快速幂就可以了，尽量都用int，不然会超时

#include<bits/stdc++.h>
#include <unordered_map>
using namespace std;
typedef long long ll;
const int N=10000005;
int mod;
unordered_map<int,int> P;
unordered_map<int,int> D;
bool isP[N];
int prime[N];
int cnt;
int phi[N];
int d[N];
int num[N];
long long inv=500000004;
void init()
{
    phi[1]=d[1]=1;
    for(int i=2; i<N; i++)
    {
        if(!isP[i])
        {
            prime[cnt++]=i;
            phi[i]=i-1;
            d[i]=2;
            num[i]=1;
        }
        for(int j=0; j<cnt&&(ll)i*prime[j]<N; j++)
        {
            isP[i*prime[j]]=true;
            if(i%prime[j])
            {
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
                d[i*prime[j]]=d[i]*d[prime[j]];
                num[i*prime[j]]=1;
            }
            else
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                num[i*prime[j]]=num[i]+1;
                d[i*prime[j]]=d[i]/(num[i]+1)*(num[i*prime[j]]+1);
                break;
            }
        }
    }
    for(int i=1; i<N; i++)
    {
        phi[i]=(phi[i]+phi[i-1])%mod;
        d[i]=((ll)i*d[i]+d[i-1])%mod;
    }
}
int Sum(int x)
{
    if(x<N)return phi[x];
    if(P[x])return P[x];
    int ans=(x+1ll)*x/2%mod;
    for(int i=2,last; i<=x; i=last+1)
    {
        last=x/(x/i);
        ans=(ans-(last-i+1ll)*Sum(x/i)%mod+mod)%mod;
    }
    ans=((ll)ans+mod)%mod;
    P[x]=ans;
    return ans;
}
int Id(int n)
{
    if(n<N) return d[n];
    else if(D[n]) return D[n];
    int ans=0;
    for(int i=1,last; i<=n; i=last+1)
    {
        last=n/(n/i);
        ans=(ans+(last+i)*(last-i+1ll)/2%mod*((n/i*(n/i+1ll)/2)%mod)%mod)%mod;
    }
    D[n]=ans;
    return ans;
}
long long quickmod(long long a,long long b,long long p)
{
    long long ans=1;
    while(b)
    {
        if(b%2==1)
            ans=ans*a%p;
        a=a*a%p;
        b=b/2;
    }
    return ans;
}
int main()
{
    int n,m,p;
    scanf("%d%d%d",&n,&m,&p);
    mod=p-1;
    init();
    int ans=0;
    for(int i=1,last; i<=n; i=last+1)
    {
        last=n/(n/i);
        ans=(ans+(Id(last)-Id(i-1)+mod)%mod*((2ll*Sum(n/i)%mod-1ll+mod)%mod)%mod+mod)%mod;
    }
    printf("%lld\n",quickmod(m,ans,p));
    return 0;
}
/*
1000000000 999999997 98765431
*/