【BZOJ4815】 [Cqoi2017]小Q的表格

题意其实很简单

f(a,b)=f(b,a)

b*f(a,a+b)=(a+b)*f(a,b)

然后我们发现各个格子之间的影响是相互的，然后对于任意f(a,b)（不妨设a>b，等于号最后说）总可以被f(a-b,b)影响（算算就好了），然后你发现这tm就是手动的辗转相除法，然后你就get到了两个格子相互之间可以影响的条件：gcd(a,b)相同！！

然后自然而然地想到了对于gcd(a,b)不同的分开处理，每次只要记录下修改后的f(a,a)的值，因为gcd(a,b)相同的和就是f(gcd(a,b),gcd(a,b))乘上一个系数

考虑这个系数，对于要处理前k行的条件下（gcd不是1的的话只要k除一下gcd得到新的k就好了）

系数为：

sum[k]=∑i=1k∑j=1k[gcd(i,j)=1]i⋅j=∑i=1k∑j=1ki⋅j∑d|gcd(i,j)μ(d)=∑d=1k((⌊nd⌋+1)⋅⌊nd⌋2)2

然后我们发现这个sum数组我们好像可能大概八成不会直接线性筛搞出来（反正我不会），所以然后我们玩个差分：
对于随便一个d，它对于k=id~(i+1)d-1的贡献是相同的，然后对于每个相同的贡献都可以视作一个区间加操作，然后差分一下，就变成了两个单点加操作，因为

o(∑i=1nni)=o(nlogn)

所以我们就学会了o(nlogn)的预处理方式！！然而n=4000000，亲测tle了。。。。。。
所以我们一定要找到一个o(n)的预处理姿势！
依然是将原数组差分，考虑得到的数组就变成了一个更漂亮的形式：

sum′[k]=∑d|kμ(d)⋅k3d

(怎么推得，其实前面相邻的那两个平方的差值就是某东西的立方值，反正还是挺好推的)
然后考虑更简单的问题

sum′′[k]=∑d|kμ(d)⋅kd

这显然是积性的（这看不出来我也没办法了）
所以考虑

sum′′[pq]=sum′′[pq−1]⋅p

（p为一个质数，from q234rty）
然后我们就学会了o(n)筛出sum[n]啦！
之后的修改已经转为了修改（前面应该说过了吧）f(gcd(a,b),gcd(a,b))
然后用某个东西维护f(gcd(a,b),gcd(a,b))的前缀和（用什么最后再说），之后询问k行k列的和就变成了

∑i=1ksum[⌊ki⌋]⋅f(i,i)

对于

⌊ki⌋

相同的分块，每次用某个东西（同前面）查询f(i,i)的前缀和就好了
最后考虑一下用什么东西。
我用树状数组开o2+特殊卡常技巧卡过去了（洛谷上），然后后来又听了q234rty的建议试了试分块，似乎一个点能快300ms+。
如果你卡在了tle的话，请看我的分块部分（依旧from q234rty）！

#pragma GCC optimize("O2")#include <bits/stdc++.h>#define N 4000009#define gc getchar()#define mod 1000000007#define ll long longusing namespace std;ll m,n,x,y,a,k,val[N],sum[N],bit[N],s[N],pf[N],lf[N];ll pri[N],pd[N],mu[N],cnt,Max[N];ll gcd(ll x,ll y){    if (y==0) return x;    else return gcd(y,x%y);}void add(ll x,ll y){    ll t=(x>>11)+1;    for (;x<(t<<11);x++) s[x]+=y;    for (;t<=(n>>11);t++) bit[t]+=y;}ll qry(ll x){    return (s[x]+bit[x>>11])%mod;}ll read(){    char ch;    ll x=1;    while (ch=gc,ch<'0'||ch>'9') if (ch=='-') x=-1;    ll s=ch-'0';    while (ch=gc,ch>='0'&&ch<='9') s=s*10+ch-'0';    return s*x;}int main(){    m=read(),n=read();    mu[1]=pd[1]=sum[1]=1;    for (ll i=1;i<=n;i++) pf[i]=i*i%mod,lf[i]=pf[i]*i%mod;    for (ll i=2;i<=n;i++)    {        if (!pd[i])        {            mu[i]=-1,pri[++cnt]=i;            sum[i]=i-1;            Max[i]=i;        }        for (ll j=1;j<=cnt&&pri[j]*i<=n;j++)        {            pd[i*pri[j]]=1;            if (i%pri[j]==0)            {                Max[pri[j]*i]=Max[i]*pri[j];                mu[i*pri[j]]=0;                if (Max[pri[j]*i]==pri[j]*i)                    sum[i*pri[j]]=(sum[i]*pri[j])%mod;                else                    sum[i*pri[j]]=sum[i/Max[i]]*sum[Max[pri[j]*i]]%mod;                break;            }            sum[i*pri[j]]=sum[i]*sum[pri[j]]%mod;            mu[i*pri[j]]=-mu[i];            Max[i*pri[j]]=pri[j];        }    }    for (ll i=1;i<=n;i++)        sum[i]=(sum[i]*pf[i]%mod+sum[i-1])%mod;    //for (ll i=1;i<=n;i++) cout<<sum[i]<<endl;    for (ll i=1;i<=n;i++)        s[i]=s[i-1]+(val[i]=pf[i]);    for (ll i=1;i<=m;i++)    {        x=read(),y=read(),a=read(),k=read();        ll Gcd=gcd(x,y),ans=0;        add(Gcd,(a/(x*y/Gcd/Gcd)-val[Gcd])%mod);        val[Gcd]=a/(x*y/Gcd/Gcd)%mod;        ll last=0;        for (ll j=1,l;j<=k;j=l+1)        {            l=k/(k/j);//[j,l]            ll t=qry(l);            ans=ans+sum[k/j]*(t-last)%mod;            last=t;        }        printf("%lld\n",(ans%mod+mod)%mod);    }    return 0;}