The 2019 ACM-ICPC China Shannxi Provincial Programming Contest B. Product(杜教筛+约数)

题意

给你$n(n\le 10^9)$，$m(m\le 2\times 10^9)$，$p(p\le 2\times 10^9)$，$p$为质数求
$$\prod _{i=1}^n\prod _{j=1}^n\prod _{k=1}^n{m}^{gcd(i,j)[k|gcd(i,j)]}$$
模$p$的答案

思路

我们先只考虑指数部分和欧拉降幂那么有
$$\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^{n}gcd(i,j)[k|gcd(i,j)]modp-1$$
可以认为$$\sum _{k=1}^{n}gcd(i,j)[k|gcd(i,j)]$$
计算的是$gcd(i,j)\ast gcd(i,j)的约数个数$，我们用$\sigma (d)$表示$d$的约数个数
那么原式就等于
$$\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)\sigma (gcd(i,j))$$
令$d=gcd(i,j)$得
$$\sum _{d=1}^{n}d\sigma (d)\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==d]$$

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

$$=\sum _{i=last+1}^n(Id(last)-Id(i-1))(2\ast \Phi (\lfloor \frac{n}{i}\rfloor )-1)$$
$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
*/