BZOJ4815: [Cqoi2017]小Q的表格

BZOJ4815

b×f(a,a+b)=(a+b)∗f(a,b) 很像辗转相减法。。 那么每次修改点 (a,b) 的值，会修改所有满足 gcd(i,j)==gcd(a,b) 的点 (i,j) 的值。 记 d=gcd(a,b) ，那么 fi,j=x∗i/d∗j/d (gcd(i,j)==d

那么可以转化为，每次修改对角线上的值。记 Numi 为 fi,i

ans=∑d=1nNumd∑i=1n∑j=1n[gcd(i,j)==d]i/d∗j/d ans=∑d=1nNumd∑I=1n/d∑j=1n/d[gcd(i,j)==1]i∗j 令 G(n)=∑nI=1∑nj=1[gcd(i,j)==1]i∗j 根据 BZOJ2154: Crash的数字表格 的做法 G(n)=∑i=1nμ(i)∗i2∗S2(n/i) S(i)=i(i+1)2 nn√ 求 G(n) 无法接受，复杂度还是高了一些。 考虑递推。当且仅当 n|d 时， ⌊nd⌋ 比 ⌊n−1d⌋ 大1 那么 G(n)=G(n−1)+∑i|nμ(i)i2(S2(⌊ni⌋)−S2(⌊n−1i⌋) 令 H(n)=∑i|nμ(i)i2(S2(⌊ni⌋)−S2(⌊n−1i⌋)=n3∑i|nμ(i)1i 令 h(n)=∑i|nμ(i)1i 发现 h(n) 是一个积性函数。 证明一下： 对于两个互质的数 x,y 记 x 的因子集为A，大小为 s1 记 y 的因子集为B，大小为 s2 那么 x∗y 的因子集 C 为A∗B，大小为 s1∗s2 举个例子 x=10,y=21 A=1,2,5,10  s1=4 B=1,3,7,21  s2=4 那么 C=⎛⎝⎜⎜⎜1372126144251535105103070210⎞⎠⎟⎟⎟

显然 μ(a)1a∗μ(b)1b=μ(a∗b)1a∗b 那么相当于 h(xy) 为 A 中每个μ(a)∗1a乘以 B 中每个μ(b)∗1b 然后就可以线性筛出 h(n)

ans=∑d=1nNumdG(n/d)

然后就可以下底函数分块来求了。但是有修改操作,如果用树状数组，修改复杂度为 O(logn) ，但是一次查询的复杂度为 O(lognn√) ，会 T 掉。 那么可以考虑牺牲一下修改的复杂度，用分块来维护前缀和，修改变成O(n√)，查询变成 O(1) 了。

然后这个题终！于！做！完！了！ 19s 算是给卡过去了。。好像 S2(n)=∑ni=1φ(i)∗i2 但是不理解为什么啊 QAQ ，如有神犇愿意解答感激不尽！

【代码】

#include <cstdio> #include <iostream> #include <algorithm> #include <cstring> #include <cmath> #include <queue> #define N 4000005 #define M 600005 #define Mod 1000000007 #define INF 0x7fffffff using namespace std; typedef long long ll; typedef unsigned long long ull; const ull base=31; ll read() { ll x=0,f=1;char ch=getchar(); while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();} while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();} return x*f; } int m,n,block,Mx; int Prime[N],bl[N],R[N];bool Not_Prime[N]={1,1}; ll f[N],h[N],g[N],Inv[N],sum[2005][2005],Sum[2005]; ll SSqr(int x){ return 1LL*x*x%Mod*x%Mod; } ll gcd(ll a,ll b){ return b==0?a:gcd(b,a%b); } void Pre() { Inv[1]=h[1]=g[1]=1;block=sqrt(m); for(register int i=2;i<=m;i++) { Inv[i]=(Mod-Mod/i)*Inv[Mod%i]%Mod; if(!Not_Prime[i]) Prime[++Prime[0]]=i,h[i]=(1-Inv[i]+Mod)%Mod; for(register int j=1;j<=Prime[0]&&i*Prime[j]<=m;j++) { Not_Prime[i*Prime[j]]=1; if(i%Prime[j]==0) { h[i*Prime[j]]=h[i]; break; } h[i*Prime[j]]=h[i]*h[Prime[j]]%Mod; } g[i]=(g[i-1]+SSqr(i)*h[i]%Mod)%Mod; } for(register int i=1;i<=m;i++) { f[i]=1LL*i*i; bl[i]=(i-1)/block+1; Sum[bl[i]]=(Sum[bl[i]]+f[i])%Mod; sum[bl[i]][i-R[bl[i]-1]]=(sum[bl[i]][i-R[bl[i]-1]-1]+f[i])%Mod; if(i==n||(i%block==0)) R[bl[i]]=i; } Mx=bl[m]; for(int i=1;i<=Mx;i++) Sum[i]=(Sum[i]+Sum[i-1])%Mod; } void Update(ll d) { int Bl=bl[d]; for(int i=d;i<=R[Bl];i++) sum[Bl][i-R[Bl-1]]=(sum[Bl][i-1-R[Bl-1]]+f[i])%Mod; Sum[Bl]=(Sum[Bl-1]+sum[Bl][R[Bl]-R[Bl-1]])%Mod; for(int i=Bl+1;i<=Mx;i++) Sum[i]=Sum[i-1]+sum[i][R[i]-R[i-1]]; } ll Get_Sum(int x,int y) { int sx=bl[x],sy=bl[y]; if(sx==sy) return (sum[sx][y-R[sx-1]]-sum[sx][x-1-R[sx-1]]+Mod)%Mod; return ((sum[sx][R[sx]-R[sx-1]]-sum[sx][x-1-R[sx-1]]+sum[sy][y-R[sy-1]]+Sum[sy-1]-Sum[sx])%Mod+Mod)%Mod; } void Solve(int x) { int pos; ll ans=0; for(int i=1;i<=x;i=pos+1) { pos=(x/(x/i)); ans=(ans+Get_Sum(i,pos)*g[x/i]%Mod)%Mod; } printf("%lld\n",ans); } int main() { n=read(),m=read(); Pre(); for(int i=1;i<=n;i++) { static int a,b,k;ll x; a=read(),b=read(),x=read(),k=read();x%=Mod; int GCD=gcd(a,b); f[GCD]=x*Inv[a/GCD]%Mod*Inv[b/GCD]%Mod; Update(GCD); Solve(k); } return 0; }