[BZOJ4815][CQOI2017]小Q的表格数论+分块

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=4815

题目中所给条件中的$(a,a+b)$和$(a,b)$的关系很瞩目。

然后大家都知道$(a,b)=(a,a-b)=(a,a+b)$，于是观察（猜）一下这个表格与gcd的关系。

可以发现每次修改$(a,b)$会影响到所有$(i,j)=(a,b)$的点，并且关系为$$f(i,j)=\frac{i}{a}*\frac{j}{b}*f(a,b)$$

所以只需要知道$f(d,d)$的值记为$f(d)$，就能推出其他的值。

然后慢慢推推推大概可以推到这一步$$ans=\sum_{d=1}^nf(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}(i,j)[(i,j)==1]$$

可以发现这个式子中$i$和$j$是对称的$$S(\frac{n}{d})=\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}(i,j)[(i,j)==1]$$

不妨先设$i>j$，于是我们有$$S′(n)=\sum_{i=1}^n\frac{φ(i)*i^{2}}{2}$$

由于$i$与$j$对称，所以有$$S(n)=2*S′(n)=\sum_{i=1}^nφ(i)*i^{2}$$

所以最终的答案就变成了$$ans=\sum_{d=1}^nf(d)S(\frac{n}{d})$$

我们记录$f$的前缀和，并且分块维护这个数列，而$S$很明显是可以预处理出来的。

询问了$m$次，于是总体复杂度应该是$O(m\sqrt{n})$

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
const int mod=1e9+;
int inline readint(){
int Num;char ch;
while((ch=getchar())<''||ch>'');Num=ch-'';
while((ch=getchar())>=''&&ch<='') Num=Num*+ch-'';
return Num;
}
ll inline readll(){
ll Num;char ch;
while((ch=getchar())<''||ch>'');Num=ch-'';
while((ch=getchar())>=''&&ch<='') Num=Num*+ch-'';
return Num;
}
void outint(int x){
if(x>=) outint(x/);
putchar(x%+'');
}
int inline gcd(int x,int y){
return !y?x:gcd(y,x%y);
}
int n,m;
int phi[],p[],cnt=;
int la,blk,add[];
int f[];
bool vis[];
void sieve(int n){
for(int i=;i<=n;i++){
if(!vis[i]){
p[++cnt]=i;
phi[i]=i-;
}
for(int j=;p[j]*i<=n;j++){
vis[p[j]*i]=true;
if(i%p[j]==){
phi[p[j]*i]=phi[i]*p[j];
break;
}
phi[p[j]*i]=phi[i]*(p[j]-);
}
phi[i]=(1LL*i*i%mod*phi[i]+phi[i-])%mod;
f[i]=(1LL*i*i+f[i-])%mod;
}
}
void modify(int x,int ad){
int l=(x-)/blk+,
r=min(n,l*blk);
for(int i=l+;i<=la;i++) add[i]=(add[i]+ad)%mod;
for(int i=x;i<=r;i++) f[i]=(f[i]+ad)%mod;
}
int inline qry(int x){
return x?(f[x]+add[(x-)/blk+])%mod:;
}
int main(){
m=readint();
n=readint();
f[]=phi[]=;
blk=(int)sqrt(n);
la=(n-)/blk+;
sieve(n);
for(int i=;i<=m;i++){
int a=readint(),
b=readint(),
g=gcd(a,b),
ans=;
ll x=readll();
int k=readint();
x=x/(1LL*(a/g)*(b/g))%mod;
modify(g,((x-qry(g)+qry(g-))%mod+mod)%mod);
for(int j=,now;j<=k;j=now+){
now=k/(k/j);
ans=(ans+1LL*(qry(now)-qry(j-)+mod)%mod*phi[k/j])%mod;
}
outint(ans);
putchar('\n');
}
return ;
}