This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub mugen1337/procon
#include "Math/ManhattanSquareSum.hpp"
#include "../type/modint.hpp" #include "../SegmentTree/SegmentTree.hpp" using mint=ModInt<998244353>; /* https://yukicoder.me/problems/no/1649 式変形して展開すると,sum(sum(|x_i - x_j| |y_i - y_j|))が計算できればよいとわかる. これはxが大きい順から処理することでx側の絶対値を外し,y側をセグメント木で管理し 平面走査のような感じで解ける */ mint ManhattanSquareSum(vector<ll> x,vector<ll> y){ int N=(int)x.size(); mint ret=0; { mint sx=0,sy=0,sxx=0,syy=0; for(int i=0;i<N;i++){ ret+=(mint)x[i]*x[i]*i+sxx-sx*x[i]*2; ret+=(mint)y[i]*y[i]*i+syy-sy*y[i]*2; sx+=x[i],sy+=y[i]; sxx+=(mint)x[i]*x[i]; syy+=(mint)y[i]*y[i]; } } { vector<int> idx(N); iota(begin(idx),end(idx),0); sort(begin(idx),end(idx),[&](int i,int j){return x[i]>x[j];}); auto nx=x,ny=y; for(int i=0;i<N;i++) nx[i]=x[idx[i]],ny[i]=y[idx[i]]; swap(x,nx); swap(y,ny); } { vector<int> yid(N),yp(N); iota(begin(yid),end(yid),0); sort(begin(yid),end(yid),[&](int i,int j){return y[i]<y[j];}); for(int i=0;i<N;i++) yp[yid[i]]=i; using P=pair<mint,mint>; auto segf=[&](P a,P b){ a.first+=b.first; a.second+=b.second; return a; }; SegmentTree<P> sega(N,segf,P(0,0)),segb(N,segf,P(0,0)); for(int i=0;i<N;i++){ sega.set(yp[i],P(y[i],1)); segb.set(i,P(0,0)); } sega.build(); segb.build(); for(int i=0;i<N;i++){ int j=yp[i]; auto l=sega.query(0,j); auto r=sega.query(j+1,N); mint add=mint(y[i])*l.second-l.first+r.first-mint(y[i])*r.second; l=segb.query(0,j); r=segb.query(j+1,N); mint sub=mint(y[i])*l.second-l.first+r.first-mint(y[i])*r.second; ret+=(add-sub)*x[i]*2; sega.update(j,P(0,0)); segb.update(j,P(y[i],1)); } } return ret; }
#line 1 "type/modint.hpp" template<int Mod> struct ModInt{ int x; ModInt():x(0){} ModInt(int y): x (y >= 0 ? y % Mod : (Mod - (-y) % Mod) % Mod){} ModInt(long long y){ if (y >= 0) { x = (int)(y % (ll)(Mod)); } else { int tmp = (int)((-y) % (ll)Mod); x = (Mod - tmp) % Mod; } } ModInt &operator+=(const ModInt &p){ if((x += p.x) >= Mod) x -= Mod; return *this; } ModInt &operator-=(const ModInt &p){ if((x += Mod - p.x) >= Mod) x -= Mod; return *this; } ModInt &operator*=(const ModInt &p){ x = (int)(1ll * x * p.x % Mod); return *this; } ModInt &operator/=(const ModInt &p){ (*this) *= p.inverse(); return *this; } ModInt operator-()const{return ModInt(-x);} ModInt operator+(const ModInt &p)const{return ModInt(*this)+=p;} ModInt operator-(const ModInt &p)const{return ModInt(*this)-=p;} ModInt operator*(const ModInt &p)const{return ModInt(*this)*=p;} ModInt operator/(const ModInt &p)const{return ModInt(*this)/=p;} bool operator==(const ModInt &p)const{return x==p.x;} bool operator!=(const ModInt &p)const{return x!=p.x;} ModInt inverse()const{ int a = x, b = Mod ,u = 1, v = 0, t; while(b>0){ t=a/b; swap(a-=t*b,b);swap(u-=t*v,v); } return ModInt(u); } ModInt pow(long long n)const{ ModInt ret(1),mul(x); while(n>0){ if(n&1) ret*=mul; mul*=mul;n>>=1; } return ret; } friend ostream &operator<<(ostream &os,const ModInt &p){return os<<p.x;} friend istream &operator>>(istream &is,ModInt &a){long long t;is>>t;a=ModInt<Mod>(t);return (is);} static int get_mod(){return Mod;} }; #line 1 "SegmentTree/SegmentTree.hpp" template<typename Monoid> struct SegmentTree{ using F=function<Monoid(Monoid,Monoid)>; private: int sz; vector<Monoid> seg; Monoid query(int a,int b,int k,int l,int r){ if(a<=l and r<=b) return seg[k]; if(b<=l or r<=a) return M0; Monoid L=query(a,b,2*k,l,(l+r)/2); Monoid R=query(a,b,2*k+1,(l+r)/2,r); return f(L,R); } template<typename C> int find_first(int a,const C &check,Monoid &acc,int k,int l,int r){ if(k>=sz){ acc=f(acc,seg[k]); if(check(acc)) return -1; else return k-sz; } int m=(l+r)/2; if(m<=a) return find_first(a,check,acc,2*k+1,m,r); if(a<=l and check(f(acc,seg[k]))){ acc=f(acc,seg[k]); return -1; } int x=find_first(a,check,acc,2*k+0,l,m); if(x>=0) return x; return find_first(a,check,acc,2*k+1,m,r); } template<typename C> int find_last(int b,const C &check,Monoid &acc,int k,int l,int r){ if(k>=sz){ acc=f(acc,seg[k]); if(check(acc)) return -1; else return k-sz+1;//ここはfalse, +1した位置はtrue } int m=(l+r)/2; if(b<=m) return find_last(b,check,acc,2*k,l,m); if(r<=b and check(f(acc,seg[k]))){ acc=f(acc,seg[k]); return -1; } int x=find_last(b,check,acc,2*k+1,m,r); if(x>=0) return x; return find_last(b,check,acc,2*k,l,m); } public: F f; Monoid M0;// モノイドの元 SegmentTree(int n, F f_, Monoid M0_) : f(f_), M0(M0_) { sz=1; while(sz<n)sz<<=1; seg.assign(2*sz,M0); } void set(int k,Monoid x){ seg[k+sz]=x; } void build(){ for(int k=sz-1;k>0;k--) seg[k]=f(seg[2*k],seg[2*k+1]); } void update(int k,Monoid x){ k+=sz; seg[k]=x; k>>=1; for(;k;k>>=1) seg[k]=f(seg[2*k],seg[2*k+1]); } Monoid query(int a,int b){ return query(a,b,1,0,sz); } Monoid operator[](const int &k)const{ return seg[k+sz]; } // http://codeforces.com/contest/914/submission/107505449 // max x, check(query(a, x)) = true template<typename C> int find_first(int a,const C &check){ Monoid val=M0; return find_first(a,check,val,1,0,sz); } // http://codeforces.com/contest/914/submission/107505582 // min x, check(query(x, b)) = true template<typename C> int find_last(int b,C &check){ Monoid val=M0; return find_last(b,check,val,1,0,sz); } }; #line 3 "Math/ManhattanSquareSum.hpp" using mint=ModInt<998244353>; /* https://yukicoder.me/problems/no/1649 式変形して展開すると,sum(sum(|x_i - x_j| |y_i - y_j|))が計算できればよいとわかる. これはxが大きい順から処理することでx側の絶対値を外し,y側をセグメント木で管理し 平面走査のような感じで解ける */ mint ManhattanSquareSum(vector<ll> x,vector<ll> y){ int N=(int)x.size(); mint ret=0; { mint sx=0,sy=0,sxx=0,syy=0; for(int i=0;i<N;i++){ ret+=(mint)x[i]*x[i]*i+sxx-sx*x[i]*2; ret+=(mint)y[i]*y[i]*i+syy-sy*y[i]*2; sx+=x[i],sy+=y[i]; sxx+=(mint)x[i]*x[i]; syy+=(mint)y[i]*y[i]; } } { vector<int> idx(N); iota(begin(idx),end(idx),0); sort(begin(idx),end(idx),[&](int i,int j){return x[i]>x[j];}); auto nx=x,ny=y; for(int i=0;i<N;i++) nx[i]=x[idx[i]],ny[i]=y[idx[i]]; swap(x,nx); swap(y,ny); } { vector<int> yid(N),yp(N); iota(begin(yid),end(yid),0); sort(begin(yid),end(yid),[&](int i,int j){return y[i]<y[j];}); for(int i=0;i<N;i++) yp[yid[i]]=i; using P=pair<mint,mint>; auto segf=[&](P a,P b){ a.first+=b.first; a.second+=b.second; return a; }; SegmentTree<P> sega(N,segf,P(0,0)),segb(N,segf,P(0,0)); for(int i=0;i<N;i++){ sega.set(yp[i],P(y[i],1)); segb.set(i,P(0,0)); } sega.build(); segb.build(); for(int i=0;i<N;i++){ int j=yp[i]; auto l=sega.query(0,j); auto r=sega.query(j+1,N); mint add=mint(y[i])*l.second-l.first+r.first-mint(y[i])*r.second; l=segb.query(0,j); r=segb.query(j+1,N); mint sub=mint(y[i])*l.second-l.first+r.first-mint(y[i])*r.second; ret+=(add-sub)*x[i]*2; sega.update(j,P(0,0)); segb.update(j,P(y[i],1)); } } return ret; }