题目描述
小Q是个程序员。
作为一个年轻的程序员，小Q总是被老C欺负，老C经常把一些麻烦的任务交给小Q来处理。每当小Q不知道如何解决时，就只好向你求助。
为了完成任务，小Q需要列一个表格，表格有无穷多行，无穷多列，行和列都从1开始标号。为了完成任务，表格里面每个格子都填了一个整数，为了方便描述，小Q把第a行第b列的整数记为f(a,b)。为了完成任务，这个表格要满足一些条件：
(1)对任意的正整数a,b，都要满足f(a,b)=f(b,a)；
(2)对任意的正整数a,b，都要满足b×f(a,a+b)=(a+b)×f(a,b)。
为了完成任务，一开始表格里面的数很有规律，第a行第b列的数恰好等于a×b，显然一开始是满足上述两个条件的。为了完成任务，小Q需要不断的修改表格里面的数，每当修改了一个格子的数之后，为了让表格继续满足上述两个条件，小Q还需要把这次修改能够波及到的全部格子里都改为恰当的数。由于某种神奇的力量驱使，已经确保了每一轮修改之后所有格子里的数仍然都是整数。为了完成任务，小Q还需要随时获取前k行前k列这个有限区域内所有数的和是多少，答案可能比较大，只需要算出答案mod1,000,000,007之后的结果。
题解
这题思路太神了。
乍一看需要维护一个二维的东西看起来很麻烦，但实际上是可以转化为一维上的问题的。
先来看第一个限制，关于对角线对称的数字相等，这样我们的棋盘大小变为原来的一半。
然后看第二个限制，挪一下位置，变成了
f(a,a+b)/(a+b)=f(a,b)/b
f(a,a+b)/((a+b)*a)=f(a,b)/(a*b)
我们发现了规律，f(a,b)是和a*b正相关的，而且这种关系非常像辗转相减，那么一直减到最后就会变成f(gcd,gcd)。
于是我们得出了结论，每个位置(a,b)的值之和gcd(a,b)相关。
于是它变成了一个序列上的问题，每次修改只需要修改一个数。
然后考虑答案。
ans=∑f(g,g)∑∑i*j*(gcd(i,j)==g) i,j,g<=k
后面的那个东西∑∑i*j*(gcd(i,j)==g)经过推导后可以得到∑i*i*Φ(i) i<=k/g
然后就可以除法分块了。
但我们还需要支持单点修改，区间求和，可以用树状数组，但是这题的数据范围比较特别，所以用分块搞就可以了，根号修改，O(1)查询。
代码
#include<iostream> #include<cstdio> #include<cmath> #define N 4000009 using namespace std; typedef long long ll; const int mod=1e9+7; ll n,m,sum[N],pre[N],phi[N],prime[N],f[N],val[N],n1,be[N]; bool vis[N]; inline ll rd(){ ll x=0;char c=getchar();bool f=0; while(!isdigit(c)){if(c=='-')f=1;c=getchar();} while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();} return f?-x:x; } inline void work(ll x,ll y){ pre[x]=(pre[x]+y)%mod; for(int i=x+1;be[i]==be[i-1];++i)pre[i]=(pre[i]+y)%mod; for(int i=be[x];i<=be[n];++i)sum[i]=(sum[i]+y)%mod; } ll gcd(ll a,ll b){return b?gcd(b,a%b):a;} inline ll query(int x){return (pre[x]+sum[be[x]-1])%mod;} inline void prework(){ phi[1]=1; for(int i=2;i<=n;++i){ if(!vis[i])prime[++prime[0]]=i,phi[i]=i-1; for(int j=1;j<=prime[0]&&i*prime[j]<=n;++j){ vis[i*prime[j]]=1; if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;} phi[i*prime[j]]=phi[i]*phi[prime[j]]; } } for(int i=1;i<=n;++i)f[i]=(f[i-1]+phi[i]*val[i])%mod; } int main(){ m=rd();n=rd();n1=sqrt(n); for(int i=1;i<=n;++i){ be[i]=(i-1)/n1+1;val[i]=1ll*i*i%mod; (sum[be[i]]+=val[i])%=mod; if(be[i]==be[i-1])pre[i]=(pre[i-1]+val[i])%mod; else pre[i]=val[i]; } for(int i=2;i<=be[n];++i)(sum[i]+=sum[i-1])%=mod; prework();ll a,b,x,k; while(m--){ a=rd();b=rd();x=rd();k=rd();ll g=gcd(a,b);ll ans=0,r; x=x/((a/g)*(b/g));x%=mod;work(g,(x-val[g]+mod)%mod);val[g]=x; for(ll l=1;l<=k;l=r+1){ r=k/(k/l); ans+=(query(r)-query(l-1))*f[k/l]%mod; ans=(ans%mod+mod)%mod; } printf("%lld ",ans); } return 0; }