LOJ#3020. 「CQOI2017」小 Q 的表格
首先要清楚这道题的修改是啥意思。
题目里给的限制都特别像辗转相减，于是可以往这个方向靠。
不难发现，对于任意的 \(f(a,b),f(c,d)\),如果有 \(\gcd(a,b)=\gcd(c,d)\) ，那么 \(\frac{f(a,b)}{ab}=\frac{f(c,d)}{cd}\)，也就是说 \(f(a,b)=x -> f(c,d)=cd\frac{x}{ab}\)。
可以设一个函数 \(g(x)\) 表示所有的满足 \(\gcd(a,b)=g\) 的 \(f(a,b)=abg(x)\)。
此时可以列出式子：
\[\begin{aligned}
ans&=\sum^{k}_{d=1}d^2g(d)\sum^{\lfloor\frac{k}{d}\rfloor}_{i=1}\sum^{\lfloor\frac{k}{d}\rfloor}_{j=1}[\gcd(i,j)=1]ij\\
&=\sum^{k}_{T=1}\binom{\lfloor\frac{k}{T}\rfloor}{2}^2\sum_{d|T}d^2g(d)\mu(\frac{T}{d})\\
\end{aligned}
\]
如果你得到了上面这个式子并尝试化简，那么恭喜你，你掉进坑列了
或许有一些大佬真的可以优化到可以接受的复杂度，但是我的水平最多只能优化到 \(O(mk^{\frac{3}{4}})\) ，是过不了的。
实际上，这个式子的正确方向应该是这样的：
\[\begin{aligned}
ans&=\sum^{k}_{d=1}d^2g(d)\sum^{\lfloor\frac{k}{d}\rfloor}_{i=1}\sum^{\lfloor\frac{k}{d}\rfloor}_{j=1}[\gcd(i,j)=1]ij\\
&=\sum^{k}_{d=1}g(d)(-1+2\sum^{\lfloor\frac{k}{d}\rfloor}_{i=1}i\sum^{i}_{j=1}[\gcd(i,j)=1]j)\\
\end{aligned}
\]
因为 \((a,b) = (a,a-b)\) ,所以所有和 \(i\) 互质的数可以分成 \(\frac{\varphi(i)}{2}\) 组,每组的和为 \(a\) 。
于是有 \(\sum^{i}_{j=1}[\gcd(i,j)=1]=\frac{i\varphi(i)}{2}\) ,同时 \(i=1\) 的时候要特判。
回到前面那个式子,因为外面有个 \(\times 2\),并且把这个 \(2\) 抵消之后刚好把 \(1\) 的特判给去掉了,可以进一步优化到：
\[\begin{aligned}
ans&=\sum^{k}_{d=1}d^2g(d)\sum^{\lfloor\frac{k}{d}\rfloor}_{i=1}i^2\varphi(i)\\
\end{aligned}
\]
后面那东西可以提前预处理，直接用数论分块做，同时用分块维护单点修改区间查询，可以做到 \(O(m\sqrt{n})\) 的复杂度。
code
#include<cstdio>
#include<iostream>
#include<vector>
#include<algorithm>
#include<cstring>
#include<queue>
#include<cassert>
using namespace std;
#define ll long long
#define ri int
#define pii pair<int,int>
const ll mod=1e9+7;
ll add(ll x,ll y){return (x+=y)<mod?x:x-mod;}
ll dec(ll x,ll y){return (x-=y)<0?x+mod:x;}
ll ksm(ll d,ll t,ll res=1){for(;t;t>>=1,d=d*d%mod) if(t&1) res=res*d%mod;return res;}
const int B=2048;
const int MAXN=4e6+10*B+1;
int n,m,phi[MAXN],prim[MAXN],sum[MAXN],flag[MAXN],f[MAXN];
void sieve(){
phi[1]=1;
for(ri i=2;i<=n;++i){
if(!flag[i]) prim[++prim[0]]=i,phi[i]=i-1;
for(ri j=1;j<=prim[0]&&i*prim[j]<=n;++j){
flag[i*prim[j]]=1;
if(i%prim[j]==0){
phi[i*prim[j]]=phi[i]*prim[j];
break;
}
phi[i*prim[j]]=phi[i]*(prim[j]-1);
}
}
}
struct Block{
int rsum[MAXN/B+10],bsum[MAXN];
void insert(int x,int w){
int l=x/B*B,r=l+B-1;
for(ri i=x;i<=r;++i) bsum[i]=add(bsum[i],w);
for(ri i=x/B;i<=n/B;++i) rsum[i]=add(rsum[i],w);
}
int ask(int x){
// int tot=0;
// for(ri i=1;i<=x;++i) tot=add(tot,1ll*f[i]*i*i%mod);
if(x/B){
// assert(tot==add(bsum[x],rsum[x/B-1]));
return add(bsum[x],rsum[x/B-1]);
}
else{
// assert(tot==bsum[x]);
return bsum[x];
}
}
int query(int l,int r){return dec(ask(r),ask(l-1));}
void init(){
for(ri i=1;i<=n;++i) bsum[i]=1ll*i*i%mod,rsum[i/B]=add(rsum[i/B],1ll*i*i%mod);
for(ri i=0;i<=n/B;++i){
int l=i*B,r=l+B-1;
for(ri j=l+1;j<=r;++j) bsum[j]=add(bsum[j-1],bsum[j]);
if(i) rsum[i]=add(rsum[i],rsum[i-1]);
}
}
}T;
int gcd(int x,int y){for(;y;x%=y,swap(x,y));return x;}
int main(){
// freopen("rand.in","r",stdin);
// freopen("1.out","w",stdout);
ios::sync_with_stdio(false);
cin>>m>>n;sieve();
for(ri i=1;i<=n;++i) sum[i]=add(sum[i-1],1ll*i*i%mod*phi[i]%mod),f[i]=1;T.init();
// for(ri i=1;i<=n;++i) printf("%d ",phi[i]);puts("");
while(m--){
int x,y,k,g;ll w,d;
cin>>x>>y>>w>>k;
g=gcd(x,y);
d=dec(w%mod*ksm(1ll*x*y%mod,mod-2)%mod,f[g]);
f[g]=add(f[g],d);
T.insert(g,d*g%mod*g%mod);
int ans=0;
for(ri l=1,r;l<=k;l=r+1){
r=k/(k/l);
// printf("[%d %d %d]\n",l,r,T.query(l,r));
ans=add(ans,1ll*T.query(l,r)*sum[k/l]%mod);
}
cout<<ans<<endl;
}
}
