题解 P5518 【[MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题】

前言

请确保你在做此题之前做过

P1390 公约数的和

P2303 [SDOI2012]Longge的问题

P1829 [国家集训队]Crash的数字表格 / JZPTAB

P3704 [SDOI2017]数字表格

并有一定的反演芝士已获得更好的做题体验

由于是断断续续做的可能柿子写的有不对的地方看出来的人可以指出来，谢谢~

简化题面

求

[mathcal{prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}left(frac{ ext{lcm}(i,j)}{gcd(i,k)} ight)^{f(type)}} ]

对于(mathcal{f(x)})有

[f(0)=1 ]

[f(1)=i imes j imes k ]

[f(2)=gcd(i,j,k) ]

对于(type)仅有(0,1,2)三种取值

对于数据

[mathcal{1leq A,B,Cleq 10^5 10^7 leq p leq 1.05 imes 10^9 pin { prime} T =70} ]

解题思路

第一眼看就是大毒瘤题

首先对原柿子进行展开

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} left(frac{i imes j}{gcd(i,j) imes gcd(i,k)} ight)^{f(type)} ]

很显然的可以拆成两个柿子

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} i^{f(type)} ]

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} gcd(i,j)^{f(type)} ]

所以变成了对于六个子问题的化简，慢慢梳理

Type=0

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} i ]

[=prod_{i=1}^{A}prod_{j=1}^{B}i^C ]

[=(prod_{i=1}^{A}i)^{C imes B} ]

直接预处理阶乘就可以求得

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} gcd(i,j) ]

[=(prod_{i=1}^{A}prod_{j=1}^{B} gcd(i,j))^C ]

考虑像数字表格一样进行处理

[=(prod_{d=1} d^{sum_{i=1}^{A/d}sum_{j=1}^{B/d} [gcd(i,j)=1]})^C ]

[=(prod_{d=1} d^{sum_{p=1}mu(p)leftlfloorfrac{A}{dp} ight floorleftlfloorfrac{B}{dp} ight floor})^C ]

[=(prod_{Q=1}(prod_{d|Q}d^{mu(frac{Q}{d})})^{leftlfloorfrac{A}{Q} ight floorleftlfloorfrac{B}{Q} ight floor})^{C} ]

预处理然后分块做即可

Type=1

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} i^{i imes j imes k} ]

[=prod_{i=1}^{A}prod_{j=1}^{B} i^{i imes j imes sum_{k=1}^{C}k} ]

[=prod_{i=1}^{A} i^{i imes sum_{j=1}^{B}j imes sum_{k=1}^{C}k} ]

质数部分处理一个等差数列求，并且预处理(i*i)的阶乘直接求就可以

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} gcd(i,j)^{i imes j imes k} ]

第一眼看真不知道从哪里下手，思考之后考虑先把没用的东西拿出来

[=prod_{i=1}^{A}prod_{j=1}^{B} gcd(i,j)^{i imes j imes sum_{k=1}^{C}k} ]

然后你会感觉这个柿子有点眼熟，很像Crash的数字表格，于是我们按照当时的想法进行计算

[egin{aligned}=(prod_{d=1}d^{d^2sum_{i=1}^{A/d}sum_{j=1}^{B/d} i imes j imes [gcd(i, j)=1]})^{sum_{k=1}^{C}k}end{aligned} ]

[egin{aligned}=(prod_{d=1}d^{sum_{p=1}d^2mu(p)p^2sum(leftlfloorfrac{A}{dp} ight floor)sum(leftlfloorfrac{B}{dp} ight floor)})^{sum_{k=1}^{C}k}end{aligned} ]

[egin{aligned}=(prod_{Q=1}((prod_{d|Q} d^{mu(frac{Q}{d})})^{Q^2})^{sum(leftlfloorfrac{A}{Q} ight floor)sum(leftlfloorfrac{B}{Q} ight floor)})^{sum_{k=1}^{C}k}end{aligned} ]

显然的和(type0)一样进行预处理然后分块做即可

Type=2

这道题最大的亮点最大的毒瘤的地方

老套路枚举(gcd)

[egin{aligned}prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} i^{gcd(i,j,k)}end{aligned} ]

[egin{aligned}=prod_{d=1}prod_{i=1}^A i^{d imes sum_{j=1}^Bsum_{k=1}^C[gcd(i,j,k)=d]}end{aligned} ]

[egin{aligned}=prod_{d=1}prod_{i=1}^{(leftlfloorfrac{A}{d} ight floor)} i imes d ^{d imes sum_{j=1}^{(leftlfloorfrac{B}{d} ight floor)}sum_{k=1}^{(leftlfloorfrac{C}{d} ight floor)}[gcd(i,j,k)=1]}end{aligned} ]

[=prod_{d=1}prod_{i=1}^{(leftlfloorfrac{A}{d} ight floor)} (i imes d) ^{d imes sum_{p=1}mu(p)(leftlfloorfrac{B}{dp} ight floor)(leftlfloorfrac{C}{dp} ight floor)} ]

考虑把指数上的(sum)拿到下面来

[=prod_{d=1}prod_{p=1}prod_{i=1}^{(leftlfloorfrac{A}{dp} ight floor)} (i imes d imes p) ^{d imesmu(p)(leftlfloorfrac{B}{dp} ight floor)(leftlfloorfrac{C}{dp} ight floor)} ]

[=prod_{Q=1}(prod_{d|Q}(prod_{i=1}^{(leftlfloorfrac{A}{Q} ight floor)}i imes Q)^{d imes mu(leftlfloorfrac{Q}{d} ight floor)})^{(leftlfloorfrac{B}{Q} ight floor)(leftlfloorfrac{C}{Q} ight floor)} ]

考虑再把(prod)拿上去

[=prod_{Q=1}((prod_{i=1}^{(leftlfloorfrac{A}{Q} ight floor)}i imes Q)^{d imes sum_{d|Q}mu(leftlfloorfrac{Q}{d} ight floor)})^{(leftlfloorfrac{B}{Q} ight floor)(leftlfloorfrac{C}{Q} ight floor)} ]

然后你会发现(d imes sum_{d|Q}mu(leftlfloorfrac{Q}{d} ight floor))非常的眼熟，(id*mu=varphi)

[=prod_{Q=1}((prod_{i=1}^{(leftlfloorfrac{A}{Q} ight floor)}i imes Q)^{varphi(Q)})^{(leftlfloorfrac{B}{Q} ight floor)(leftlfloorfrac{C}{Q} ight floor)} ]

预处理(varphi)的前缀和，用之前求得阶乘乱搞就可以了

[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C} gcd(i,j)^{gcd(i,j,k)} ]

这个题最恶心人的柿子，还是考虑枚举(gcd)，先考虑枚举(gcd(i,j))

[=prod_{d=1}^{A} d^{sum_{i=1}^Asum_{j=1}^B [gcd(i,j)=d] imes (sum_{k=1}^{C} gcd(d,k))} ]

[=prod_{d=1}d^{sum_{i=1}^{(leftlfloorfrac{A}{d} ight floor)}sum_{j=1}^{(leftlfloorfrac{B}{d} ight floor)}[gcd(i,j)=1] imes (sum_{k=1}^{C} gcd(d,k))} ]

[=prod_{d=1}d^{(sum_{p=1}mu(p)leftlfloorfrac{A}{dp} ight floorleftlfloorfrac{B}{dp} ight floor) imes (sum_{k=1}^{C} gcd(d,k))} ]

[=prod_{Q=1}(prod_{d|Q}d^{mu(leftlfloorfrac{Q}{d} ight floor) imes (sum_{k=1}^{C} gcd(d,k))})^{leftlfloorfrac{A}{Q} ight floorleftlfloorfrac{B}{Q} ight floor} ]

对(sum_{k=1}^{C} gcd(d,k))继续进行反演

你会发现这不就是(Longge)的问题里面的那个柿子吗

[sum_{d'|d}d'sum_{k=1}^{C/d'} [gcd(frac{d}{d'},k) = 1] ]

[sum_{d'|d}d'sum_{p'|d,d'}mu(p')leftlfloorfrac{C}{d'p'} ight floor ]

化简后可得到

[sum_{T|d}leftlfloorfrac{C}{T} ight floorvarphi(T) ]

然后套回去之后变成了

[=prod_{Q=1}(prod_{d|Q}d^{mu(leftlfloorfrac{Q}{d} ight floor) imes (sum_{T|d}leftlfloorfrac{C}{T} ight floorvarphi(T))})^{leftlfloorfrac{A}{Q} ight floorleftlfloorfrac{B}{Q} ight floor} ]

化简到这里复杂度是(O(nlogn))的，我试了试貌似过不了。。。然后我就呆住了，这怎么继续反演

看了看你谷的题解之后我悟了

把(d)拆成(frac{d}{T})和(T)来做

把指数上的(sum)拿下来

[=prod_{Q=1}(prod_{d|Q}prod_{T|d}d^{mu(leftlfloorfrac{Q}{d} ight floor) imes (leftlfloorfrac{C}{T} ight floorvarphi(T))})^{leftlfloorfrac{A}{Q} ight floorleftlfloorfrac{B}{Q} ight floor} ]

首先对于(frac{d}{T})我们有

[prod_{Q=1}(prod_{d|Q}prod_{T|d}(frac{d}{T})^{mu(leftlfloorfrac{Q}{d} ight floor) imes (leftlfloorfrac{C}{T} ight floorvarphi(T))})^{leftlfloorfrac{A}{Q} ight floorleftlfloorfrac{B}{Q} ight floor} ]

对于(Q)和(d)我们同时乘上(T)可以得到

[=prod_{T=1}^x prod_{Q=1}^{leftlfloor{frac{x}{T}} ight floor} prod_{d|Q} d^{mu(frac{QT}{dT})varphi(T)leftlfloor{frac{z}{T}} ight floor{leftlfloor{frac{x}{QT}} ight floor}{leftlfloor{frac{y}{QT}} ight floor}}]

[=prod_{T=1}^x prod_{Q=1}^{leftlfloor{frac{x}{T}} ight floor} prod_{d|Q} d^{mu(frac{Q}{d})varphi(T)leftlfloor{frac{z}{T}} ight floor{leftlfloor{frac{x}{QT}} ight floor}{leftlfloor{frac{y}{QT}} ight floor}}]

然后里面的部分我们在(type0)中已经预处理过了，直接白嫖

对于(T)的话

[prodlimits_{Q=1}^xprod_{d|Q}prod_{T|d}T^{mu(frac{Q}{d})varphi(p)leftlfloor{frac{z}{T}} ight floor{leftlfloor{frac{x}{Q}} ight floor}{leftlfloor{frac{y}{Q}} ight floor}} ]

枚举(T)

[=prod_{T=1}^xprod_{Q=1}^{leftlfloor{frac{x}{T}} ight floor}prod_{d|Q}p^{mu(frac{QT}{dT})varphi(T)leftlfloor{frac{z}{T}} ight floor{leftlfloor{frac{x}{QT}} ight floor}{leftlfloor{frac{y}{QT}} ight floor}} ]

然后把(prod)移动到指数上

[=prod_{T=1}^xT^{sumlimits_{Q=1}^{leftlfloor{frac{x}{T}} ight floor}sumlimits_{d|Q}mu(frac{Q}{d})varphi(T)leftlfloor{frac{z}{T}} ight floor{leftlfloor{frac{x}{QT}} ight floor}{leftlfloor{frac{y}{QT}} ight floor}} ]

然后整理一下可以得到

[=prod_{T=1}^xT^{varphi(T)leftlfloor{frac{z}{T}} ight floorsumlimits_{Q=1}^{leftlfloor{frac{x}{Q}} ight floor}{leftlfloor{frac{x}{QT}} ight floor}{leftlfloor{frac{y}{QT}} ight floor}sumlimits_{d|Q} mu(frac{Q}{d})}]

然后最妙的地方是(sumlimits_{d|Q}mu(frac{Q}{d}))可以替换成(e)，也就是([Q=1])，所以当且仅当(Q=1)的时候才会有贡献，所以就变成了

[=prod_{T=1}^xT^{varphi(T){leftlfloor{frac{x}{T}} ight floor}{leftlfloor{frac{y}{T}} ight floor}{leftlfloor{frac{z}{T}} ight floor}} ]

然后这一部分就可以算了，直接分块

细节处理

预处理阶乘逆元记得倒叙，其他的就没了

(Code)

#include<bits/stdc++.h>
    
#define LL long long
    
#define _ 0

#define R register

#define int long long

using namespace std;
    
/*Grievous Lady*/
    
template <typename T> void read(T & t){
    t = 0;int f = 1;char ch = getchar();
    while(ch < '0' || ch > '9'){if(ch == '-')f =- 1;ch = getchar();}
    do{t = (t << 1) + (t << 3) + ch - '0';ch = getchar();}while(ch >= '0' && ch <= '9');t *= f;
}
    
const int kato = 1e5 + 10;

const int atri = 1e5 + 10;

int a , b , c , t , mod;

int prime[kato] , cnt , phi[kato] , mu[kato] , fac[kato] , Fac[kato] , type0[kato] , type1[kato];

int phi_sum[kato] , fac_0[kato] , fac_1[kato] , fac_0_inv[kato] , fac_1_inv[kato] , inv[kato];

bool ispri[kato];

inline LL quick_pow(LL a , LL b){
    LL res = 1;
    for(; b ; b >>= 1 , a = a * a % mod){
        if(b & 1){
            res = res * a % mod;
        }
    }
    return res % mod;
}

// void put(int x){
// 	if(x < 0) putchar('-') , x =- x;
// 	if(x >= 10) put(x / 10);
// 	putchar((x % 10) + 48);
// }

inline void pri(){
    for(int i = 1;i < atri;i ++){
        fac[i] = fac[i - 1] * i % mod;
        Fac[i] = Fac[i - 1] * quick_pow(i , i) % mod;
    }
    for(int i = 2;i < atri;i ++){
        if(!ispri[i]){
            prime[++ cnt] = i;
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for(int j = 1;j <= cnt && i * prime[j] < atri;j ++){
            ispri[i * prime[j]] = 1;
            if(i % prime[j] == 0){
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i = 1;i < atri;i ++){
        phi_sum[i] = (phi_sum[i - 1] + phi[i]) % (mod - 1);
        type0[i] = fac_0[i] = fac_1[i] = fac_0_inv[i] = fac_1_inv[i] = type1[i] = 1;
    }
    for(int i = 1;i < atri;i ++){
        if(mu[i] == 1){
            for(int j = i , cnt = 1;j < atri;j += i , cnt ++){
                type0[j] = type0[j] * cnt % mod;
                type1[j] = type1[j] * cnt % mod;
            }
        }
        if(mu[i] == -1){
            for(int j = i , cnt = 1;j < atri;j += i , cnt ++){
                type0[j] = type0[j] * inv[cnt] % mod;
                type1[j] = type1[j] * inv[cnt] % mod;
            }
        }
    }
    // cout << "cnt = " << cnt << ' ' << "primemax = " << prime[cnt] <<'
';
    for(int i = 1;i < atri;i ++){
        // if(i == 46341) cout << i << '
';
        fac_0[i] = fac_0[i - 1] * type0[i] % mod;
        type1[i] = quick_pow(type1[i] , i * i % (mod - 1));
        fac_1[i] = fac_1[i - 1] * type1[i] % mod;
        // if(i > 46340)cout << "i = " << i << " " << "fac_0[i] = " << fac_0[i] << ' ' << "fac_0_inv[i] = " << fac_0_inv[i] << ' ' << "fac_1[i] = " << fac_1[i] << ' ' << "fac_1_inv[i] = " << fac_1_inv[i] << '
';
    }
    fac_0_inv[kato - 1] = quick_pow(fac_0[kato - 1] , mod - 2);
    fac_1_inv[kato - 1] = quick_pow(fac_1[kato - 1] , mod - 2);
    // cout << fac_0_inv[kato - 1] << ' ' << fac_1_inv[kato - 1];
    for(int i = kato - 2 ; i ; i --){
        fac_0_inv[i] = fac_0_inv[i + 1] * type0[i + 1] % mod;
    }
    for(int i = kato - 2 ; i ; i --){
        fac_1_inv[i] = fac_1_inv[i + 1] * type1[i + 1] % mod;
    }
    // cout << "QAQ" << '
';
}

inline int get_type2_(int x , int y){
    int ans = 1;
    for(int l = 1 , r = 0;l <= min(x , y);l = r + 1){
        r = min(x / (x / l) , y / (y / l));
        ans = ans * quick_pow(1LL * fac_0[r] * fac_0_inv[l - 1] % mod , 1LL * (x / l) * (y / l) % (mod - 1)) % mod; 
    }
    return ans % mod;
}

inline int get_sum(int a){
    return ((a * (a + 1))>> 1) % (mod - 1);
}

inline int get_type0_f1(int x , int y , int z){
    return quick_pow(fac[x] , y * z % (mod - 1)) % mod;
}

inline int get_type0_f2(int x , int y , int z){
    // cout << x << ' ' << y << ' ' << z << '
';
    int ans = 1;
    for(int l = 1 , r = 0;l <= min(x , y);l = r + 1){
        r = min(x / (x / l) , y / (y / l));
        // cout << fac_0[r] << ' ' << fac_0_inv[l - 1] << '
' ;
        ans = ans * quick_pow(fac_0[r] * fac_0_inv[l - 1] % mod , (x / l) * (y / l) % (mod - 1)) % mod; 
    }
    // cout << ans << '
';
    return quick_pow(ans , z);
}

inline int get_type1_f1(int x , int y , int z){
    return quick_pow(Fac[x] , get_sum(y) * get_sum(z) % (mod - 1)) % mod;
}

inline int get_type1_f2(int x , int y, int z){
    int ans = 1;
    for(int l = 1 , r = 0;l <= min(x , y);l = r + 1){
        r = min(x / (x / l) , y / (y / l));
        ans = ans * quick_pow(fac_1[r] * fac_1_inv[l - 1] % mod , get_sum(x / l) * get_sum(y / l) % (mod - 1)) % mod;
    }
    return quick_pow(ans , get_sum(z));
}

inline int get_type2_f1(int x , int y , int z){
    int ans = 1;
    for(int l = 1 , r = 0;l <= min(x , min(y , z));l = r + 1){
        r = min(x / (x / l) , min(y / (y / l) , z / (z / l)));
        ans = ans * quick_pow(fac[x / l] % mod , (phi_sum[r] - phi_sum[l - 1] + mod - 1) % (mod - 1) * (y / l) % (mod - 1) * (z / l) % (mod - 1)) % mod;
    }
    return ans % mod;
}

inline int get_type2_f2(int x , int y , int z){
    int ans = 1;
    for(int l = 1 , r = 0;l <= min(x , min(y , z));l = r + 1){
        r = min(x / (x / l) , min(y / (y / l) , z / (z / l)));
        ans = ans * quick_pow(get_type2_(x / l , y / l) , (phi_sum[r] - phi_sum[l - 1] + mod - 1) * (z / l) % (mod - 1)) % mod;
    }
    return ans % mod;
}

inline int get_type0(){
    int res1 = get_type0_f1(a , b , c) * get_type0_f1(b , a , c) % mod;
    int res2 = get_type0_f2(a , b , c) * get_type0_f2(a , c , b) % mod;
    // cout << res1 << ' ' << res2 << '
';
    return (res1 * quick_pow(res2 , mod - 2)) % mod;
}

inline int get_type1(){
    int res1 = get_type1_f1(a , b , c) * get_type1_f1(b , a , c) % mod;
    int res2 = get_type1_f2(a , b , c) * get_type1_f2(a , c , b) % mod;
    return (res1 * quick_pow(res2 , mod - 2)) % mod;
}

inline int get_type2(){
    int res1 = get_type2_f1(a , b , c) * get_type2_f1(b , a , c) % mod;
    int res2 = get_type2_f2(a , b , c) * get_type2_f2(a , c , b) % mod;
    return (res1 * quick_pow(res2 , mod - 2)) % mod;
}

inline int Ame_(){
    read(t) , read(mod);
    inv[1] = 1;
	for(int i = 2;i < kato;i ++){
        inv[i] = (mod - mod / i) * inv[mod % i] % mod;
    }
    mu[1] = phi[1] = phi_sum[0] = fac[0] = fac[1] = Fac[0] = Fac[1] = fac_0[0] = fac_0[1] = fac_1[0] = fac_1[1] = fac_0_inv[0] = fac_1_inv[0] = 1;
    pri();
    // cout << "QAQ" << '
';
    // for(int i = 1;i <= 10;i ++){
    //     cerr << i << ' ' << mu[i] << ' ' << phi[i] << '
';
    // }
    // cout << type0[1] << ' ' << type0[2] << ' ' << type0[3] << '
';
    // for(int i = 1;i <= 10;i ++){
    //     cerr << fac_0_inv[i] << ' ' << fac_1_inv[i] << '
';
    // }
    for(; t --> 0 ;){
        read(a) , read(b) , read(c);
        printf("%lld %lld %lld
" , get_type0() , get_type1() , get_type2());
        // put(get_type0());cout << ' ';put(get_type1());cout << ' ';put(get_type2());cout << '
';
    }
    return ~~(0^_^0);
}
    
int Ame__ = Ame_();
    
signed main(){;}

/*
3 998244853
1 1 1
2 2 2
3 3 3
*/

/*
1 1 1
16 4096 16
180292630 873575259 180292630
*/