動的平方分割(sumlt)
要素の追加削除ができる平方分割について考える
Nを最大入力件数、bs(バケットサイズ)を√Nで固定、バケットへの要素の追加削除がO(√N)、バケットの追加削除がO(N)とする
添字0,1,2のバケット3個が使われているとする
→ (0,1,2)
要素の追加でバケット1が満杯になったら、新バケット3を追加してそこに1の後ろ半分の要素を移す
→ (0,1,3,2)
どのバケットもbsの1/3は使用されているものとする(1/4とかでもよい)
要素の削除によってバケットの使用量が1/3bs未満になったら、左右のバケットの使用量が1/3bsより多ければ要素をもらう、1/3bs以下なら2つのバケットを1つにマージする
バケットの追加によってできた2つのバケットはともに使用量が1/2bsなので、次のバケット追加までに要素の追加に対して1/2bs回の余裕が、次のバケット削除までに要素の削除に対して1/6bs回の余裕ができる。また、バケットの削除によってできた1つのバケットは使用量が2/3bsなので、次のバケットの追加削除までに要素の追加削除に対して1/3bs回の余裕ができる。このことから、バケットの追加削除にO(N)かけても、要素の追加削除はならし計算量というのでO(√N)といえそうなのでここではこれを動的平方分割とする。
静的WM(ウェーブレット行列)でできることは平方分割でO(√n)でできる
要素数nの要素列を平方分割する
bs=√n (nが平方数でない場合はbsなどを適当に調整する)
要素列(1次元)と平方分割(2次元)の表記を小文字と大文字で区別する
aやAなどのサイズの宣言は、a[n]、A[bs][bs]、となる
A[i][j] = Aのi-thバケットのj-th要素 = a[i*bs+j]
a[i] = A[i/bs][i mod bs]
a : 要素数nの要素列
sa : aをソート(安定ソート)した列
b[i] : sa[i]のソート前の位置 ( a[b[i]] = sa[i] )
C[i][j] : Aの0-thからi-thバケットまでの要素のうちSAの0-thからj-thバケットまでに含まれる要素の個数
d[i] : a[i]のソート後の位置 ( a[i] = sa[d[i]] )
rankとselect
saとbを二分探索するとできる
quantile
quantile(l,r,k) : aの[l,r)をソートしたときのk番目に位置する要素
クエリ範囲[l,r)は、A[L],A[L+1],...,A[R-1]バケットと、その左右からはみ出した部分(両端と呼ぶことにする)からなるとする
C[i]などのバケットをベクトルとして考えてC[i]+C[j]などをベクトルどうしの足し算とする
サイズ√nの作業領域wを0で初期化
for each x in d列の両端部分{
w[x/bs]+=1
}
for (i=1;i<bs;++i){
w[i]+=w[i-1]
}
w+=C[R-1]-C[L-1]
とするとw[i]は(r-l)個のうちSA[i]までに含まれる要素数になるのでk番目がSAのどのバケットにあるか分かるのでそこから探索する
ranklt(rank less than)
ranklt(l,r,x) : aの[l,r)に含まれるx未満の要素の個数
[l,r)はバケットちょうど何個分としてよいので(両端部分やクエリ範囲が1バケット内に含まれる場合は全探索する)、クエリ範囲はA[L],A[L+1],...,A[R-1]バケットとする
SAでx以上の最小値が存在するバケットをSA[i]とする
ranklt(l,r,x) = C[R-1][i-1] - C[L-1][i-1] + SA[i]内でx未満かつbが[l,r)を満たす要素の個数
Cの作り方の一例を示す(Cは0で初期化しておく)
// step1
for (i=0;i<n;++i){
j=b[i]/bs
k=i/bs
C[j][k]+=1
}
// step2(バケット間累積和？)
for (i=1;i<bs;++i){
for (j=0;j<bs;++j){
C[i][j]+=C[i-1][j]
}
}
// step3(バケット内累積和？)
for (i=0;i<bs;++i){
for (j=1;j<bs;++j){
C[i][j]+=C[i][j-1]
}
}
step1のC[j][k]+=1をC[j][k]+=sa[i]とするなどの多少の変更でrankltのコードはsumltになる
sumlt(l,r,x) : l<=i<r かつ a[i]<x である a[i]の和
また、aの付随列a2があるときstep1のC[j][k]+=1をC[j][k]+=a2[b[i]]とするなどの多少の変更でrankltのコードはsumlt2になる
sumlt2(l,r,x) : l<=i<r かつ a[i]<x である a2[i]の和
(ここでいう「+」とか「和」というのは、(整数,足し算)に限らず、ソートできる集合Sがあって(S,Sの上の2項演算)がアーベル群の場合の2項演算とその計算結果のこととしてたぶん一般化できる)
各バケットを再度平方分割することで、rankltはO(pow(n,1/4))に、sumltも付随列のsumltの考え方を使うとO(pow(n,1/4))になる(pow(x,y)=xのy乗)
平方分割は何がうれしいのか
静的動的共通の話として、領域計算量ないし実際のメモリ使用量が他のデータ構造より小さくなる可能性がある(WMは補助データを使うときのメモリの大きさが実用上の問題になりうる)
複数のバケットを同時に処理したり、バケット(ベクトル)単位での処理が効率的にできるようになると強くなる可能性がある
動的平方分割によるsumlt
動的な場合、バケットの添字と位置は冒頭の(0,1,3,2)のように必ずしも一致しないが面倒なのでここでは一致しているものとする
insert(k,x)で要素xがaのk番目に追加されたとする
追加先をA[I]とSA[J]バケットとする(Jは(x,I)をキーにした二分探索で見つける)
(xの追加位置がSA[J][K]のとき、IをB[J][K]に追加する)
このとき静的な場合のC(rankltをsumltに変更したC)で考えると次のように更新する必要がありO(n)かかる
I<=i<bs かつ J<=j<bs について C[i][j]+=x
BIT(Binary Indexed Tree)のノードがBITを持つ2次元BITを使う
C[i]というバケットをBITのiノード上に置いて、そのバケットもBITで表現するみたいなこと
for (i=I;i<bs;i=iの親){
for (j=J;j<bs;j=jの親){
BIT[i][j]+=x
}
}
O(n)版では更新開始バケットの右側にある全バケットを更新しているが、BIT版は開始バケット(ノード)の先祖しか更新しないのでO(pow(log(n),2))になる
A側のバケットや(SAとB)側のバケットが追加削除されたらこのBITもO(N)で再構築する
静的な場合のCがBITになっているところ以外は動的な場合も考え方はだいたい同じ
バケットサイズをpow(N,2/3)、バケットの個数をpow(N,1/3)として各バケットを動的平方分割するとinsert、remove、sumltがならし計算量というのでO(pow(N,1/3))になるかもしれない
下のコードについて
動的平方分割による(整数,足し算)の場合のsumltで入力値(int)の和がlong longを超えない範囲で使える
上のCの定義のAとSAを入れ換えた定義にもとづいたBITを使っている(考え方に大差はないが再構築時のランダムアクセスを減らせるため)
#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <algorithm>
#define ll long long
struct bucketa{
int *a,n,hd;
};
struct bucketb{
int *b,*sa,n,hd;
ll w;
};
class BinaryIndexedTree{
private:
int *BIT,BITsize;
public:
int calc_mem(int n){
for (BITsize=2;BITsize<n;BITsize*=2);
return BITsize+1;
}
void set_mem(int *p){
BIT=p;
}
void add(int i,int x){
for (;i<=BITsize;i+=i&-i) BIT[i]+=x;
}
void addl(int i){
int j=++i+1;
int k=i+(i&-i);
for (;i!=k;i+=i&-i) ++BIT[i];
for (;j!=k;j+=j&-j) --BIT[j];
}
int sum(int i){
int s=0;
for (;i;i-=i&-i) s+=BIT[i];
return s;
}
int bsearch(int *y){
int i=BITsize>>1;
int x=*y;
for (int t=i>>1;t;t>>=1){
if (x<=BIT[i]) i-=t;
else x-=BIT[i],i+=t;
}
if (BIT[i]<x) *y=x-BIT[i],++i;
else *y=x;
return i;
}
void build(bucketa *a,int *idxa,int na){
calc_mem(na);
for (int i=0;i<na;++i) BIT[i+1]=a[idxa[i]].n;
for (int i=na+1;i<=BITsize;++i) BIT[i]=0;
for (int i=1;i<BITsize;++i) BIT[i+(i&-i)]+=BIT[i];
}
};
class BinaryIndexedTree2d{
private:
ll *a0;
int db,na,nb;
void BIT_build(ll *b){
int i,j;
for (i=1;i<na;++i){
j=i+(i&-i);
if (j<=na) b[j]+=b[i];
}
}
void csum(ll* const p,ll* const q,const int n_){
for (int i=0;i<n_;++i) p[i]+=q[i];
}
void csum_s(ll* const p,ll* const q,const int n_){
for (int i=0;i<n_;++i) p[i]-=q[i];
}
public:
int calc_mem(int argna){
db=argna+1;
return db*db;
}
void calc_mem(int argna,int argnb){
na=argna;
nb=argnb;
}
void set_mem(ll *p){
a0=p;
na=nb=1;
}
void add(int ib,int i,ll x){
for (ll *b=a0+(++ib)*db;ib<=nb;ib+=ib&-ib){
for (int j=i+1;j<=na;j+=j&-j) b[j]+=x;
b+=(ib&-ib)*db;
}
}
void addl_a(int ib,int i,ll x){
int rb=++ib+1;
int j,k=ib+(ib&-ib);
if (k>nb) k=nb+1;
for (ll *b=a0+ib*db;ib<k;ib+=ib&-ib){
for (j=i+1;j<=na;j+=j&-j) b[j]+=x;
b+=(ib&-ib)*db;
}
for (ll *b=a0+rb*db;rb<k;rb+=rb&-rb){
for (j=i+1;j<=na;j+=j&-j) b[j]-=x;
b+=(rb&-rb)*db;
}
}
void addl_b(int ib,int i,ll x){
++i;
int j,k=i+(i&-i);
if (k>na) k=na+1;
for (ll *b=a0+(++ib)*db;ib<=nb;ib+=ib&-ib){
for (j=i;j<k;j+=j&-j) b[j]+=x;
for (j=i+1;j<k;j+=j&-j) b[j]-=x;
b+=(ib&-ib)*db;
}
}
ll sum(int ib,int i){
ll s=0LL;
for (ll *b=a0+(++ib)*db;ib;ib-=ib&-ib){
for (int j=i+1;j;j-=j&-j) s+=b[j];
b-=(ib&-ib)*db;
}
return s;
}
void rebuild(bucketb *b,int *idxb,int argna,int argnb,int argbsz){
int i,j,t,u,ib2,*p,*q,h,d;
ll *r;
calc_mem(argna,argnb);
memset(a0+db,0,db*nb*8);
d=db;
r=a0+(d+1);
for (i=0;i<nb;++i,r+=d){
ib2=idxb[i];
u=b[ib2].n;
h=b[ib2].hd;
t=(argbsz-h<=u)? argbsz-h:u;
u-=t;
p=b[ib2].b;
q=b[ib2].sa;
for (j=0;j<u;++j) r[p[j]]+=(ll)q[j];
p+=h;
q+=h;
for (j=0;j<t;++j) r[p[j]]+=(ll)q[j];
}
for (i=1;i<nb;++i){
j=i+(i&-i);
if (j<=nb) csum(a0+(j*d),a0+(i*d),d);
}
r=a0+d;
for (i=1;i<=nb;++i,r+=d) BIT_build(r);
}
void rebuild_b_ins(bucketb *b,int *idxb,int argna,int argnb,int argbsz,const int pos){
int i,j,t,u,ib2,*p,*q,h;
ll *r;
const int d=db;
for (i=nb-1;0<i;--i){
j=i+(i&-i);
if (pos<j && j<=nb) csum_s(a0+(j*d),a0+(i*d),na+1);
}
r=a0+nb*d;
t=(na+1)*8;
for (i=nb;pos<i;--i,r-=d) memcpy(r+d,r,t);
memset(a0+(pos+1)*d,0,t);
r=a0+((pos+1)*d+1);
ib2=idxb[pos];
u=b[ib2].n;
h=b[ib2].hd;
t=(argbsz-h<=u)? argbsz-h:u;
u-=t;
p=b[ib2].b;
q=b[ib2].sa;
for (j=0;j<u;++j) r[p[j]]+=(ll)q[j];
p+=h;
q+=h;
for (j=0;j<t;++j) r[p[j]]+=(ll)q[j];
BIT_build(a0+(pos+1)*d);
csum_s(a0+(pos+2)*d,a0+(pos+1)*d,na+1);
calc_mem(argna,argnb);
for (i=1;i<nb;++i){
j=i+(i&-i);
if (pos<j && j<=nb) csum(a0+(j*d),a0+(i*d),na+1);
}
}
void rebuild_b_rem(int argna,int argnb,const int pos){
int i,j,t;
ll *r;
for (i=nb-1;0<i;--i){
j=i+(i&-i);
if (pos<j && j<=nb) csum_s(a0+(j*db),a0+(i*db),na+1);
}
csum(a0+(pos+1)*db,a0+(pos+2)*db,na+1);
r=a0+(pos+2)*db;
t=(na+1)*8;
for (i=nb;pos+1<i;--i,r+=db) memcpy(r,r+db,t);
calc_mem(argna,argnb);
for (i=1;i<nb;++i){
j=i+(i&-i);
if (pos<j && j<=nb) csum(a0+(j*db),a0+(i*db),na+1);
}
}
void rebuild_a_ins(bucketb *b,int *idxb,int argna,int argnb,const int pos){
int i,j,t,u;
ll *q,*r,s;
const int d=db;
u=(na-pos)*8;
r=a0+d;
for (i=nb;i;--i){
for (j=na-1;0<j;--j){
t=j+(j&-j);
if (t<=na) r[t]-=r[j];
}
memmove(r+(pos+2),r+(pos+1),u);
r[pos+1]=0LL;
r+=d;
}
r=a0+d;
for (i=1;i<=nb;++i,r+=d){
j=idxb[i-1];
s=(r[pos+1]+=b[j].w);
r[pos+2]-=s;
if ((i+(i&-i))<=nb){
q=r+(i&-i)*d;
q[pos+1]+=s;
}
}
r=a0+d;
for (i=nb;i;--i){
for (j=1;j<argna;++j){
t=j+(j&-j);
if (t<=argna) r[t]+=r[j];
}
r+=d;
}
calc_mem(argna,argnb);
}
void rebuild_a_rem(int argna,int argnb,const int pos){
int i,j,t,u;
ll *r;
u=(na-(pos+2))*8;
r=a0+db;
for (i=nb;i;--i){
for (j=na-1;0<j;--j){
t=j+(j&-j);
if (t<=na) r[t]-=r[j];
}
r[pos+1]+=r[pos+2];
memmove(r+(pos+2),r+(pos+3),u);
r+=db;
}
r=a0+db;
for (i=nb;i;--i){
for (j=1;j<argna;++j){
t=j+(j&-j);
if (t<=argna) r[t]+=r[j];
}
r+=db;
}
calc_mem(argna,argnb);
}
};
class BucketDecomposition{
private:
int *idxa,*idxb,bsz,na,nb,ui_a,mi_a,ui_b,mi_b,crem;
ll *pseg;
bucketa *a;
bucketb *b;
BinaryIndexedTree2d c;
BinaryIndexedTree BITa;
void tu(int h,int n_,int *t,int *u){
*t=(bsz-h<=n_)? bsz-h:n_;
*u=n_-*t;
}
void popr_shiftr(int *p,int l,int r){
if (l<r){
memmove(p+l+1,p+l,(r-l)*4);
}else if (l>r){
int t=p[bsz-1];
memmove(p+l+1,p+l,(bsz-1-l)*4);
memmove(p+1,p,r*4);
p[0]=t;
}
}
void popl_shiftl(int *p,int l,int r){
if (l<r){
memmove(p+l,p+l+1,(r-l)*4);
}else if (l>r){
int t=p[0];
memmove(p,p+1,r*4);
memmove(p+l,p+l+1,(bsz-1-l)*4);
p[bsz-1]=t;
}
}
void array_insert(int *p,int i,int n,int x){
popr_shiftr(p,i,n);
p[i]=x;
}
void array_remove(int *p,int i,int n){
popl_shiftl(p,i,n-1);
}
void bucketa_insert(int ib,int i,int x){
int l=a[ib].hd;
int r=a[ib].n++;
if (i+i<r){
a[ib].hd=l=(l? l-1:bsz-1);
r=(l+i<bsz)? l+i:l+i-bsz;
popl_shiftl(a[ib].a,l,r);
a[ib].a[r]=x;
}else{
r=(l+r<bsz)? l+r:l+r-bsz;
l=(l+i<bsz)? l+i:l+i-bsz;
popr_shiftr(a[ib].a,l,r);
a[ib].a[l]=x;
}
}
int bucketa_remove(int ib,int i){
int l,r,x;
l=a[ib].hd;
r=--a[ib].n;
if (i+i<=r){
r=(l+i<bsz)? l+i:l+i-bsz;
a[ib].hd=(l+1<bsz)? l+1:0;
x=a[ib].a[r];
popr_shiftr(a[ib].a,l,r);
}else{
r=(l+r<bsz)? l+r:l+r-bsz;
l=(l+i<bsz)? l+i:l+i-bsz;
x=a[ib].a[l];
popl_shiftl(a[ib].a,l,r);
}
return x;
}
void bucketb_insert(int ib,int i,int xb,int xsa){
int l=b[ib].hd;
int r=b[ib].n++;
if (i+i<r){
b[ib].hd=l=(l? l-1:bsz-1);
r=(l+i<bsz)? l+i:l+i-bsz;
popl_shiftl(b[ib].b,l,r);
b[ib].b[r]=xb;
popl_shiftl(b[ib].sa,l,r);
b[ib].sa[r]=xsa;
}else{
r=(l+r<bsz)? l+r:l+r-bsz;
l=(l+i<bsz)? l+i:l+i-bsz;
popr_shiftr(b[ib].b,l,r);
b[ib].b[l]=xb;
popr_shiftr(b[ib].sa,l,r);
b[ib].sa[l]=xsa;
}
}
void bucketb_remove(int ib,int i){
int l=b[ib].hd;
int r=--b[ib].n;
if (i+i<=r){
r=(l+i<bsz)? l+i:l+i-bsz;
b[ib].hd=(l+1<bsz)? l+1:0;
popr_shiftr(b[ib].b,l,r);
popr_shiftr(b[ib].sa,l,r);
}else{
r=(l+r<bsz)? l+r:l+r-bsz;
l=(l+i<bsz)? l+i:l+i-bsz;
popl_shiftl(b[ib].b,l,r);
popl_shiftl(b[ib].sa,l,r);
}
}
void bucketa_merge(int lb,int ib){
int i,j,t,u,*p,*q,nl;
nl=a[lb].n;
i=a[lb].hd;
j=(i+nl)%bsz;
q=a[ib].a;
if (i<=j){
p=a[lb].a;
memmove(p,p+i,nl*4);
p+=nl;
a[lb].hd=0;
}else{
p=a[lb].a+j;
}
tu(a[ib].hd,a[ib].n,&t,&u);
memcpy(p,q+a[ib].hd,t*4);
memcpy(p+t,q,u*4);
}
void bucketb_merge(int lb,int ib){
int i,j,t,u,*p,*q,nl;
nl=b[lb].n;
i=b[lb].hd;
j=(i+nl)%bsz;
if (i<=j){
p=b[lb].b;
q=b[lb].sa;
memmove(p,p+i,nl*4);
memmove(q,q+i,nl*4);
p+=nl;q+=nl;
b[lb].hd=0;
}else{
p=b[lb].b+j;
q=b[lb].sa+j;
}
tu(b[ib].hd,b[ib].n,&t,&u);
memcpy(p,b[ib].b+b[ib].hd,t*4);
memcpy(p+t,b[ib].b,u*4);
memcpy(q,b[ib].sa+b[ib].hd,t*4);
memcpy(q+t,b[ib].sa,u*4);
}
int b_bsearch_pos(const int x,const int pos,int *k){
int i,j,m,l,r,ib,*p,*q,ret;
l=0;r=nb-1;
while (l<=r){
m=(l+r)>>1;
ib=idxb[m];
i=b[ib].hd;
j=(i+b[ib].n-1);
if (j>=bsz) j-=bsz;
p=b[ib].sa;
q=b[ib].b;
if ((x<p[i]) || (p[i]==x && pos<=q[i])){
r=m-1;
}else if ((p[j]<x) || (p[j]==x && pos>q[j])){
l=m+1;
}else{
l=m;break;
}
}
if (nb<=l) l=nb-1;
ret=l;
ib=idxb[l];
l=0;r=b[ib].n-1;
p=b[ib].sa;
q=b[ib].b;
i=b[ib].hd;
while (l<=r){
m=(l+r)>>1;
j=(i+m<bsz)? i+m:i+m-bsz;
if ((x<p[j]) || (p[j]==x && pos<=q[j])){
r=m-1;
}else{
l=m+1;
}
}
*k=l;
return ret;
}
int b_bsearch(const int x){
int i,j,m,l,r,ib,*p;
l=0;r=nb-1;
while (l<=r){
m=(l+r)>>1;
ib=idxb[m];
p=b[ib].sa;
i=b[ib].hd;
j=(i+b[ib].n-1);
if (j>=bsz) j-=bsz;
if (x<=p[i]){
r=m-1;
}else if (p[j]<x){
l=m+1;
}else{
l=m;
break;
}
}
if (l>=nb) l=nb-1;
return l;
}
int b_bsearch2(int ib,const int x){
int i,j,l,r,m,*p;
l=0;r=b[ib].n-1;
p=b[ib].sa;
i=b[ib].hd;
while (l<=r){
m=(l+r)>>1;
j=(i+m<bsz)? i+m:i+m-bsz;
if (x<=p[j]){
r=m-1;
}else{
l=m+1;
}
}
return l;
}
int getnewindex(int flga){
int i;
if (flga){
if (ui_a) i=ui_a,ui_a=a[i].n;
else i=++mi_a;
a[i].n=a[i].hd=0;
}else{
if (ui_b) i=ui_b,ui_b=b[i].n;
else i=++mi_b;
b[i].n=b[i].hd=0;
}
return i;
}
void a_rebuild_ins(const int pos,int ib){
int i,j,k,t,u,ib2,pos2,newib,*p,*q,hf;
for (i=0;i<nb;++i){
ib2=idxb[i];
b[ib2].w=0LL;
tu(b[ib2].hd,b[ib2].n,&t,&u);
q=b[ib2].b;
p=q+b[ib2].hd;
for (;u;++q,--u) if (pos<=*q) ++*q;
for (;t;++p,--t) if (pos<=*p) ++*p;
}
hf=bsz/2;
newib=getnewindex(1);
j=(a[ib].hd+hf)%bsz;
t=(j<=hf)? hf:bsz-j;
u=hf-t;
p=a[ib].a;
q=a[newib].a;
memcpy(q,p+j,t*4);
memcpy(q+t,p,u*4);
a[ib].n=a[newib].n=hf;
j=a[ib].hd;
for (i=0;i<hf;++i){
u=p[j];
pos2=b_bsearch_pos(u,pos+1,&k);
ib2=idxb[pos2];
t=(b[ib2].hd+k)%bsz;
b[ib2].b[t]=pos;
b[ib2].w+=(ll)u;
j=(j+1)%bsz;
}
array_insert(idxa,pos+1,na,newib);
++na;
BITa.build(a,idxa,na);
}
void b_rebuild_ins(int pos2,int ib2){
int i,newib,hf,*p,*q,t,u;
newib=getnewindex(0);
hf=bsz/2;
b[ib2].n=b[newib].n=hf;
i=(b[ib2].hd+hf)%bsz;
t=(i<=hf)? hf:bsz-i;
u=hf-t;
p=b[ib2].b;
q=b[newib].b;
memcpy(q,p+i,t*4);
memcpy(q+t,p,u*4);
p=b[ib2].sa;
q=b[newib].sa;
memcpy(q,p+i,t*4);
memcpy(q+t,p,u*4);
array_insert(idxb,pos2+1,nb,newib);
++nb;
}
void a_rebuild_rem(const int pos,int ib){
int i,t,u,ib2,*p,*q,lb;
lb=idxa[pos-1];
u=a[ib].n;
bucketa_merge(lb,ib);
a[lb].n+=u;
for (i=0;i<nb;++i){
ib2=idxb[i];
tu(b[ib2].hd,b[ib2].n,&t,&u);
q=b[ib2].b;
p=q+b[ib2].hd;
for (;t;++p,--t) if (pos<=*p) --*p;
for (;u;++q,--u) if (pos<=*q) --*q;
}
array_remove(idxa,pos,na);
--na;
BITa.build(a,idxa,na);
a[ib].n=ui_a;ui_a=ib;
}
void b_rebuild_rem(int pos2,int ib2){
int t,u,lb2;
lb2=idxb[pos2-1];
t=b[lb2].n;
u=b[ib2].n;
bucketb_merge(lb2,ib2);
b[lb2].n+=u;
array_remove(idxb,pos2,nb);
--nb;
b[ib2].n=ui_b;ui_b=ib2;
}
void a_movel(int pos,int ib){
int i,x,pos2,ib2,lb;
lb=idxa[pos-1];
i=a[ib].hd;
x=a[ib].a[i];
a[ib].hd=(i+1<bsz)? i+1:0;
--a[ib].n;
i=(a[lb].hd+a[lb].n)%bsz;
a[lb].a[i]=x;
++a[lb].n;
pos2=b_bsearch_pos(x,pos,&i);
ib2=idxb[pos2];
i=(b[ib2].hd+i)%bsz;
b[ib2].b[i]=pos-1;
BITa.addl(pos-1);
c.addl_b(pos2,pos-1,(ll)x);
}
void b_movel(int pos2,int ib2){
int i,pos,lb2,x;
lb2=idxb[pos2-1];
i=b[ib2].hd;
x=b[ib2].sa[i];
pos=b[ib2].b[i];
b[ib2].hd=(i+1<bsz)? i+1:0;
--b[ib2].n;
i=(b[lb2].hd+b[lb2].n)%bsz;
b[lb2].b[i]=pos;
b[lb2].sa[i]=x;
++b[lb2].n;
c.addl_a(pos2-1,pos,(ll)x);
}
ll a_sumlt(int ib,int s,int u,const int x){
int i,t,*p,*q;
ll ret=0LL;
i=(a[ib].hd+s)%bsz;
t=(bsz-i<=u)? bsz-i:u;
u-=t;
q=a[ib].a;
p=q+i;
for (i=0;i<u;++i){
if (q[i]<x) ret+=(ll)q[i];
}
for (i=0;i<t;++i){
if (p[i]<x) ret+=(ll)p[i];
}
return ret;
}
ll b_sumlt(int ib,int x,const int lpos,const int rpos){
int i,j,t,u,*p,*q;
ll ret=0LL;
u=b_bsearch2(ib,x);
j=b[ib].hd;
t=(bsz-j<=u)? bsz-j:u;
u-=t;
p=b[ib].b;
q=b[ib].sa;
for (i=0;i<u;++i){
if (lpos<p[i] && p[i]<rpos) ret+=(ll)q[i];
}
p+=j;
q+=j;
for (i=0;i<t;++i){
if (lpos<p[i] && p[i]<rpos) ret+=(ll)q[i];
}
return ret;
}
public:
BucketDecomposition(int n){
int i,*p,mbit,mseg;
bsz=(int)(sqrt((double)n))*4;
crem=0.35*bsz;
na=n/crem;
if (crem*na<n) ++na;
i=bsz*na+(bsz+bsz)*na+na*2;
mbit=BITa.calc_mem(na);
mseg=c.calc_mem(na);
a=new bucketa[na];
b=new bucketb[na];
idxa=new int [i+mbit];
pseg=new ll[mseg];
memset(idxa,0,(i+mbit)*4);
memset(pseg,0,mseg*8);
p=idxa+na;
idxb=p;p+=na;
BITa.set_mem(p);p+=mbit;
for (i=0;i<na;++i){
a[i].a=p;p+=bsz;
a[i].n=a[i].hd=0;
}
for (i=0;i<na;++i){
b[i].b=p;p+=bsz;
b[i].sa=p;p+=bsz;
b[i].n=b[i].hd=0;
}
c.set_mem(pseg);
na=nb=1;
ui_a=mi_a=ui_b=mi_b=0;
}
~BucketDecomposition(){
delete [] a;delete [] b;delete [] idxa;delete [] pseg;
}
// i番目の要素を返す
int access(int i){
int j,k=i+1;
j=BITa.bsearch(&k)-1;
j=idxa[j];
k=(a[j].hd+k-1)%bsz;
return a[j].a[k];
}
// i番目にxを追加
void insert(int i,int x){
int t,pos,ib,pos2,ib2;
t=i;
pos=BITa.bsearch(&t)-1;
BITa.add(pos+1,1);
ib=idxa[pos];
bucketa_insert(ib,t,x);
pos2=b_bsearch_pos(x,pos,&t);
ib2=idxb[pos2];
bucketb_insert(ib2,t,pos,x);
c.add(pos2,pos,(ll)x);
if (a[ib].n==bsz){
if (pos && a[idxa[pos-1]].n+1<bsz){
a_movel(pos,ib);
}else{
a_rebuild_ins(pos,ib);
c.rebuild_a_ins(b,idxb,na,nb,pos);
}
}
if (b[ib2].n==bsz){
if (pos2 && b[idxb[pos2-1]].n+1<bsz){
b_movel(pos2,ib2);
}else{
b_rebuild_ins(pos2,ib2);
c.rebuild_b_ins(b,idxb,na,nb,bsz,pos2);
}
}
//if (t) c.rebuild(b,idxb,na,nb,bsz);
}
// i番目の要素を削除
int remove(int i){
int t,pos,ib,pos2,ib2,x;
t=i+1;
pos=BITa.bsearch(&t)-1;
BITa.add(pos+1,-1);
ib=idxa[pos];
x=bucketa_remove(ib,t-1);
pos2=b_bsearch_pos(x,pos,&t);
ib2=idxb[pos2];
bucketb_remove(ib2,t);
c.add(pos2,pos,(ll)-x);
if (a[ib].n<crem){
if (pos+1<na && (a[idxa[pos+1]].n>crem)){
a_movel(pos+1,idxa[pos+1]);
}else if (pos+1<na){
a_rebuild_rem(pos+1,idxa[pos+1]);
c.rebuild_a_rem(na,nb,pos);
}else if (a[ib].n==0 && na>1){
a[ib].n=ui_a;ui_a=ib;
array_remove(idxa,pos,na);
--na;
}
}
if (b[ib2].n<crem){
if (pos2+1<nb && (b[idxb[pos2+1]].n>crem)){
b_movel(pos2+1,idxb[pos2+1]);
}else if (pos2+1<nb){
b_rebuild_rem(pos2+1,idxb[pos2+1]);
c.rebuild_b_rem(na,nb,pos2);
}else if (b[ib2].n==0 && nb>1){
b[ib2].n=ui_b;ui_b=ib2;
array_remove(idxb,pos2,nb);
--nb;
}
}
return x;
}
// [l,r)に含まれるx未満の要素の和
ll sum_less_than(int l,int r,int x){
int k,k2,lpos,rpos,pos2,ib2,lb,rb;
ll s;
k=l+1;
lpos=BITa.bsearch(&k)-1;
lb=idxa[lpos];
--k;
k2=r;
rpos=BITa.bsearch(&k2)-1;
rb=idxa[rpos];
pos2=b_bsearch(x);
ib2=idxb[pos2];
if (lpos+1<rpos){
s=a_sumlt(lb,k,a[lb].n-k,x)+a_sumlt(rb,0,k2,x)+b_sumlt(ib2,x,lpos,rpos);
if (pos2>0) s+=c.sum(pos2-1,rpos-1)-c.sum(pos2-1,lpos);
}else if (lpos<rpos){
s=a_sumlt(lb,k,a[lb].n-k,x)+a_sumlt(rb,0,k2,x);
}else{
s=a_sumlt(lb,k,r-l,x);
}
return s;
}
};
int main(void){
int i,j,t,x,l,r;
ll s;
const int N=200000; // 最大入力件数
const int M=N; // 入力値の絶対値の最大値
BucketDecomposition b(N);
// N要素追加(100回に1回は負数か0を混ぜる)
for (i=0;i<N;++i){
j= (double)rand() / (RAND_MAX + 1) * (i+1);
x= (double)rand() / (RAND_MAX + 1) * (M);
if (i%100==0) x=-x;
b.insert(j,x);
}
// sumlt10万回
for (i=0;i<100000;++i){
l= (double)rand() / (RAND_MAX + 1) * (N);
r= (double)rand() / (RAND_MAX + 1) * (N);
if (l>r) t=l,l=r,r=t;
x= (double)rand() / (RAND_MAX + 1) * (M);
s=b.sum_less_than(l,r,x);
}
// N要素削除
for (i=N;0<i;--i){
j= (double)rand() / (RAND_MAX + 1) * (i);
x=b.remove(j);
}
return 0;
}
