【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]=\sum_{i=1}^{k} \sum_{j=1}^{k}[gcd(i,j)=1]i·j \\=\sum_{i=1}^{k} \sum_{j=1}^{k}i·j\sum_{d|gcd(i,j)}\mu(d)\\=\sum_{d=1}^{k}(\frac{(\left \lfloor \frac{n}{d} \right \rfloor+1)·\left \lfloor \frac{n}{d} \right \rfloor}{2})^{2}$$
然后我们发现这个sum数组我们好像可能大概八成不会直接线性筛搞出来（反正我不会），所以然后我们玩个差分：
对于随便一个d，它对于k=id~(i+1)d-1的贡献是相同的，然后对于每个相同的贡献都可以视作一个区间加操作，然后差分一下，就变成了两个单点加操作，因为
$$o(\sum_{i=1}^{n}\frac{n}{i})=o(nlogn)$$
所以我们就学会了o(nlogn)的预处理方式！！然而n=4000000，亲测tle了。。。。。。
所以我们一定要找到一个o(n)的预处理姿势！
依然是将原数组差分，考虑得到的数组就变成了一个更漂亮的形式：
$$sum'[k]=\sum_{d|k}\frac{\mu(d)·k^{3}}{d}$$
(怎么推得，其实前面相邻的那两个平方的差值就是某东西的立方值，反正还是挺好推的)
然后考虑更简单的问题
$$sum''[k]=\sum_{d|k}\frac{\mu(d)·k}{d}$$
这显然是积性的（这看不出来我也没办法了）
所以考虑
$$sum''[p^{q}]=sum''[p^{q-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列的和就变成了
$$\sum_{i=1}^{k}sum[\left \lfloor \frac{k}{i} \right \rfloor]·f(i,i)$$
对于 $$\left \lfloor \frac{k}{i} \right \rfloor$$
相同的分块，每次用某个东西（同前面）查询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 long
using 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;
}