LOJ#3020. 「CQOI2017」小 Q 的表格
阿新 • • 發佈：2022-02-19
...
首先要清楚這道題的修改是啥意思。
題目裡給的限制都特別像輾轉相減，於是可以往這個方向靠。
不難發現，對於任意的 \(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)\)
於是有 \(\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;
}
}