【BZOJ4815】[CQOI2017]小Q的表格（莫比烏斯反演，分塊）

題面

BZOJ
洛谷

題解

神仙題啊。
首先\(f(a,b)=f(b,a)\)告訴我們矩陣只要算一半就好了。
接下來是\(b*f(a,a+b)=(a+b)*f(a,b)\)
這個式子怎麼看呢？
\[\begin{aligned}b*f(a,a+b)&=(a+b)*f(a,b)\\\frac{f(a,a+b)}{a+b}&=\frac{f(a,b)}{b}\\\frac{f(a,a+b)}{a*(a+b)}&=\frac{f(a,b)}{a*b}\end{aligned}\]
這個式子很明顯類似於輾轉相減的式子，如果我們令\(c=(a+b)\)

，那麼等式就可以寫成
\[\frac{f(a,c)}{a*c}=\frac{f(a,c-a)}{a*(c-a)}\]
那麼進一步，意味著我們可以寫成輾轉相除的式子。
\[\frac{f(a,b)}{a*b}=\frac{f(a,a\%b)}{a*(a\%b)}\]
所以，類似於求\(gcd\)的終止狀態，我們可以利用這個寫出一個等式：
\[\frac{f(a,b)}{a*b}=\frac{f(gcd(a,b),gcd(a,b)}{gcd^2(a,b)}\]
為了方便後面直接令\(d=gcd(a,b)\)，即等式為
\[f(a,b)=\frac{a*b}{d^2}f(d,d)\]
那麼，不難發現，單次的修改只會對於\(d\)相等的位置產生影響。
考慮\(ans\)是個啥玩意
\[\begin{aligned}ans&=\sum_{i=1}^k\sum_{j=1}^kf(i,j)\\&=\sum_{d=1}^kf(d,d)\sum_{i=1}^k\sum_{j=1}^k\frac{i*j}{d^2}[gcd(i,j)=d]\\&=\sum_{d=1}^kf(d,d)\sum_{i=1}^{k/d}\sum_{j=1}^{k/d}ij[gcd(i,j)=1]\end{aligned}\]
後面一半的式子可以用莫比烏斯反演直接計算值。
當然了，也可以這樣子推:
\[\begin{aligned}\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=1]&=2*\sum_{i=1}^ni\sum_{j=1}^ij[gcd(i,j)=1]\\&=2*\sum_{i=1}^n\frac{i\varphi(i)}{2}\\&=\sum_{i=1}^ni^2\varphi(i)\end{aligned}\]
後面那一步化簡簡單證明一下，假設\(x\)與\(n\)互質，那麼\(n-x\)也與\(n\)互質，因此所有與\(n\)互質的數的和是兩兩配對的。
也就是後面這一部分是可以提前預處理出來的，因為是\(k/d\)，也就是等價於可以數論分塊，所以我們現在唯一要做的就是維護\(f(i,i)\)的字首和。
而修改操作顯然是把當前位置會影響的位置的值全部對應的擴倍。也就是會影響到\(f(gcd,gcd)\)這個位置的值。我們要找個方法快速計算\(f(i,i)\)的字首和。
考慮資料範圍，操作次數很少，但是操作的範圍很大。單次詢問的時候我們有一個\(\sqrt k\)的數論分塊的複雜度，也就是\(2*10^3\)，操作次數有\(10^4\)，所以我們顯然不能帶\(log\)的計算字首和。那麼考慮分塊，將計算字首和的複雜度將至\(O(1)\)，而修改複雜度變為\(O(\sqrt n)\)。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;
#define ll long long
#define MOD 1000000007
#define MAX 4000400
#define BLK 2200
inline ll read()
{
    ll x=0;bool t=false;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=true,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return t?-x:x;
}
int fpow(int a,int b)
{
    int s=1;
    while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}
    return s;
}
int phi[MAX],pri[MAX],tt,s[MAX];
bool zs[MAX];
void Pre(int n)
{
    phi[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!zs[i])pri[++tt]=i,phi[i]=i-1;
        for(int j=1;j<=tt&&i*pri[j]<=n;++j)
        {
            zs[i*pri[j]]=true;
            if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1);
            else{phi[i*pri[j]]=phi[i]*pri[j];break;}
        }
    }
    for(int i=1;i<=n;++i)s[i]=1ll*i*i%MOD*phi[i]%MOD;
    for(int i=1;i<=n;++i)s[i]=(s[i]+s[i-1])%MOD;
}
int pre[BLK],ps[BLK],bs[BLK][BLK],val[MAX];
int b[MAX],tot,blk;
int n,m;
int Query(int x){if(!x)return 0;return (pre[b[x]-1]+bs[b[x]][x-blk*(b[x]-1)])%MOD;}
void Modify(int x,int w)
{
    val[x]=w;int p=x-blk*(b[x]-1);
    for(int i=x;b[i]==b[x];++i,++p)
        bs[b[x]][p]=(bs[b[x]][p-1]+val[i])%MOD;
    ps[b[x]]=bs[b[x]][p-1];
    for(int i=1;i<=tot;++i)pre[i]=(pre[i-1]+ps[i])%MOD;
}
int Calc(int k)
{
    int ret=0;
    for(int i=1,j;i<=k;i=j+1)
    {
        j=k/(k/i);
        ret=(ret+1ll*(Query(j)-Query(i-1)+MOD)*s[k/i])%MOD;
    }
    return ret;
}
int main()
{
    m=read();n=read();Pre(n);blk=sqrt(n);
    for(int i=1;i<=n;++i)b[i]=(i-1)/blk+1;
    tot=b[n];
    for(int i=1;i<=n;++i)val[i]=1ll*i*i%MOD;
    for(int i=1;i<=tot;++i)
        for(int j=1,a=(i-1)*blk+1;j<=blk&&a<=n;++j,++a)
            ps[i]=(ps[i]+val[a])%MOD,bs[i][j]=(bs[i][j-1]+val[a])%MOD;
    for(int i=1;i<=tot;++i)pre[i]=(pre[i-1]+ps[i])%MOD;
    while(m--)
    {
        int a=read(),b=read(),x=read()%MOD,k=read();
        int d=__gcd(a,b);
        Modify(d,1ll*x*d%MOD*d%MOD*fpow(1ll*a*b%MOD,MOD-2)%MOD);
        printf("%d\n",Calc(k));
    }
    return 0;
}