【CQOI2017】bzoj4815 小Q的表格

根据辗转相减可以知道，相互影响的就是$\gcd$相同的那些数。又根据条件$2$可以知道，所有$\gcd$相同的数的比例是不会改变的，又因为最开始$a(x,y)=xy$是一组解，因此每次操作都相当于把所有$\gcd$相同的数修改为$a(x,y)=kxy$。因此，对于修改我们只需要维护$f(d)$表示对于所有$\gcd(x,y)=d$的$(x,y)$有$a(x,y)=f(d)*xy$。
接下来考虑询问，也就是求

$$\begin{align} & \sum_{i=1}^n\sum_{j=1}^nf(\gcd(i,j))*ij \\ = & \sum_{d=1}^nf(d)\sum_{i=1}^n\sum_{j=1}^nij*[\gcd(i,j)=d] \\ = & \sum_{d=1}^nf(d)*d^2\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ij*e(\gcd(i,j)) \\ = & \sum_{d=1}^nf(d)*d^2\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ij\sum_{x|i,x|j} \mu(x) \\ = & \sum_{d=1}^nf(d)*d^2\sum_{x=1}^{\lfloor\frac nd\rfloor}\mu(x)x^2\sum_{i=1}^{\lfloor\frac n{dx}\rfloor} \sum_{j=1}^{\lfloor\frac n{dx}\rfloor}ij \end{align}$$
记
$$s(n)=\sum_{i=1}^ni=\frac{n(n+1)}2$$
$$g(n)=\sum_{i=1}^n\mu(i)*i^2*s^2(\lfloor\frac ni\rfloor)$$
如果能够预处理出 $g$，就可以下底函数分块来求解原问题。但是如果直接依次求是$O(n\sqrt n)$的，无法接受。考虑递推。注意到 $\lfloor\frac nd\rfloor$比 $\lfloor\frac{n-1}d\rfloor$大 $1$，当且仅当$d|n$，也就是 $n$除以$d$余数为零。因此可以得到
$$\begin{align} g(n) & =g(n-1)+\sum_{i|n}\mu(i)*i^2(s^2(\frac ni)-s^2(\frac ni-1)) \\ & = g(n-1)+n^3\sum_{i|n}\mu(i)*\frac 1i \end{align}$$
其中 $h(n)=\sum_{i|n}\mu(i)*\frac 1i$为积性函数可以线性筛，于是我们就可以 $O(n)$预处理出 $g$。
还有一个问题，在对答案进行下底函数分块的时候需要维护$f(i)*i^2$的前缀和，如果直接用树状数组的话，修改是 $O(\log n)$的，询问是 $O(\sqrt n\log n)$的，无法接受。因此可以分块，牺牲修改的复杂度来降低询问的复杂度，单次修改 $O(\sqrt n)$，单次询问 $O(1)$。这样最终的复杂度是 $O(n+m\sqrt n)$。
好像有一点卡常，考虑到比较小的 $d=\gcd(x,y)$被调用的次数比较多，可以在分块的时候维护后缀和，这样修改次数可以减少。

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define LL long long
const int p=1000000007,maxn=4000010,maxT=2010;
LL f[maxn],g[maxn],h[maxn],inv[maxn],
sum[maxT],sumin[maxT][maxT];
int prm[maxn],have[maxn],bel[maxn],L[maxT],R[maxT],
q,n,T,tot;
int gcd(int x,int y)
{
    return y?gcd(y,x%y):x;
}
void pause()
{
    int x;
    x=1;
}
void init()
{
    inv[1]=g[1]=h[1]=1;
    for (int i=2;i<=n;i++)
    {
        inv[i]=(p-inv[p%i])*(p/i)%p;
        if (!have[i])
        {
            prm[++tot]=i;
            h[i]=(1-inv[i]+p)%p;
        }
        for (int j=1;j<=tot&&(LL)i*prm[j]<=n;j++)
        {
            have[i*prm[j]]=1;
            if (i%prm[j]==0)
            {
                h[i*prm[j]]=h[i];
                break;
            }
            else h[i*prm[j]]=h[i]*h[prm[j]]%p;
        }
        g[i]=(g[i-1]+(LL)i*i%p*i%p*h[i])%p;
    }
    T=sqrt(n);
    for (int i=1;i<=T;i++)
    {
        L[i]=R[i-1]+1;
        R[i]=i==T?n:L[i]+T-1;
        for (int j=L[i];j<=R[i];j++)
        {
            f[j]=1;
            bel[j]=i;
            sumin[i][j-L[i]+1]=(sumin[i][j-L[i]]+(LL)j*j)%p;
        }
    }
    for (int i=T;i;i--)
        sum[i]=(sum[i+1]+sumin[i][R[i]-L[i]+1])%p;
}
LL getsum(int l,int r)
{
    //if (l==10&&r==13) pause();
    //LL ret1;
    int x=bel[l],y=bel[r];
    if (x==y) return (sumin[x][r-L[x]+1]-sumin[x][l-L[x]]+p)%p;
    return (sum[x+1]-sum[y]+p+sumin[y][r-L[y]+1]+sumin[x][R[x]-L[x]+1]-sumin[x][l-L[x]]+p)%p;
    /*LL ret=0;
    for (int i=l;i<=r;i++) ret=(ret+f[i]*i%p*i%p)%p;
    if (ret!=ret1) pause();
    return ret;*/
}
void update(int d)
{
    for (int j=d;j<=R[bel[d]];j++)
        sumin[bel[d]][j-L[bel[d]]+1]=(sumin[bel[d]][j-L[bel[d]]]+f[j]*j%p*j)%p;
    for (int j=bel[d];j;j--)
        sum[j]=(sum[j+1]+sumin[j][R[j]-L[j]+1])%p;
}
void solve(int m)
{
    LL ans=0;
    int u;
    for (int j=1;j<=m;j=u+1)
    {
        u=m/(m/j);
        ans=(ans+getsum(j,u)*g[m/j])%p;
    }
    printf("%lld\n",ans);
}
int main()
{
    /*freopen("table.in","r",stdin);
    freopen("table.out","w",stdout);*/
    int xx,yy,m,d;
    LL x;
    scanf("%d%d",&q,&n);
    init();
    while (q--)
    {
        scanf("%d%d%lld%d",&xx,&yy,&x,&m);
        x%=p;
        d=gcd(xx,yy);
        f[d]=(LL)x*inv[xx]%p*inv[yy]%p;
        update(d);
        solve(m);
    }
}