https://blog.sengxian.com/solutions/bzoj-2877
注意二维线段树的upd()也是一个O(log n)的函数（pushdown()应该也是但没写过）。
1 #include<cstdio> 2 #include<algorithm> 3 #define Ls (x<<1) 4 #define Rs (Ls|1) 5 #define Lson Ls,L,mid 6 #define Rson Rs,mid+1,R 7 #define lson ls[x],L,mid 8 #define rson rs[x],mid+1,R 9 #define rep(i,l,r) for (int i=(l); i<=(r); i++) 10 typedef long long ll; 11 using namespace std; 12 13 const int N=500010,M=10000010; 14 ll c,sm[M],a[N]; 15 int n,m,x,y,T,nd,op,x1,x2,y1,y2,rt[N<<2],ls[M],rs[M]; 16 17 int F(int x,int y){ return (x-1)*m+y; } 18 ll gcd(ll a,ll b){ return b ? gcd(b,a%b) : a; } 19 20 void mdfy(int &x,int L,int R,int pos,ll k){ 21 if (!x) x=++nd; 22 if (L==R){ sm[x]+=k; return; } 23 int mid=(L+R)>>1; 24 if (pos<=mid) mdfy(lson,pos,k); else mdfy(rson,pos,k); 25 sm[x]=gcd(sm[ls[x]],sm[rs[x]]); 26 } 27 28 void upd(int &x,int L,int R,int pos,int lp,int rp){ 29 if (!x) x=++nd; 30 if (L==R){ sm[x]=gcd(sm[lp],sm[rp]); return; } 31 int mid=(L+R)>>1; 32 if (pos<=mid) upd(lson,pos,ls[lp],ls[rp]); else upd(rson,pos,rs[lp],rs[rp]); 33 sm[x]=gcd(sm[ls[x]],sm[rs[x]]); 34 } 35 36 void mdfx(int x,int L,int R,int x1,int y1,ll k){ 37 if (L==R){ mdfy(rt[x],0,m,y1,k); return; } 38 int mid=(L+R)>>1; 39 if (x1<=mid) mdfx(Lson,x1,y1,k); else mdfx(Rson,x1,y1,k); 40 upd(rt[x],0,m,y1,rt[Ls],rt[Rs]); 41 } 42 43 ll quey(int x,int L,int R,int l,int r){ 44 if (L==l && r==R) return sm[x]; 45 int mid=(L+R)>>1; 46 if (r<=mid) return quey(lson,l,r); 47 else if (l>mid) return quey(rson,l,r); 48 else return gcd(quey(lson,l,mid),quey(rson,mid+1,r)); 49 } 50 51 ll quex(int x,int L,int R,int l,int r,int y1,int y2){ 52 if (L==l && r==R) return quey(rt[x],0,m,y1,y2); 53 int mid=(L+R)>>1; 54 if (r<=mid) return quex(Lson,l,r,y1,y2); 55 else if (l>mid) return quex(Rson,l,r,y1,y2); 56 else return gcd(quex(Lson,l,mid,y1,y2),quex(Rson,mid+1,r,y1,y2)); 57 } 58 59 int main(){ 60 freopen("chess.in","r",stdin); 61 freopen("chess.out","w",stdout); 62 scanf("%d%d%d%d%d",&n,&m,&x,&y,&T); int r1=(n+1)*4+1,r2=r1+1; 63 rep(i,1,n) rep(j,1,m) scanf("%lld",&a[F(i,j)]); 64 rep(i,1,n-1) rep(j,1,m-1) mdfx(1,0,n,i,j,a[F(i+1,j+1)]-a[F(i+1,j)]-a[F(i,j+1)]+a[F(i,j)]); 65 rep(i,1,n-1) mdfy(rt[r1],0,n,i,a[F(i+1,y)]-a[F(i,y)]); 66 rep(i,1,m-1) mdfy(rt[r2],0,m,i,a[F(x,i+1)]-a[F(x,i)]); 67 while (T--){ 68 scanf("%d%d%d%d%d",&op,&x1,&y1,&x2,&y2); 69 if (op==0){ 70 ll t1=(x1+x2) ? quey(rt[r1],0,n,x-x1,x+x2-1) : 0; 71 ll t2=(y1+y2) ? quey(rt[r2],0,m,y-y1,y+y2-1) : 0; 72 ll t3=((x1+x2)&&(y1+y2)) ? quex(1,0,n,x-x1,x+x2-1,y-y1,y+y2-1) : 0; 73 printf("%lld ",abs(gcd(gcd(t1,t2),gcd(t3,a[F(x,y)])))); 74 }else{ 75 scanf("%lld",&c); 76 mdfx(1,0,n,x1-1,y1-1,c); mdfx(1,0,n,x1-1,y2,-c); 77 mdfx(1,0,n,x2,y1-1,-c); mdfx(1,0,n,x2,y2,c); 78 if (y1<=y && y2>=y) mdfy(rt[r1],0,n,x1-1,c),mdfy(rt[r1],0,n,x2,-c); 79 if (x1<=x && x2>=x) mdfy(rt[r2],0,m,y1-1,c),mdfy(rt[r2],0,m,y2,-c); 80 if (x1<=x && x2>=x && y1<=y && y2>=y) a[F(x,y)]+=c; 81 } 82 } 83 return 0; 84 }