旧试题解题报告

_LHF_ · 2021-03-03 12:24:42 · 题解

警告：此题卡常

花了我一天半的时间才弄好。

做法：纯正的莫比乌斯反演。

\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk)

先考虑后面那个函数

d(ijk)=\sum_{x|i}\sum_{y|j}\sum_{z|k}[x\perp y][y\perp z][x\perp z]

证明，见 [SDOI2015约数个数和] 好了，带回去得到原式得到

\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x|i}\sum_{y|j}\sum_{z|k}[x\perp y][y\perp z][x\perp z]

然后枚举 x,y,z 得到

\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C\lfloor\frac{A}{x}\rfloor\lfloor\frac{B}{y}\rfloor\lfloor\frac{C}{z}\rfloor[x\perp y][y\perp z][x\perp z]

考虑拎出来其中一个进行反演，得到

\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C\sum_{d|x,d|y}\mu(d)\lfloor\frac{A}{x}\rfloor\lfloor\frac{B}{y}\rfloor\lfloor\frac{C}{z}\rfloor[y\perp z][x\perp z]

钦定 A\le B

\sum_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum_{d=1}^A\mu(d)\sum_{x|d}^A\lfloor\frac{A}{x}\rfloor[y\perp z]\sum_{y|d}^B\lfloor\frac{B}{y}\rfloor[x\perp z]

设 p=\lfloor\frac{A}{d}\rfloor,q=\lfloor\frac{B}{d}\rfloor

\sum_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum_{d=1}^A\mu(d)\sum_{x=1}^p\lfloor\frac{A}{x}\rfloor[yd\perp z]\sum_{y|d}^B\lfloor\frac{B}{q}\rfloor[xd\perp z]

\sum_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum_{d=1}^A\mu(d)[d\perp z]\sum_{x=1}^p\lfloor\frac{A}{x}\rfloor[y\perp z]\sum_{y|d}^B\lfloor\frac{B}{q}\rfloor[x\perp z]

设 f(x,y)=\sum_{i=1}^x\lfloor\frac{x}{i}\rfloor[i\perp y]

\sum_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum_{d=1}^A\mu(d)[d\perp z]f(p,z)f(q,z)

再设 g(x,y)=\sum_{i=1}^x\mu(i)[i\perp y]

\sum_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum_{l,r}(g(r,z)-g(l-1,z))f(p,z)f(q,z)

其中后面是一个数论分块。

现在考虑计算 f(x,y) 和 g(x,y) ，设 p 是 y 的一个质因数（和上文的 p 无关）

那么考虑 y→y/p 之后的答案，再用容斥原理减回去

f(x,y)=\sum_{i=1}^x\lfloor\frac{x}{i}\rfloor[i\perp y/p]-\sum_{i=1}^x\lfloor\frac{x}{i}\rfloor[i\perp y/p][i|p]

=f(x,y/p)-\sum_{i=1}^{\lfloor\frac{x}{p}\rfloor}\lfloor\frac{x}{ip}\rfloor[ip\perp y/p]

=f(x,y/p)-[p\perp y/p]\sum_{i=1}^{\lfloor\frac{x}{p}\rfloor}\lfloor\frac{x}{ip}\rfloor[i\perp y/p]

=f(x,y/p)-[p\perp y/p]f(\lfloor\frac{x}{p}\rfloor,y/p)

然后是 g(x,y) ，同样的套路，设 p 是 y 的一个质因数（也和上文的 p 无关）

g(x,y)=\sum_{i=1}^x\mu(i)[i\perp y/p]-\sum_{i=1}^x\mu(i)[i\perp y/p][i|p]

=g(x,y/p)-\sum_{i=1}^{\lfloor\frac{x}{p}\rfloor}\mu(ip)[ip\perp y/p]

然后根据 \mu(ab)=\mu(a)\mu(b)[a\perp b] 以及 \mu(p)=-1 可得

=g(x,y/p)+\sum_{i=1}^{\lfloor\frac{x}{p}\rfloor}\mu(i)[i\perp p][ip\perp y/p]

=g(x,y/p)+[p\perp y/p]\sum_{i=1}^{\lfloor\frac{x}{p}\rfloor}\mu(i)[i\perp p][i\perp y/p]

=g(x,y/p)+[p\perp y/p]\sum_{i=1}^{\lfloor\frac{x}{p}\rfloor}\mu(i)[i\perp y]

=g(x,y/p)+[p\perp y/p]g(\lfloor\frac{x}{p}\rfloor,y)

然后…… 可以打出一种暴力，时间复杂度 O(n\sqrt n) ，这里假设 n=\max(A,B,C)空间复杂度也是 O(n\sqrt n) ，可以拿到 70 分。

#include<cstdio>
#include<algorithm>
#define N 100010
using namespace std;
int pp[N],C,lx[N];
const int mod=1e9+7;
struct node{
    int f[640][N],id[N],d[N],len,l,r,i,j,x;
    void work(int n,int*h)
    {
        len=0;
        for(l=1;l<=n;l=r+1)
        {
            r=n/(n/l);
            id[r]=++len;
            d[len]=r;
        }
        for(i=1;i<=len;i++)
        {
            x=d[i];
            f[i][1]=(h[x]+mod)%mod;
            for(j=2;j<=C;j++)
                f[i][j]=(f[i][j/lx[j]]-pp[j]*f[id[x/lx[j]]][j/lx[j]]+mod)%mod;
        }
    }
    inline int calc(int x,int y)
    { return x?f[id[x]][y]:0; }
}f1,f2;
struct node2{
    int f[1270][N],id[N],d[N],len,l,r,i,j,x;
    void work(int n,int m,int*h)
    {
        len=0;
        x=min(n,m);
        for(l=1;l<=x;l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            id[r]=++len;
            d[len]=r;
        }
        for(i=1;i<=len;i++)
        {
            x=d[i];
            f[i][1]=(h[x]+mod)%mod;
            for(j=2;j<=C;j++)
                f[i][j]=(f[i][j/lx[j]]+pp[j]*f[id[x/lx[j]]][j])%mod;
        }
    }
    inline int calc(int x,int y)
    { return x?f[id[x]][y]:0; }
}g;
int vis[N],p[N],mu[N],h[N],T,id[N],d[N],n,A,B,x,k[N],len;
long long ans,s1,s2,s3,s4,s5,s6,s7,s8;
signed main()
{
    int i,j,l,r;
    n=N-10;
    mu[1]=k[1]=h[1]=1;
    for(i=2;i<=n;i++)
    {
        h[i]++;
        if(!vis[i]) p[++len]=i,mu[i]=-1,lx[i]=i;
        for(j=1;j<=len&&i*p[j]<=n;j++)
        {
            lx[i*p[j]]=p[j];
            vis[i*p[j]]=1;
            if(i%p[j]==0) break;
            mu[i*p[j]]=-mu[i];
        }
        for(j=i;j<=n;j+=i) h[j]++;
        k[i]=k[i-1]+mu[i];
        pp[i]=(i%(1ll*lx[i]*lx[i])!=0);
    }
    for(i=2;i<=n;i++) h[i]+=h[i-1];
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d",&A,&B,&C);
        ans=s1=s2=s3=s4=s5=s6=s7=s8=0;
        f1.work(A,h);
        f2.work(B,h);
        x=min(A,B);
        g.work(A,B,k);
        x=min(A,B);
        //循环展开
        for(i=1;i+8<=C;i+=8)
        {
            for(l=1;l<=min(A,B);l=r+1)
            {
                r=min(A/(A/l),B/(B/l));
                s1+=1ll*(C/i)*(g.calc(r,i)-g.calc(l-1,i)+mod)%mod*f1.calc(A/r,i)%mod*f2.calc(B/r,i)%mod;
                s2+=1ll*(C/(i+1))*(g.calc(r,i+1)-g.calc(l-1,i+1)+mod)%mod*f1.calc(A/r,i+1)%mod*f2.calc(B/r,i+1)%mod;
                s3+=1ll*(C/(i+2))*(g.calc(r,i+2)-g.calc(l-1,i+2)+mod)%mod*f1.calc(A/r,i+2)%mod*f2.calc(B/r,i+2)%mod;
                s4+=1ll*(C/(i+3))*(g.calc(r,i+3)-g.calc(l-1,i+3)+mod)%mod*f1.calc(A/r,i+3)%mod*f2.calc(B/r,i+3)%mod;
                s5+=1ll*(C/(i+4))*(g.calc(r,i+4)-g.calc(l-1,i+4)+mod)%mod*f1.calc(A/r,i+4)%mod*f2.calc(B/r,i+4)%mod;
                s6+=1ll*(C/(i+5))*(g.calc(r,i+5)-g.calc(l-1,i+5)+mod)%mod*f1.calc(A/r,i+5)%mod*f2.calc(B/r,i+5)%mod;
                s7+=1ll*(C/(i+6))*(g.calc(r,i+6)-g.calc(l-1,i+6)+mod)%mod*f1.calc(A/r,i+6)%mod*f2.calc(B/r,i+6)%mod;
                s8+=1ll*(C/(i+7))*(g.calc(r,i+7)-g.calc(l-1,i+7)+mod)%mod*f1.calc(A/r,i+7)%mod*f2.calc(B/r,i+7)%mod;
                s1%=mod,s2%=mod,s3%=mod,s4%=mod,s5%=mod,s6%=mod,s7%=mod,s8%=mod;
            }
        }
        for(;i<=C;i++)
        {
            for(l=1;l<=min(A,B);l=r+1)
            {
                r=min(A/(A/l),B/(B/l));
                ans+=1ll*(C/i)*(g.calc(r,i)-g.calc(l-1,i)+mod)%mod*f1.calc(A/r,i)%mod*f2.calc(B/r,i)%mod;
                ans%=mod;
            }
        }
        printf("%lld\n",(ans+s1+s2+s3+s4+s5+s6+s7+s8)%mod);
    }
}

然后你就发现最后几个点挂掉了。

空间太大了啊啊啊啊啊啊啊！！！

考虑再优化。

如果 y 含有某个质数的平方，那么 p 我们可以取该质数，那么 f(x,y)=f(x,y/p)

因为这样 [p\perp y/p]=0

所以我们可以弄一个 DFS ，来枚举每一个不含某质数的平方因子的数（也可以理解为 \mu(x)\not=0 的数），不仅时间减少了，空间更是少了不少。

#include<cstdio>
#define N 101010
#define ri register int
using namespace std;
const int mod=1e9+7;
int vis[N],p[N],mu[N],T,n,A,B,x,len;
long long g[N][12],f[N][12];
int k[N],h[N],lw[N],s[N],ans,C;
inline int min(int a,int b){return a<b?a:b;}
void dfs(ri x,ri i,ri k)
{
    ri ss=0,l,r;
    for(l=1;l<=A;l=r+1)
    {
        r=min(A/(A/l),B/(B/l));
        ss=(ss+(g[r][k]-g[l-1][k]+mod)*f[A/r][k]%mod*f[B/r][k]%mod)%mod;
    }
    ans=(ans+1ll*ss*s[x]%mod)%mod;
    for(i=i+1;1ll*p[i]*x<=C&&i<=len;i++)
    {
        for(l=1;l<=A;l=r+1)
        {
            r=min(A/(A/l),B/(B/l));
            g[r][k+1]=(g[r][k]+g[r/p[i]][k+1])%mod;
            f[A/r][k+1]=(f[A/r][k]-f[A/r/p[i]][k]+mod)%mod;
            f[B/r][k+1]=(f[B/r][k]-f[B/r/p[i]][k]+mod)%mod;
        }
        dfs(x*p[i],i,k+1);
    }
}
int main()
{
    n=N-1010;
    mu[1]=k[1]=h[1]=lw[1]=1;
    int i,j,l,r;
    for(i=2;i<=n;i++)
    {
        if(!vis[i]) p[++len]=lw[i]=i,mu[i]=-1,h[i]=2;
        for(j=1;j<=len&&i*p[j]<=n;j++)
        {
            vis[i*p[j]]=1;
            lw[i*p[j]]=lw[i]*p[j];
            h[i*p[j]]=h[i]*2;
            if(i%p[j]==0)
            {
                h[i*p[j]]-=h[i/p[j]];
                lw[i*p[j]]=lw[i];
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
        k[i]=k[i-1]+mu[i];
    }
    for(i=2;i<=n;i++) h[i]+=h[i-1];
    scanf("%d",&T);
    while(T--)
    {
        ans=0;
        scanf("%d%d%d",&A,&B,&C);
        if(A>B) i=A,A=B,B=i;
        for(l=1,r;l<=A;l=r+1)
        {
            r=min(A/(A/l),B/(B/l));
            g[r][1]=k[r];
            f[A/r][1]=h[A/r];
            f[B/r][1]=h[B/r];
        }
        for(i=1;i<=C;i++) s[i]=0;
        for(i=1;i<=C;i++) s[lw[i]]+=1ll*C/i,s[lw[i]]%=mod;
        dfs(1,0,1);
        printf("%lld\n",ans);
    }
}

卡常卡到极致了，但你发现还是TLE了……

某大佬告诉我，如果 a,b\in[0,p) ，那么可以直接自定义一个取模函数 int Mod(int a){return a<p?a:a-p;}

然后就 AC 了。

代码不放了（中国武功讲究点到为止）

旧试题 解题报告

警告：此题卡常

旧试题解题报告