NKOJ 4040 (CQOI 2017) 小Q的表格（莫比乌斯反演+分块+递推+线性筛/欧拉函数+分块+线性筛）

P4040小Q 的表格

问题描述

---

题目给出了一个有规律的表格，因此我们先随便修改一个数找一下所有被修改的数之间有没有什么规律，很容易发现好像被修改的数的行号和列号的gcd是一样的，于是我们考虑证明，实际上我们的修改过程和辗转相减的过程是一样的，因此很容易得证。
接着我们来考虑gcd一样的这些格子的数有什么特点，容易发现他们的倍数关系是固定的，等于行号列号乘积之商，所以我们用A[d]表示（d，d）这个格子的值，列出求和式

Ans=∑d=1kA[d]∑i=1k∑j=1k(gcd(i,j)==d)i∗jd2

接下来就有两种做法了。

---

做法一：

上面的式子看起来很眼熟，于是我们将他变形

Ans=∑d=1kA[d]∑i=1⌊kd⌋∑j=1⌊kd⌋(gcd(i,j)==1)i∗j

于是莫比乌斯反演一下，

f(x)=∑i=1⌊kd⌋∑j=1⌊kd⌋(gcd(i,j)==x)i∗j

F(x)=∑i=1⌊kd⌋∑j=1⌊kd⌋(x|gcd(i,j))i∗j

Ans=∑d=1kA[d]∗f(1)=∑d=1kA[d]∑x=1⌊kd⌋u(x)∗x2∗[⌊kdx⌋∗(⌊kdx⌋+1)2]2

现在，我们希望把A[d]后面的一坨解决掉，于是，我们令

g(i)=∑x=1iu(x)∗x2∗[⌊ix⌋∗(⌊ix⌋+1)2]2

Ans=∑d=1kA[d]∗g(⌊kd⌋)

现在我们发现，如果能预处理出g(i)的值，那么可以利用分块的方法在O(n√)的复杂度内处理每一次询问，问题在于如何预处理，直接算显然不靠谱，我们考虑递推。可以发现，如果⌊ix⌋和⌊i−1x⌋是一样的，那么这一项就被减掉了，简单推到一下可以得到

g(i)−g(i−1)=∑x|iu(x)∗i3x

不妨令右边这一坨为h(i)
因此我们有

g(i)−g(i−1)=h(i)

然后h(i)显然是积性的，可以用线性筛搞出来，因此我们可以在O(n)的复杂度内搞出我们想要的g(i)，之后就可以愉快的分块处理了。

最后还需要注意，分块时我们需要A[i]的前缀和，然而用树状数组维护查询是O(logn)的，肯定要T，
因此我们需要用到分块前缀和这种黑科技，它可以实现O(n√)插入，O(1)查询，方法就是将前缀和数组分块，一般分成n√块，然后我们每次更新前缀和相当于将一个区间整体增加一个d，因此我们借鉴线段树的lazy思想，用一个数组记录下整块的更新，然后暴力更新零碎的部分，时间复杂度不超过O(n√)，查询的时候只需要返回两个数组的和即可。

附上代码

#include<stdio.h>#include<iostream>#include<algorithm>#include<cstring>#include<cmath>#define ll long long#define N 4444444using namespace std;const ll mod=1000000007LL;ll A[N],S[N],Q[N],sn,ssn,n,m;ll pri[N],G[N],F[N],tot;bool mark[N];ll gcd(ll a,ll b){return b==0?a:gcd(b,a%b);}void Pretreat(){    ll i,j;    for(i=1;i<=n;i++)    {        A[i]=i*i%mod;        S[i]=(S[i-1]+A[i])%mod;    }    F[1]=1;G[1]=1;    for(i=2;i<=n;i++)    {        if(!mark[i])pri[++tot]=i,F[i]=i*i%mod*(i-1)%mod;        for(j=1;j<=tot&&pri[j]*i<=n;j++)        if(i%pri[j])mark[i*pri[j]]=1,F[i*pri[j]]=F[i]*F[pri[j]]%mod;        else {mark[i*pri[j]]=1;F[i*pri[j]]=F[i]*(pri[j]*pri[j]%mod)%mod*pri[j]%mod;}        G[i]=(G[i-1]+F[i])%mod;    }}void MD(ll x,ll d){    ll i,l=(x-1)/sn+1,r=l*sn;    for(i=x;i<=r&&i<=n;i++)S[i]=(S[i]+d)%mod;    for(i=l+1;i<=ssn;i++)Q[i]=(Q[i]+d)%mod;}ll GS(ll x){    ll l=(x-1)/sn+1;    if(x==0)return 0;    return (S[x]+Q[l])%mod;}int main(){    ll i,j,k,t,a,b,x,y,ans;    scanf("%lld%lld",&m,&n);    Pretreat();sn=sqrt(n);ssn=(n-1)/sn+1;    for(i=1;i<=m;i++)    {        scanf("%lld%lld%lld%lld",&a,&b,&x,&k);        y=gcd(a,b);        x=x/(a*b/y/y)%mod;        MD(y,x-A[y]);        A[y]=x;ans=0;        for(j=1;j<=k;j=t+1)        {            t=k/(k/j);            ans=(ans+G[k/j]*(GS(t)-GS(j-1))%mod)%mod;        }        printf("%lld\n",(ans+mod)%mod);    }}

---

做法二：

做法二非常的厉害，我们先回顾一下做法一中

Ans=∑d=1kA[d]∑i=1⌊kd⌋∑j=1⌊kd⌋(gcd(i,j)==1)i∗j

然后，厉害的来了，根据欧拉函数的性质，我们有

∑i=1ni∗[gcd(i,n)==1]=n∗φ(n)2

注意到i和j是相互独立的，于是套用公式，注意到上面公式中i是小于n的，但是我们要求的式子中，i和j并没有大小关系，于是需要乘上2，于是我们可以得到

Ans=∑d=1kA(d)∑i=1⌊kd⌋i2φ(i)

于是，后面一坨显然是积性的，之后同样分块处理即可。
关于这个公式的证明，注意到如果a和n互质，那么n-a和n也一定互质，因为假设p同时整除n和n-a，那么他一定整除a，矛盾。

以下代码引用自F.F.Chopin，同时不得不说，做法二比做法一跑得要快，因为预处理取模较少，更快。

#include<stdio.h>#include<cmath>#define LL long longconst int MAXN=4000005;const LL mod=1000000007;using namespace std;LL m,n,prime[MAXN>>3],phi[MAXN],tot,g[MAXN],w[MAXN],f[MAXN],tg[MAXN];bool ntpr[MAXN];LL a,b,x,k,d,nn;LL gcd(LL a,LL b){    if(!b)return a;    return gcd(b,a%b);}void getphi(int n){    phi[1]=g[1]=f[1]=w[1]=1;    for(LL i=2; i<=n; i++)    {        if(!ntpr[i])            prime[++tot]=i,phi[i]=i-1;        for(int j=1; j<=tot&&i*prime[j]<=n; j++)        {            ntpr[i*prime[j]]=1;            if(!(i%prime[j]))            {                phi[i*prime[j]]=phi[i]*prime[j];                break;            }            phi[i*prime[j]]=phi[i]*(prime[j]-1);        }        f[i]=i*i%mod;        g[i]=(g[i-1]+f[i]*phi[i]%mod)%mod;        w[i]=(w[i-1]+f[i])%mod;    }}void calc(LL k){    LL ans=0,i,ls;    for(i=1; i<=k; i=ls+1)    {        ls=k/(k/i);        ans=((w[ls]-w[i-1]+tg[(ls+nn-1)/nn]-tg[(i+nn-2)/nn])*g[k/i]%mod+ans)%mod;    }    printf("%lld\n",(ans+mod)%mod);}int main(){    scanf("%lld%lld",&m,&n);    getphi(n);    nn=sqrt(n);    LL c,v;    while(m--)    {        scanf("%lld%lld%lld%lld",&a,&b,&x,&k);        d=gcd(a,b);        c=f[d],f[d]=x/(a/d)/(b/d)%mod;        for(int i=d; i<=(nn-1+d)/nn*nn; i++)            w[i]=(w[i]+f[d]-c)%mod;        for(int i=(d+nn-1)/nn+1; i<=(n+nn-1)/nn; i++)            tg[i]=(tg[i]+f[d]-c)%mod;        calc(k);    }}