我想时间|

Sya_Resory

园龄:6年1个月粉丝:3关注:8

【题解】Luogu P5518 题解

Description

令:

f(type)={1type=0ijktype=1gcd(i,j,k)type=2

给定 A,B,C,求:

anstype=i=1Aj=1Bk=1C(lcm(i,j)gcd(i,k))f(type)

多测,对给定质数 p 取模。

1A,B,C105107p1.05×109,pPT=70

Solution

首先对原式化简一下:

anstype=i=1Aj=1Bk=1C(lcm(i,j)gcd(i,k))f(type)=i=1Aj=1Bk=1C(ijgcd(i,j)gcd(i,k))f(type)

于是原式可以分为两个部分:

F(A,B,C)=i=1Aj=1Bk=1Cif(type)G(A,B,C)=i=1Aj=1Bk=1Cgcd(i,j)f(type)

显然有:

anstype=F(A,B,C)F(B,A,C)G(A,B,C)G(A,C,B)

接下来就开始推柿子吧!


type=0

F(A,B,C)=i=1Aj=1Bk=1Ci=i=1Aj=1BiC=i=1AiBC=(A!)BC

预处理阶乘,每次计算快速幂即可。

G(A,B,C)=i=1Aj=1Bk=1Cgcd(i,j)=(i=1Aj=1Bgcd(i,j))C

然后开始莫反:

i=1Aj=1Bgcd(i,j)=d=1i=1Aj=1Bd[gcd(i,j)=d]=d=1i=1Adj=1Bdd[gcd(i,j)=1]=d=1di=1A/dj=1B/d[gcd(i,j)=1]

把指数提出来推一下:

i=1Adj=1Bd[gcd(i,j)=1]=e=1i=1Adj=1Bd[ei][ej]μ(e)=e=1i=1Adej=1Bdeμ(e)=e=1AdeBdeμ(e)

带回去就有:

i=1Aj=1Bgcd(i,j)=d=1de=1AdeBdeμ(e)=t=1dtdAtBtμ(td)=t=1(dtdμ(td))AtBtG(A,B,C)=i=1Aj=1Bk=1Cgcd(i,j)=(t=1(dtdμ(td))AtBt)C

最里面括号内的东西可以用类似埃氏筛的方法预处理一下,即对于每个 d,将它的贡献累乘到它的每个倍数上。预处理完成后就可以快乐的整除分块了。


type=1

定义:S(n)=i=1ni=n(n+1)2

然后开始推柿子:

F(A,B,C)=i=1Aj=1Bk=1Ciijk=(i=1Aj=1Biij)S(C)=(i=1Aii)S(B)S(C)

预处理 i=1Aii 即可。

G(A,B,C)=i=1Aj=1Bk=1Cgcd(i,j)ijk=(i=1Aj=1Bgcd(i,j)ij)S(C)i=1Aj=1Bgcd(i,j)ij=d=1i=1Aj=1Bd[gcd(i,j)=d]ij=d=1i=1Adj=1Bdd[gcd(i,j)=1]ijd2=d=1dd2i=1A/dj=1B/d[gcd(i,j)=1]ij

把指数拎出来:

i=1Adj=1Bd[gcd(i,j)=1]ij=e=1i=1Adj=1Bd[ei][ej]μ(e)ij=e=1i=1Adej=1Bdee2μ(e)ij=e=1e2μ(e)(i=1Adej=1Bdeij)=e=1e2μ(e)S(Ade)S(Bde)

带回去:

i=1Aj=1Bgcd(i,j)ij=d=1dd2e=1e2μ(e)S(Ade)S(Bde)=e=1d=1d(de)2μ(e)S(Ade)S(Bde)=t=1dtdt2μ(td)S(At)S(Bt)=t=1(dtdμ(td))t2S(At)S(Bt)G(A,B,C)=(i=1Aj=1Bgcd(i,j)ij)S(C)=((t=1(dtdμ(td))t2)S(At)S(Bt))S(C)

括号里的东西其实就是第一问预处理的那个东西的 t2 次方,预处理一下即可整除分块。


type=2

毒瘤。

F(A,B,C)=i=1Aj=1Bk=1Cigcd(i,j,k)=i=1Aij=1Bk=1Cgcd(i,j,k)=d=1i=1Aij=1Bk=1C[gcd(i,j,k)=d]d=d=1i=1Ad(id)j=1B/dk=1C/d[gcd(i,j,k)=1]d=d=1e=1i=1Ade(ide)BdeCdeμ(e)d=t=1dti=1At(it)BtCtμ(td)d=t=1(dt(At!tAt)μ(td)d)BtCt=t=1((At!tAt)dtμ(td)d)BtCt=t=1((At!tAt)φ(t))BtCt

这玩意有点复杂,把它拆成两个部分:

f0(A,B,C)=t=1((At!)BtCt)φ(t)f1(A,B,C)=t=1(tφ(t))AtBtCt

f1 可以整除分块了,把 f2 先放着,看 G

G(A,B,C)=i=1Aj=1Bk=1Cgcd(i,j)gcd(i,j,k)=d=1i=1Aj=1Bk=1Cd[gcd(i,j)=d]gcd(d,k)=d=1di=1A/dj=1B/d[gcd(i,j)=1]k=1Cgcd(d,k)=d=1de=1AdeBdeμ(e)k=1Cgcd(d,k)=t=1(dtdμ(td)k=1Cgcd(d,k))AtBt

指数上的东西:

k=1Cgcd(d,k)=fdk=1Cf[gcd(df,k)=1]f=fdgdfCfgμ(g)f=udφ(u)Cu

然后带回去:

G(A,B,C)=t=1(dtdμ(td)udφ(u)Cu)AtBt

这玩意不怎么好预处理。考虑把底数的 d 拆成 udu

g0(A,B,C)=t=1(dtud(du)μ(td)φ(u)Cu)AtBtg1(A,B,C)=t=1(dtuduμ(td)φ(u)Cu)AtBt

分别推一下:

g0(A,B,C)=t=1(dtud(du)μ(td)φ(u)Cu)AtBt=t=1(utudt(du)μ(td)φ(u)Cu)AtBt=t=1(ut(dtudμ(td))φ(u)Cu)AtBt=u=1(t=1(dtdμ(td))AtuBtu)φ(u)Cu

括号里的东西预处理过了,两次整除分块可以搞定。

g1(A,B,C)=t=1(dtuduμ(td)φ(u)Cu)AtBt=u=1t=1(dtuμ(td)φ(u)Cu)AtuBtu=u=1ut=1dtμ(td)φ(u)CuAtuBtu=u=1uφ(u)Cut=1AtuBtudtμ(td)=u=1uφ(u)Cut=1AtuBtu[t=1]=u=1uφ(u)AuBuCu

注意到,g1 和之前的 f1 其实是一样的,而且它们的值不受 A,B,C 顺序的影响,因此可以把它们约掉。


卡常小技巧(不过应该众所周知):把预处理的东西的逆元同时预处理出来可以大幅减小常数

code:

#include <cstdio>

const int N = 1e5,maxn = N + 5;

int T,p,phi,A,B,C,ta,tb,tc;
int pri[maxn],mu[maxn],Inv[maxn],f[maxn],Phi[maxn];
int fac[maxn],fd[maxn],facInv[maxn],fdInv[maxn];
int Fac[maxn],Fd[maxn],FdInv[maxn],ps[maxn];
bool isp[maxn];

inline int read() {
#define gc c = getchar()
    int d = 0,f = 0,gc;
    for(;c < 48 || c > 57;gc) f |= (c == '-');
    for(;c > 47 && c < 58;gc) d = (d << 1) + (d << 3) + (c ^ 48);
#undef gc
    return f ? -d : d;
}

inline int Mul(int a,int b,int mod) { return 1LL * a * b % mod; }
inline int Add(int a,int b,int mod) { a += b; return a > mod ? a - mod : a; }
inline int max(int a,int b) { return a > b ? a : b; }
inline int min(int a,int b) { return a < b ? a : b; }
inline int min(int a,int b,int c) { return min(min(a,b),c); }

inline int fpow(int a,int b,int mod) {
    int res = 1;
    for(;b;a = Mul(a,a,mod),b >>= 1) if(b & 1) res = Mul(res,a,mod);
    return res;
}

inline int GetInv(int a) { return a <= N ? Inv[a] : fpow(a,p - 2,p); }

namespace Type0 {
    inline int F(int A,int B,int C) { return fpow(fac[A],Mul(B,C,phi),p); }

    inline int G(int A,int B,int C) {
        int res,ans = 1;
        for(int r,l = 1;l <= A && l <= B;l = r + 1) {
            ta = A / l,tb = B / l; r = min(A / ta,B / tb);
            res = Mul(fd[r],fdInv[l - 1],p);
            ans = Mul(ans,fpow(res,Mul(ta,tb,phi),p),p);
        }
        return fpow(ans,C,p);
    }

    inline int solve() {
        int res = Mul(F(A,B,C),F(B,A,C),p);
        int inv = Mul(G(A,B,C),G(A,C,B),p);
        return Mul(res,GetInv(inv),p);
    }
}

namespace Type1 {
    inline int S(int n) { return 1LL * n * (n + 1) / 2 % phi; }

    inline int F(int A,int B,int C) { return fpow(Fac[A],Mul(S(B),S(C),phi),p); }

    inline int G(int A,int B,int C) {
        int res,ans = 1;
        for(int r,l = 1;l <= A && l <= B;l = r + 1) {
            ta = A / l,tb = B / l;
            r = min(A / ta,B / tb);
            res = Mul(Fd[r],FdInv[l - 1],p);
            ans = Mul(ans,fpow(res,Mul(S(ta),S(tb),phi),p),p);
        }
        return fpow(ans,S(C),p);
    }

    inline int solve() {
        int res = Mul(F(A,B,C),F(B,A,C),p);
        int inv = Mul(G(A,B,C),G(A,C,B),p);
        return Mul(res,GetInv(inv),p);
    }
}

namespace Type2 {
    inline int F(int A,int B,int C) { 
        int exp,res,ans = 1;
        for(int r,l = 1;l <= min(A,B,C);l = r + 1) {
            ta = A / l,tb = B / l,tc = C / l;
            r = min(A / ta,B / tb,C / tc);
            exp = Add(ps[r],phi - ps[l - 1],phi);
            res = fpow(fac[ta],Mul(tb,tc,phi),p);
            ans = Mul(ans,fpow(res,exp,p),p);
        }
        return ans;
    }

    inline int g0(int A,int B) {
        int res,ans = 1;
        for(int r,l = 1;l <= A && l <= B;l = r + 1) {
            ta = A / l,tb = B / l;
            r = min(A / ta,B / tb);
            res = Mul(fd[r],fdInv[l - 1],p);
            ans = Mul(ans,fpow(res,Mul(ta,tb,phi),p),p);
        }
        return ans;
    }

    inline int G(int A,int B,int C) {
        int res,ans = 1;
        for(int r,l = 1;l <= min(A,B,C);l = r + 1) {
            ta = A / l,tb = B / l,tc = C / l;
            r = min(A / ta,B / tb,C / tc);
            res = fpow(g0(ta,tb),tc,p);
            ans = Mul(ans,fpow(res,Add(ps[r],phi - ps[l - 1],phi),p),p);
        }
        return ans;
    }

    inline int solve() {
        int res = Mul(F(A,B,C),F(B,A,C),p);
        int inv = Mul(G(A,B,C),G(A,C,B),p);
        return Mul(res,GetInv(inv),p);
    }
}

inline void Init() {
    mu[1] = Phi[1] = f[1] = Inv[0] = Inv[1] = 1;
    fac[0] = fac[1] = Fac[0] = facInv[0] = 1;
    fd[0] = Fd[0] = fdInv[0] = FdInv[0] = 1;
    for(int i = 2;i <= N;i ++) {
        f[i] = 1,fac[i] = Mul(i,fac[i - 1],p);
        if(!isp[i]) pri[++ pri[0]] = i,mu[i] = -1,Phi[i] = i - 1;
        for(int j = 1;j <= pri[0] && i * pri[j] <= N;j ++) {
            isp[i * pri[j]] = true;
            if(!(i % pri[j])) {
                mu[i * pri[j]] = 0;
                Phi[i * pri[j]] = Phi[i] * pri[j];
                break;
            }
            mu[i * pri[j]] = -mu[i];
            Phi[i * pri[j]] = Phi[i] * Phi[pri[j]];
        }
    }
    facInv[N] = fpow(fac[N],p - 2,p);
    for(int i = N - 1;i;i --) {
        facInv[i] = Mul(facInv[i + 1],i + 1,p);
        Inv[i + 1] = Mul(facInv[i + 1],fac[i],p);
    }
    for(int i = 1;i <= N;i ++) {
        ps[i] = Add(Phi[i],ps[i - 1],phi);
        Fac[i] = Mul(Fac[i - 1],fpow(i,i,p),p);
        fd[i] = Mul(f[i],fd[i - 1],p);
        Fd[i] = Mul(fpow(f[i],Mul(i,i,phi),p),Fd[i - 1],p);
        fdInv[i] = GetInv(fd[i]),FdInv[i] = GetInv(Fd[i]);
        if(!mu[i]) continue;
        for(int j = 1;i * j <= N;j ++)
           f[i * j] = Mul(f[i * j],mu[i] == 1 ? j : Inv[j],p);
    }
}

int main() {
    T = read(),p = read();
    phi = p - 1; Init();
    for(;T;T --) {
        A = read(),B = read(),C = read();
        printf("%d %d %d\n",Type0::solve(),Type1::solve(),Type2::solve());
    }
    return 0;
}

本文作者:Resory

本文链接:https://www.cnblogs.com/resorie/p/soln-lg-p5518.html

版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。

posted @   Sya_Resory  阅读(37)  评论(0编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起