2019年西安邀请赛B：Product 【杜教筛+降幂】

题目:

2019年西安邀请赛B：Product

题意：

给定 N，M，P，求下面式子的值最后模P：

分析：

m 的值恒定，先考虑求m最后的幂，即：

上式中 d(x) 表示 x 的约数的个数；F 函数的值是一段一段的，每一段都相等，如果知道了 g*d(g) 的前缀和就可以分块加速枚举 g； 设其前缀和 G(n) 为：

完美！！！这个求和式是可以分块快速计算的；就剩下 F 函数的值，如果用反演求，但 n 比较大，还是要用杜教筛计算莫比乌斯函数的前缀和，实际上 F(n) 也等价于：

欧拉函数前缀和是可以杜教筛快速求的；令 h(n) = (I * phi)(n)，（I 为恒等函数，即：I(x) = 1）写成狄利克雷卷积展开式有 h(n)：

 

最后一步是欧拉函数的性质，那么就有：

提出左边式子的第一项，得到：

设 S(n) 为欧拉函数的前缀和，有：

上面已经推出 h(x) = x，那么递归求S(n)即可；

降幂直接对（p-1）取模就好了；为什么呢？根据广义欧拉降幂公式可知如果 gcd(m,p) = 1时，直接取余是成立的；如果不等于1呢？只可能为：m 是 p 的倍数，此时不管幂是多少答案都是0，当然幂也不能为0，依然成立

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;
const int maxn = 1e6+56;
const int maxm = 7e5+75;
LL n,m,p,phi[maxn],prime[maxm],d[maxn];
unordered_map<LL,LL> mp;
LL qpow(LL a,LL x,LL mod){
    LL res = 1ll;
    while(x){
        if(x&1) res = res*a%mod;
        a = a*a%mod;
        x >>= 1;
    }
    return res;
}
void init(){                              //预处理G()和S()函数的小范围前缀和
    phi[1] = 1; int cnt = 0;
    for(int i = 2;i < maxn; ++i){
        if(!phi[i]) phi[i] = i-1,prime[cnt++] = i;
        for(int j = 0;j<cnt&&prime[j]*i<maxn; ++j){
            if(i%prime[j] == 0){
                phi[i*prime[j]] = phi[i]*prime[j];
                break;
            }
            else phi[i*prime[j]] = phi[i] * phi[prime[j]];
        }
    }
    for(int i = 1;i < maxn; ++i)
    for(int j = i;j < maxn; j += i)
        d[j]++;
    for(int i = 2;i < maxn; ++i){
        phi[i] += phi[i-1];
        if(phi[i] >= p) phi[i] -= p;
        d[i] = d[i-1]+d[i]*i%p;
        if(d[i] >= p) d[i] -= p;
    }
}
LL Sum(int n){                        //计算S(n)
    if(n < maxn) return phi[n];
    if(mp.count(n)) return mp[n];
    LL res = 0;
    for(int i = 2,up;i <= n; i=up+1){
        int t = n/i; up = n/t;
        res += 1ll*(up-i+1)*Sum(t)%p;
        if(res >= p) res -= p;
    }
    return mp[n] = (1ll*n*(n+1)/2-res+p)%p;
}
LL cal(int n){                        //计算G(n)
    if(n < maxn) return d[n];
    LL res = 0;
    for(int i = 1,up;i <= n; i=up+1){
        int t = n/i; up = n/t;
        res += 1ll*t*(t+1)/2%p*(1ll*(up+i)*(up-i+1)/2%p)%p;
        if(res >= p) res -= p;
    }
    return res;
}
LL solve(int n){
    LL res = 0,per = 0;
    for(int i = 1,up;i <= n; i=up+1){
        int t = n/i; up = n/t;
        LL tep = cal(up);
        res += (tep-per+p)*(Sum(t)*2-1+p)%p;
        if(res >= p) res -= p; 
        per = tep;
    }
    return res;
}
int main(){
    cin >> n >> m >> p; p--; 
    if(m % p == 0) return puts("0"),0;
    else init();
    cout << qpow(m,solve(n),p+1) << endl;
    return 0;
}