为什么我连分块都想不到啊...
题意
定义一个矩阵$f$满足
$ f(a,b)=f(b,a)$
$ b·f(a,a+b)=(a+b)·f(a,b)$
初始$ f(a,b)=ab$
有$ m$次修改，每次会将$ f(x,y)$改成$ c$，并修改这个矩阵使得仍然满足以上两条
保证修改后矩阵中每个数均为正整数
你需要在每次修改后求出$ \sum\limits_{i=1}^k\sum\limits_{j=1}^kf(i,j)$
数据范围满足$ m \leq 10000$，涉及到的矩形范围不超过$4·10^6$
$Solution $
给定的信息类似一个辗转相减的过程
和求$ gcd(x,y)$的过程非常像
我们不断辗转相减可以得到$ f(x,y)=f(d,d)*\frac{xy}{d^2}$
其中$ d=gcd(x,y)$
考虑求答案的式子的意义
$ ans=\sum\limits_{i=1}^k\sum\limits_{j=1}^kf(i,j)$
$ans=\sum\limits_{k=1}^nf(k,k)\sum\limits_{i=1}^{\frac{n}{k}}\sum\limits_{j=1}^{\frac{n}{k}}ij[gcd(i,j)=1]$
设
$ A(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^nij[gcd(i,j)=1]$
有
$A(n)=\sum\limits_{i=1}^n i^2\varphi(i)$
原理是若$ x<y,x和y互质则y-x和y互质$
这样可以化简为
$ans=\sum\limits_{k=1}^nf(k,k)A(\frac{n}{k})$
显然$ A(x)$可以预处理，然后数论分块计算
单次$ O(\sqrt{4000000})=O(2000)$
考虑修改
现在问题是：单点修改，区间求和
修改次数只有$ 10000$次，而求和次数达到$ 2*10^7$次
考虑分块
对于每个块我们维护块内前缀和以及前$ i$个块的和
这样查询是$ O(1)$的而修改是$ O(2000)$的
就以一种非常优秀的均摊复杂度解决了这道题
$ my \ code$
#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#define p 1000000007
#define rt register int
#define ll long long
using namespace std;
inline ll read(){
ll x = ; char zf = ; char ch = getchar();
while (ch != '-' && !isdigit(ch)) ch = getchar();
if (ch == '-') zf = -, ch = getchar();
while (isdigit(ch)) x = x * + ch - '', ch = getchar(); return x * zf;
}
void write(ll y){if(y<)putchar('-'),y=-y;if(y>)write(y/);putchar(y%+);}
void writeln(const ll y){write(y);putchar('\n');}
int i,j,k,m,n,x,y,z,cnt,blo;
int phi[],ss[],inv[],f[];bool pri[];
int val[],qz1[][],qz2[];
void init(int M){
//mu[1]=1;
phi[]=;blo=(int)sqrt(M);
for(rt i=;i<=M;i++){
if(!pri[i])phi[i]=i-,ss[++cnt]=i;
for(rt j=;j<=cnt&&i*ss[j]<=M;j++){
pri[i*ss[j]]=;
phi[i*ss[j]]=phi[i]*phi[ss[j]];
if(i%ss[j]==){
phi[i*ss[j]]=phi[i]*ss[j];
break;
}
}
}
for(rt i=;i<=M;i++)val[i]=(val[i-]+1ll*i*i%p*phi[i]%p)%p;
for(rt i=;i<=M;i++){
f[i]=1ll*i*i%p;
if((i-)%blo==)qz1[(i-)/blo][]=f[i];
else qz1[(i-)/blo][(i-)%blo]=(qz1[(i-)/blo][(i-)%blo-]+f[i])%p;
}
qz2[]=qz1[][blo-];
for(rt i=;i<=blo;i++)qz2[i]=(qz2[i-]+qz1[i][blo-])%p;
inv[]=inv[]=;
for(rt i=;i<=M;i++)inv[i]=1ll*inv[p%i]*(p-p/i)%p;
}
int qz(int x){
if(!x)return ;int ks=(x-)/blo;
if(ks==)return qz1[][x-];
else return (qz1[ks][(x-)%blo]+qz2[ks-])%p;
}
int calc(int n){
int ans=;
for(rt i=;i<=n;){
int R=n/(n/i);
(ans+=1ll*(qz(R)-qz(i-))*val[n/i]%p)%=p;
i=R+;
}
return (ans+p)%p;
}
void update(int x,int y){
int ks=(x-)/blo;
for(rt i=x;(i-)%blo!=||i==x;i++)(qz1[ks][(i-)%blo]+=y)%=p;
for(rt i=ks;i<=blo;i++)(qz2[i]+=y)%=p;
(f[x]+=y)%=p;
}
int main(){
n=read();m=read();
init(m);
for(rt i=;i<=n;i++){
x=read();y=read();int v=read()%p;k=read();
int gc=__gcd(x,y);
v=1ll*v*inv[x]%p*inv[y]%p*gc%p*gc%p;
update(gc,v-f[gc]);
writeln(calc(k));
}
return ;
}