稍加推导就会发现\(f(a,b)=a\cdot b\cdot h(gcd(a,b))\)。
初始时\(h(n)=1\)。
询问前\(k\)行\(k\)列时我们就反演:
\[
\begin{align}
\displaystyle
ans&=\sum_{g=1}h(g)\cdot g^2\sum_{a=1}^{\lfloor\frac{k}{g}\rfloor} \sum_{b=1}^{\lfloor\frac{k}{g}\rfloor}a\cdot b\sum_{d|a,d|b}\mu(d)\&=\sum_{g=1}h(g)\cdot g^2\sum_{d=1}^{\lfloor\frac{k}{g}\rfloor}d^2\cdot \mu(d)\cdot sum^2(\lfloor\frac{k}{gd}\rfloor)
\end{align}
\]
其中\(\displaystyle sum(n)=\sum_{i=1}^ni\)。
我们设
\[
\displaystyle
q(n)=\sum_{d=1}^n d^2\cdot \mu(d)\cdot sum^2(\lfloor\frac{n}{d}\rfloor)
\]
我们可以在\(O(NlogN)\)的时间内求出所有的\(q(i)\)。
考虑\(q(n)\)与\(q(n-1)\)的区别:只有在\(d|n\)的时候才会产生。
所以:
\[
\displaystyle
q(n)=q(n-1)+\sum_{d|n}d^2\mu(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)<<"\n";
}
return 0;
}
原文:https://www.cnblogs.com/hchhch233/p/10507687.html