P7496 干杯！再会！ 题解

· · 题解

n\le 8

直接枚举 p 然后计算方差即可。

写起来意外的麻烦

复杂度 O(n!)。

#include<bits/stdc++.h>
#define N 10
#define MOD 1000000007
using namespace std;
void read(int &x)
{
    x=0;
    char ch=getchar();
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9')
    {
        x=x*10+ch-'0';
        ch=getchar();
    }
}
int n,invn,ans;
int a[N],b[N];
int p[N],v[N];
int ksm(int x,int k)
{
    int res=1;
    while(k)
    {
        if(k&1)res=1ll*res*x%MOD;
        x=1ll*x*x%MOD;
        k>>=1; 
    }
    return res;
}
int gcd(int x,int y)
{
    if(x==0)return y;
    return gcd(y%x,x);
}
void dfs(int x)
{
    if(x>n)
    {
        int pj=0;
        for(int i=1;i<=n;i++)
        {
            int res=gcd(a[i],b[p[i]]);
            pj=(pj+res)%MOD;
        }
        pj=1ll*pj*invn%MOD;
        for(int i=1;i<=n;i++)
        {
            int res=gcd(a[i],b[p[i]]);
            ans=(ans+1ll*(res-pj)*(res-pj)%MOD)%MOD;
        }
        return;
    }
    for(int i=1;i<=n;i++)
        if(!v[i])
        {
            p[x]=i;
            v[i]=1;
            dfs(x+1);
            p[x]=0;
            v[i]=0;
        }
}
int main()
{
    scanf("%d",&n);
    invn=ksm(n,MOD-2);
    for(int i=1;i<=n;i++)
        read(a[i]);
    for(int i=1;i<=n;i++)
        read(b[i]);
    dfs(1);
    ans=1ll*ans*invn%MOD;
    printf("%d",ans);
}

n\le 100

如果你想要直接递推，比如说 dp_i 表示只考虑前 i 位的时候的答案，那么由于原题是个随机的排列，你就得知道前 i 位到底用了什么位置，必须要 2^n 来记录，肯定是做不了的。

所以我们把方差拆掉：\begin{aligned}\dfrac{1}{n}\sum\limits_{i=1}^{n}(x_i-\dfrac{S}{n})^2=\dfrac{1}{n}(\sum\limits_{i=1}^{n}x_i^2-\dfrac{S^2}{n})\end{aligned}

那么这题就变成了求两个东西：

1.数列中的数的平方的和 \begin{aligned}sum1=\sum\limits_{i=1}^{n}x_i^2\end{aligned}。

2.数列中的数的和的平方 sum2=S^2。

问题就是 $sum2$，没法独立考虑，因为这题的性质不允许我们 $dp$，所以就**直接把它拆掉**！ 具体地， $\begin{aligned}S^2=\sum\limits_{i=1}^{n}x_i^2+\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[i\ne j]x_ix_j\end{aligned}

那么我们记录两个值：

\begin{aligned}res1=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_i,b_j)^2\end{aligned} \begin{aligned}res2=\sum\limits_{i1=1}^{n}\sum\limits_{i2=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}[i1 \ne i2][j1 \ne j2]\gcd(a_{i1},b_{j1})\times \gcd(a_{i2},b_{j2})\end{aligned}

就会有 sum1=res1\times (n-1)!，sum2=res1\times (n-1)!+res2\times (n-2)!

最后答案就是 ans=\dfrac{1}{n}(sum1-\dfrac{sum2}{n})。

这样直接做就是 O(n^4) 的，可以拿 25 分。

下面这个代码真的是直接明了，就不写注释了。

#include<bits/stdc++.h>
#define N 110
#define MOD 1000000007
using namespace std;
void read(int &x)
{
    x=0;
    char ch=getchar();
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9')
    {
        x=x*10+ch-'0';
        ch=getchar();
    }
}
int n,invn,sum1,sum2,ans;
int a[N],b[N],mul[N];
int c[N][N];
int ksm(int x,int k)
{
    int res=1;
    while(k)
    {
        if(k&1)res=1ll*res*x%MOD;
        x=1ll*x*x%MOD;
        k>>=1;
    }
    return res;
}
int gcd(int x,int y)
{
    if(x==0)return y;
    return gcd(y%x,x);
}
int main()
{
    scanf("%d",&n);
    invn=ksm(n,MOD-2);
    for(int i=1;i<=n;i++)
        read(a[i]);
    for(int i=1;i<=n;i++)
        read(b[i]);
    mul[0]=1;
    for(int i=1;i<=n;i++)
        mul[i]=1ll*mul[i-1]*i%MOD;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            c[i][j]=gcd(a[i],b[j]);
            sum1=(sum1+1ll*c[i][j]*c[i][j])%MOD;
        }
    }
    sum1=1ll*sum1*mul[n-1]%MOD;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            for(int ii=1;ii<=n;ii++)
            {
                if(ii==i)continue;
                for(int jj=1;jj<=n;jj++)
                {
                    if(jj==j)continue;
                    int res1=c[i][j],res2=c[ii][jj];
                    sum2=(sum2+1ll*res1*res2)%MOD;
                }
            }
        }
    }
    sum2=1ll*sum2*mul[n-2]%MOD;
    sum2=(sum1+sum2)%MOD;
    ans=(sum1-1ll*invn*sum2%MOD+MOD)%MOD;
    ans=1ll*ans*invn%MOD;
    printf("%d",ans);
}

n\le 10^6

（O(n^2)都是些奇怪的东西，就跳过了。）

然后就是想办法加快求 res1,res2 的速度了。

首先是 res1，非常经典的 \gcd 容斥：（知道的可以跳过了）

我们记录 f_x 表示 \gcd(a,b)=x 的 (a,b) 有多少组。

也就是：\begin{aligned}f_x=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[\gcd(a_i,b_j)=x]\end{aligned}。

不过想要 \gcd(a,b)=x 的话，至少得有 x|a,x|b 吧。

也就是： \begin{aligned}f_x=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[x|a_i][x|b_j]\end{aligned}

但是这样的话如果 \gcd(a,b)=kx 也会被算进去，所以要减去 f_{2x},f_{3x},f_{4x}... 就把错的去掉了。

最终会有：\begin{aligned}f_x=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[x|a_i][x|b_j]-\sum\limits_{i=2}^{\left\lfloor\frac{n}{x}\right\rfloor}f_{ix}\end{aligned}

直接预处理 \begin{aligned}ca_x=\sum\limits_{i=1}^{n}[x|a_i]\end{aligned}，\begin{aligned}cb_x=\sum\limits_{i=1}^{n}[x|b_i]\end{aligned}，就可以从大往小求出 f_i 了。

最后 \begin{aligned}res1=\sum\limits_{i=1}^{n}f_i\times i^2\end{aligned}

复杂度为 O(\left\lfloor\dfrac{n}{1}\right\rfloor+\left\lfloor\dfrac{n}{2}\right\rfloor+...\left\lfloor\dfrac{n}{n}\right\rfloor)=O(n\log n)。

然后就是本题的难点所在，res2 该怎么求？

仔细观察原式：

\begin{aligned}res2=\sum\limits_{i1=1}^{n}\sum\limits_{i2=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}[i1 \ne i2][j1 \ne j2]\gcd(a_{i1},b_{j1})\times \gcd(a_{i2},b_{j2})\end{aligned}

其中要满足两个条件 [i1 \ne i2][j1 \ne j2]。

但是在计数的时候保证选的两个不同的数是难的，而选两个相同的数或者是任选两个数都是好做的。

所以这就启发我们去容斥！

具体来说我们先求出无限制下的答案，然后减去满足 [i1=i2]（注意这里是相等）而 (j1,j2) 无限制的答案，再减去 (i1,i2) 无限制而满足 [j1=j2] 的答案，最后还要加回来满足 [i1=i2] 且 [j1=j2] 的答案。

有点繁琐，不过意会一下应该能明白。

我们把上面四个情况顺次记作 tot1,tot2,tot3,tot4，下面我们一个个求。

tot1

先把式子写出来：

\begin{aligned}tot1=\sum\limits_{i1=1}^{n}\sum\limits_{i2=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}\gcd(a_{i1},b_{j1})\times \gcd(a_{i2},b_{j2})\end{aligned}

观察一下发现：

\begin{aligned}tot1=(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j}))^2\end{aligned}

想到之前求 res1 时用的 f_i 了吗？借用一下即可。

\begin{aligned}tot1=(\sum\limits_{i=1}^{n}f_i\times i)^2\end{aligned}

tot2,tot3

因为它俩一模一样，所以就讲讲 tot2 咋求。

\begin{aligned}tot2=\sum\limits_{i=1}^{n}\sum\limits_{j1=1}^{n}\sum\limits_{j2=1}^{n}\gcd(a_{i},b_{j1})\times \gcd(a_{i},b_{j2})\end{aligned}

（因为 i1=i2，所以就直接改成 i 了）

再观察一下可以发现 j1,j2 能用一个平方弄掉：

\begin{aligned}tot2=\sum\limits_{i=1}^{n}(\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j}))^2\end{aligned}

然后就用到我们的欧拉反演了！

当然，反演的式子很简单：\begin{aligned}n=\sum\limits_{d|n}\varphi(d)\end{aligned}。

咋用呢，其实就是把gcd换掉！

\begin{aligned}\gcd(a_i,b_j)=\sum\limits_{d|\gcd(a_i,b_j)}\varphi(d)\end{aligned}

代入即有：

tot2&= \sum\limits_{i=1}^{n}(\sum\limits_{j=1}^{n}\sum\limits_{d|\gcd(a_i,b_j)}\varphi(d))^2\\ &=\sum\limits_{i=1}^{n}(\sum\limits_{j=1}^{n}\sum\limits_{d|a_i,d|b_j}\varphi(d))^2\\ &=\sum\limits_{i=1}^{n}(\sum\limits_{d|a_i}\sum\limits_{j=1}^{n}[d|b_j]\varphi(d))^2\\ \end{aligned}

我们记录一个辅助数组 \begin{aligned}gb_x=\sum\limits_{i=1}^{n}[x|b_i]\varphi(x)\end{aligned}。

那么上式就是：\begin{aligned}tot2= \sum\limits_{i=1}^{n}(\sum\limits_{d|a_i}gb_d)^2\end{aligned}

这样就可以 O(n\log n) 啦。

如果你不用欧拉反演而是继续用之前的容斥的话，复杂度就是 O(nlog^2n) 的，不过仔细想想把这个容斥的过程预处理出来，其实就是欧拉函数！

还是挺简单的。

tot4

\begin{aligned}tot4=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j})^2=res1\end{aligned}

没了。

然后根据 res2=tot1-tot2-tot3+tot4，这题就做完了。

复杂度 O(n\log n)。

写代码的时候还要注意一个细节，就是所有的 a 我们是要存进桶中的，这样不同的数才会去计算一遍，否则同一个数计算多次的话复杂度是错的。

代码长度也就暴力的两倍

#include<bits/stdc++.h>
#define N 1000010
#define MOD 1000000007
using namespace std;
void read(int &x)
{
    x=0;
    char ch=getchar();
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9')
    {
        x=x*10+ch-'0';
        ch=getchar();
    }
}
int n,m,invn,sum1,sum2,ans;
int a[N],b[N],mul[N];
int ca[N],cb[N],f[N],phi[N],ga[N],gb[N];
int ksm(int x,int k)
{
    int res=1;
    while(k)
    {
        if(k&1)res=1ll*res*x%MOD;
        x=1ll*x*x%MOD;
        k>>=1; 
    }
    return res;
}
int gcd(int x,int y)
{
    if(x==0)return y;
    return gcd(y%x,x);
}
int main()
{
    scanf("%d",&n);
    invn=ksm(n,MOD-2);
    for(int i=1;i<=n;i++)
    {
        int x;
        read(x);
        ++a[x];//要存进桶里 
        m=max(m,x);
    }
    for(int i=1;i<=n;i++)
    {
        int x;
        read(x);
        ++b[x];
        m=max(m,x);
    }
    mul[0]=1;
    for(int i=1;i<=n;i++)
        mul[i]=1ll*mul[i-1]*i%MOD;//前缀积 
    for(int i=m;i>=1;i--)
        for(int j=1;j*i<=m;j++)
            ca[i]+=a[j*i];//把多少个a是i的倍数算进ca[i]中 
    for(int i=m;i>=1;i--)
        for(int j=1;j*i<=m;j++)
            cb[i]+=b[j*i];
    for(int i=m;i>=1;i--)
    {
        f[i]=1ll*ca[i]*cb[i]%MOD;
        for(int j=2;j*i<=m;j++)
            (f[i]+=MOD-f[j*i])%=MOD;//容斥,求出有多少gcd(a,b)=i,存进f[i] 
        sum1=(sum1+1ll*f[i]*i%MOD*i)%MOD;//其实是res1
        sum2=(sum2+1ll*f[i]*i)%MOD;//tot1
    }
    sum2=1ll*sum2*sum2%MOD;

    for(int i=1;i<=m;i++)
    {
        phi[i]=(phi[i]+i)%MOD;
        //这也是phi的一种求法，这样写主要是为了体现出其容斥的内涵 
        int tmpa=1ll*phi[i]*ca[i]%MOD,tmpb=1ll*phi[i]*cb[i]%MOD;
        (gb[i]+=tmpb)%=MOD;
        (ga[i]+=tmpa)%=MOD;
        sum2=(sum2-1ll*gb[i]*gb[i]%MOD*a[i]%MOD+MOD)%MOD;//tot2
        sum2=(sum2-1ll*ga[i]*ga[i]%MOD*b[i]%MOD+MOD)%MOD;//tot3
        for(int j=2;j*i<=m;j++)
            (phi[j*i]+=MOD-phi[i])%=MOD,(gb[j*i]+=tmpb)%=MOD,(ga[j*i]+=tmpa)%=MOD;
    }

    sum2=(sum2+sum1)%MOD;//tot4

    sum1=1ll*sum1*mul[n-1]%MOD;
    sum2=1ll*sum2*mul[n-2]%MOD;
    sum2=(sum1+sum2)%MOD;
    ans=(sum1-1ll*invn*sum2%MOD+MOD)%MOD;
    ans=1ll*ans*invn%MOD;
    printf("%d",ans);
}

n\le 5\times 10^6

没有被放进去的数据范围。

这咋整啊！卡常？

其实是高维前缀和优化。

复杂度上看，高维前缀和的复杂度是 \sum\limits_{p}\left\lfloor\dfrac{n}{p}\right\rfloor （p 是质数）=O(n\log \log n)，卡5000000还是很轻松的。

那么高维前缀和优化到底是啥呢？

其实就是我们想要求一个东西 \begin{aligned}f_i=\sum\limits_{j=1}^{ij\le n}g_{ij}\end{aligned}，其中 g 是已知的。

那么我们可以先把 j=2 的求出来，再加上 j=3 的，再加上 j=5 的......

而在加的过程中我们可以把之前算过的 f 拿来用，跟前缀和一样。而且为了不算重，枚举的 j 都是质数。

而这可以看作是一个多维的空间，其中每一维都是由一个质数 p 的倍数构成的，因此原理其实就是高维前缀和。

代码也很短，如下：

for(int i=1;i<=n;i++)f[i]=g[i];

for(int i=1;i<=p[0];i++)
    for(int j=n/p[i];j>=1;j--)
            f[j]+=f[j*p[i]];

那么我们原来的求 ca[i] 的代码就可以改掉了，但是有些是无法改的！

我们举例一下原来所有要 O(n\log n) 求的东西现在可以做到的复杂度：

1.phi[i]： 直接线性筛。O(n)

2.ca[i],cb[i]：高维前缀和。O(n\log\log n)

3.ga[i],gb[i]：高维前缀和。O(n\log\log n)

4.f[i]：因为要由自己容斥出自己，所以无法高维前缀和优化。O(n\log n)

好家伙，经典 \gcd 容斥部分的求解成为了复杂度上限。

不过我们仍然有办法！

我们想，反正都是容斥，这里也用欧拉反演不好吗？

比如说：

tot1&=(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_{i},b_{j}))^2\\ &=(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{d|a_i,d|b_j}\varphi(d))^2\\ &=(\sum\limits_{d=1}^{n}ca_dcb_d\varphi(d))^2 \end{aligned}

要用 f 求的还有一个 res1：

\begin{aligned}res1=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(a_i,b_j)^2\end{aligned}

啊这怎么欧拉反演。。。

我们参考欧拉反演的式子 n=\sum\limits_{d|n}\varphi(d)，也设计一个反演！

具体的，我们找到一个函数 \phi(x)，使得 n^2=\sum\limits_{d|n}\phi(d) 就能像上面那样反演了。

这里根据狄利克雷卷积有 id_2=\phi\ \ast\ I ，\phi=id_2\ \ast\ \mu 所以 \phi 是个积性函数，我们只要知道 \phi(p^k) 是多少即可。

我们可以轻易发现 \phi(p^k)=(p^k)^2-\phi(1)-\phi(p)-...-\phi(p^{k-1})=(p^k)^2-(p^{k-1})^2

然后就没了。

复杂度 O(n\log\log n)。

全村最快的代码

在洛谷上跑 n=5000000 只用 500ms （不计读入）。

#include<bits/stdc++.h>
#define N 5000010
#define MOD 1000000007
using namespace std;
void read(int &x)
{
    x=0;
    char ch=getchar();
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9')
    {
        x=x*10+ch-'0';
        ch=getchar();
    }
}
int f[N],phi[N],v[N],p[N];
int n,m,invn,sum1,sum2,ans;
int a[N],b[N],mul[N];
int ca[N],cb[N],ga[N],gb[N];
int ksm(int x,int k)
{
    int res=1;
    while(k)
    {
        if(k&1)res=1ll*res*x%MOD;
        x=1ll*x*x%MOD;
        k>>=1; 
    }
    return res;
}
int gcd(int x,int y)
{
    if(x==0)return y;
    return gcd(y%x,x);
}
int main()
{
    read(n);
    invn=ksm(n,MOD-2);
    for(int i=1;i<=n;i++)
    {
        int x;
        read(x);
        ++a[x];
        ++ca[x];//这边要多一步，等价于ca[x]=a[x] 
        m=max(m,x);
    }
    for(int i=1;i<=n;i++)
    {
        int x;
        read(x);
        ++b[x];
        ++cb[x];
        m=max(m,x);
    }
    mul[0]=1;
    for(int i=1;i<=n;i++)
        mul[i]=1ll*mul[i-1]*i%MOD;

    phi[1]=1;f[1]=1;
    for(int i=2;i<=m;i++)
    {
        if(!v[i])p[++p[0]]=i,phi[i]=i-1,f[i]=1ll*i*i%MOD-1;//f[i]就是我们自己求的积性函数！ 
        for(int j=1;j<=p[0]&&p[j]*i<=m;j++)
        {
            v[i*p[j]]=1;
            if(i%p[j]==0)
            {
                phi[i*p[j]]=phi[i]*p[j];
                f[i*p[j]]=1ll*f[i]*p[j]%MOD*p[j]%MOD;
                break;
            }
            else phi[i*p[j]]=phi[i]*phi[p[j]],f[i*p[j]]=1ll*f[i]*f[p[j]]%MOD;
        }
    }

    for(int i=1;i<=p[0];i++)
        for(int j=m/p[i];j>=1;j--)
            ca[j]+=ca[j*p[i]],cb[j]+=cb[j*p[i]];//高维前缀和！ 

    for(int i=m;i>=1;i--)
    {
        int tmp=1ll*ca[i]*cb[i]%MOD;
        sum1=(sum1+1ll*tmp*f[i])%MOD;
        sum2=(sum2+1ll*tmp*phi[i])%MOD;
    }
    sum2=1ll*sum2*sum2%MOD;

    for(int i=1;i<=m;i++)//初始化 
    {
        gb[i]=1ll*phi[i]*cb[i]%MOD;
        ga[i]=1ll*phi[i]*ca[i]%MOD;
    }

    for(int i=1;i<=p[0];i++)
        for(int j=1;j<=m/p[i];j++)
            (ga[j*p[i]]+=ga[j])%=MOD,(gb[j*p[i]]+=gb[j])%=MOD;//高维前缀和！ 

    for(int i=1;i<=m;i++)
    {
        sum2=(sum2-1ll*gb[i]*gb[i]%MOD*a[i]%MOD+MOD)%MOD;
        sum2=(sum2-1ll*ga[i]*ga[i]%MOD*b[i]%MOD+MOD)%MOD;
    }

    sum2=(sum2+sum1)%MOD;

    sum1=1ll*sum1*mul[n-1]%MOD;
    sum2=1ll*sum2*mul[n-2]%MOD;
    sum2=(sum1+sum2)%MOD;
    ans=(sum1-1ll*invn*sum2%MOD+MOD)%MOD;
    ans=1ll*ans*invn%MOD;
    printf("%d",ans);
}