注意二维线段树的upd()也是一个O(log n)的函数（pushdown()应该也是但没写过）。

1 #include
     
       2 #include
      
        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\n",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 }