来自FallDream的博客，未经允许，请勿转载，谢谢。
小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之后的结果。
每次修改操作把(a,b)改成x并且求前k行k列的和
操作数量m<=10000 n,k,a,b<=4*10^6 x<=10^18
首先从条件入手，发现很像辗转相除法。仔细观察发现，f(a,b)总是和f(g,g)( g=gcd(a,b) )有关系。更详细地，g(a,b)=f(g,g)*a/g*b/g
所以只需要几下主对角线的数字即可，考虑计算答案。以下的n表示询问的k,且num(x)表示f(x,x)
枚举gcd是啥
$$Ans=sum_{g=1}^{n}num(g)*sum_{i=1}^{lfloorfrac{n}{g}
floor}sum_{j=1}^{lfloorfrac{n}{g}
floor}ijg^{2}*[gcd(i,j)==1]$$
当然，把后面那一坨提出来比较舒服，发现可以用phi来化简
$$G(n)=sum_{i=1}^{n}sum_{j=1}^{n}i*j*[gcd(i,j)==1]$$
因为$$sum_{i=1}^{n}i*[gcd(i,n)==1]=frac{n*varphi(n)}{2}$$
所以$$G(n)=sum_{i=1}^{n}i^{2}varphi(i)$$
显然可以打表
然后这时候
$$Ans=sum_{g=1}^{n}num(g)*G(lfloorfrac{n}{g}
floor)$$
$lfloorfrac{n}{g}
floor$只有根号种取值，所以只需要维护前面那东西的前缀和就行了
但是每次查询必须是O(1)的，很自然想到分块维护前缀和，修改的时候直接修改gcd即可。这样就做完啦。
复杂度是$O(msqrt{n})$
强行写了一个llread返回个int查了好久错...心塞
#include<iostream> #include<cstdio> #include<cmath> #define ll long long #define MN 4000000 #define MB 2000 #define mod 1000000007 using namespace std; inline int read() { int x = 0; char ch = getchar(); while(ch < '0' || ch > '9')ch = getchar(); while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0';ch = getchar();} return x; } inline ll llread() { ll x = 0 ; char ch = getchar(); while(ch < '0' || ch > '9') ch = getchar(); while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0';ch = getchar();} return x; } int n,m,phi[MN+5],s[MN],cnt=0,la,block,add[MB+5]; int num[MN+5]; bool b[MN+5]; inline int gcd(int x,int y) {return !y?x:gcd(y,x%y);} void Modify(int x,int ad) { int bl=(x-1)/block+1,M=min(n,bl*block); for(int j=bl+1;j<=la;++j) (add[j]+=ad)%=mod; for(int j=x;j<=M;++j) (num[j]+=ad)%=mod; } int Query(int x) { if(!x) return 0; int bl=(x-1)/block+1; return (num[x]+add[bl])%mod; } int main() { m=read();n=read();num[1]=phi[1]=1;block=sqrt(n);la=(n-1)/block+1; for(int i=2;i<=n;++i) { if(!b[i]) phi[s[++cnt]=i]=i-1; for(int j=1;s[j]*i<=n;++j) { b[s[j]*i]=1; if(i%s[j]==0){ phi[s[j]*i]=phi[i]*s[j];break;} phi[s[j]*i]=phi[i]*(s[j]-1); } phi[i]=(phi[i-1]+1LL*i*i%mod*phi[i])%mod; num[i]=(num[i-1]+1LL*i*i)%mod; } for(int i=1;i<=m;++i) { int x=read(),y=read();ll X=llread();int k=read(); int g=gcd(x,y),ans=0;X/=1LL*(x/g)*(y/g);X%=mod; Modify(g,((X-Query(g)+mod)%mod+Query(g-1))%mod); for(int j=1,last;j<=k;j=last+1) { last=k/(k/j); ans=(ans+1LL*(Query(last)-Query(j-1)+mod)%mod*phi[k/j])%mod; } printf("%d ",ans); } return 0; }