【CQOI2017】小Q的表格
稍加推导就会发现(f(a,b)=acdot bcdot h(gcd(a,b)))。
初始时(h(n)=1)。
询问前(k)行(k)列时我们就反演：
[egin{align}
displaystyle
ans&=sum_{g=1}h(g)cdot g^2sum_{a=1}^{lfloorfrac{k}{g}
floor} sum_{b=1}^{lfloorfrac{k}{g}
floor}acdot bsum_{d|a,d|b}mu(d)\
&=sum_{g=1}h(g)cdot g^2sum_{d=1}^{lfloorfrac{k}{g}
floor}d^2cdot mu(d)cdot sum^2(lfloorfrac{k}{gd}
floor)
end{align}
]
其中(displaystyle sum(n)=sum_{i=1}^ni)。
我们设
[displaystyle
q(n)=sum_{d=1}^n d^2cdot mu(d)cdot sum^2(lfloorfrac{n}{d}
floor)
]
我们可以在(O(NlogN))的时间内求出所有的(q(i))。
考虑(q(n))与(q(n-1))的区别：只有在(d|n)的时候才会产生。
所以：
[displaystyle
q(n)=q(n-1)+sum_{d|n}d^2mu(d)(sum^2(frac{n}{d})-sum^2(frac{n}{d}-1))
]
于是我们就可以数论分块处理询问。
但是如果我们用树状数组处理修改，那么我们的复杂度就多了一个(log)，无法接受。所以我们牺牲修改时间来平衡询问时间。
使用分块来维护(h(g)cdot g^2)的前缀和。
代码：
#include<bits/stdc++.h>
#define ll long long
#define N 4000005
using namespace std;
inline ll Get() {ll x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}
int n,m;
ll a,b;
ll x,k;
const ll mod=1e9+7;
int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
bool vis[N];
int p[N],u[N];
ll g[N],f[N],sum[N];
ll inv[N],bel[N];
const int blk=2e3+7;
int add[blk];
void pre(int n) {
g[1]=u[1]=1;
for(int i=2;i<=n;i++) {
if(!vis[i]) {
g[i]=mod+1-inv[i];
p[++p[0]]=i;
}
for(int j=1;j<=p[0]&&1ll*i*p[j]<=n;j++) {
vis[i*p[j]]=1;
if(i%p[j]==0) {
g[i*p[j]]=g[i];
break;
}
g[i*p[j]]=1ll*(mod+1-inv[p[j]])*g[i]%mod;
}
}
for(int i=1;i<=n;i++) {
g[i]=(1ll*g[i]*i%mod*i%mod*i%mod+g[i-1])%mod;
}
}
void Add(int v,int x) {
int b=bel[v];
for(int i=(b-1)*blk+1;bel[i]==b;i++) (f[i]+=add[b])%=mod;
add[b]=0;
for(int i=v;bel[i]==b;i++) (f[i]+=x)%=mod;
for(int i=b+1;i<=bel[n];i++) (add[i]+=x)%=mod;
}
ll query(int v) {return (f[v]+add[bel[v]])%mod;}
int solve(int n) {
ll ans=0,last=0;
for(int i=1;i<=n;i=last+1) {
last=n/(n/i);
(ans+=1ll*(query(last)-query(i-1)+mod)*g[n/i])%=mod;
}
return ans;
}
int main() {
m=Get(),n=Get();
inv[1]=inv[0]=1;
for(int i=2;i<=n;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
pre(n);
for(int i=1;i<=n;i++) bel[i]=(i-1)/blk+1;
for(int i=1;i<=n;i++) f[i]=1ll*i*i%mod;
for(int i=1;i<=n;i++) f[i]=(f[i]+f[i-1])%mod;
while(m--) {
a=Get(),b=Get(),x=Get(),k=Get();
int g=gcd(a,b);
ll now=x%mod*inv[a]%mod*inv[b]%mod*g%mod*g%mod;
now=(now-(query(g)-query(g-1))+mod)%mod;
Add(g,now);
int ans=0;
cout<<solve(k)<<"
";
}
return 0;
}