静的区間転倒数2
領域O(n)、構築O(n*sqrt(n/log(n)))、クエリO(sqrt(n*log(n)))、を示す
n : 数列長
B : 大バケットのサイズ
b : 小バケットのサイズ
n=Bb
大バケットがb個ある、小バケットがB個ある、b≦sqrt(n)≦B、となる
(Bとbは整数で、Bはbの整数倍とする)
(必要に応じて番兵を入れてnを調整する)
WM(ウェーブレット行列)を使うので数列は0以上n未満の整数値に変換する
順列に変換しても問題ないので変換後の数列は順列とする
I(X) : 数列の区間Xの転倒数
I(X,Y) : 集合{(x,y)|x>y,x∈X,y∈Y}の要素数
図のX' : クエリ区間で大バケットの左にはみ出した小バケットちょうど何個分の区間
図のY : クエリ区間で大バケットちょうど何個分の区間
I(X',Y)を前処理で全て計算してもn=Bbなので領域O(n)で保存できる(I(Y,Z')も同様)
図のX : クエリ区間で小バケットの左にはみ出した区間
I(X,Y)はWMで1要素ずつ計算する(I(Y,Z)も同様)
クエリ区間[L,R)が大バケットを1つ以上含む場合の転倒数
=I(XX')+I(Z'Z)+I(Y)+I(X',Y)+I(Y,Z')+I(X,Y)+I(Y,Z)+I(XX',Z'Z)
各項の時間計算量
I(XX') : 前処理によりO(1)
I(Z'Z) : 同上
I(Y) : 前処理によりO(1)
I(X',Y) : 前処理によりO(1)
I(Y,Z') : 同上
I(X,Y) : WMによりO(b*log(n))
I(Y,Z) : 同上
I(XX',Z'Z) : 大バケットを(値,添字)の値をキーにしてソートした列を使う尺取り法によりO(B)
クエリ区間が大バケットを含まず大バケット間の境界を跨ぐ場合の区間をXX'Z'Zとする(図からYを取り除いたもの)
I(XX'Z'Z)=I(XX')+I(Z'Z)+I(XX',Z'Z)
計算量はO(B)
クエリ区間が1つの大バケットに含まれる場合の区間をX(図のX)とする
I(X)=I(XX')-I(X')-I(X,X')
計算量はO(B)
以上によりクエリの計算量はO(max(b*log(n),B))となる
B=sqrt(n*log(n))、b=sqrt(n/log(n))、とすると冒頭の計算量を達成できる
Xをある大バケット、YをX以外の大(小)バケットとするとき、I(X,Y)やI(Y,X)の構築は数列の要素がn未満の整数値なので次のようにXに関する累積和みたいなものを作り、Xに含まれるx未満(以上)の要素数をO(1)で求められるようにするのがおそらく一番簡単である
サイズnの作業領域wを0で初期化
for each x in X{
w[x]+=1
}
for (i=1;i<n;++i){
w[i]+=w[i-1]
}このI(X,Y)やI(Y,X)をもとに、図のI(Y)やI(X',Y)やI(Y,Z')を構築する
変種1
n=Bbの制約を外して、B=sqrt(n)、b=sqrt(n)/log(n)、とすると領域O(n*log(n))、クエリO(sqrt(n))になる
変種2
前に「平方分割でrankltはO(sqrt(n))でできる」ことを書いたが、この平方分割をバケット法として一般化して考えると「バケット法でrankltはクエリ区間が大バケットちょうど何個分ならばO(b)でできる(bはn=Bbのb)」ということがいえる
さらに、このO(b)は小バケット内の全探索なのでWMを使うとここがO(log(b))になるので「バケット法とWMでrankltはクエリ区間が大バケットちょうど何個分ならばO(log(b))でできる」になる
これを使うと図のI(X,Y)とI(Y,Z)はO(b*log(n))からO(b*log(b))になるので、n=b*b*log(b)のとき、クエリがO(b*log(b))になる
(I(X,Y)とI(XX',Z'Z)の計算量が等しい ⇔ b*log(b)=B (B=n/b))
O(b*log(b))がnの式でどう書けるのか分からないので計算量は分からない
(下のコードでは、B=sqrt(n*log(n))、b=sqrt(n/log(n))、であるがこのバケット法を使っているのでWMの処理の部分はO(b*log(n))ではなくO(b*log(b))であり定数倍の改善になっている)
変種3
WMを使わなくてもバケット法でI(X,Y)とI(Y,Z)はO(b*b)なので、B=pow(n,2/3)、b=pow(n,1/3)、とすると領域O(n)、クエリO(pow(n,2/3))になる
#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <algorithm>
#define ll long long
#define ushort unsigned short
class RangeInversions2{
private:
struct bucket{
int *a;
ll *inv;
};
struct Bucket{
int *a,*sa,*ix,*w;
ll *inv,*li,*ri;
bucket *bk;
};
class wm16{
private:
ushort *a,bsz,nb,*tbl,db,M,nbit;
void add1(ushort *pa,ushort *ps,ushort i){
pa[i>>4]|=(ushort)1<<(i&15);
++ps[(i>>4)+1];
}
void wm_build(int *pb,int *wr,ushort *pa,ushort *ps){
int i,n0,n1,*tp,m,nt=nbit;
for (m=M;m;m>>=1){
n0=n1=nt;
for (i=0;i<nt;++i){
if (pb[i]&m){
add1(pa,ps,i);
--n0;
}
}
for (i=1;i<bsz;++i) ps[i]+=ps[i-1];
for (i=nt-1;0<=i;--i){
if (pb[i]&m) wr[--n1]=pb[i];
else wr[--n0]=pb[i];
}
tp=pb,pb=wr,wr=tp;
pa=ps+bsz;
ps=pa+bsz;
}
}
public:
int calc_mem(int nB_,int bsz_,int nb_){
int i;
nbit=bsz_;
bsz=(bsz_%16)? bsz_/16+1:bsz_/16;
++bsz;
nb=nb_*nB_;
for (i=0;(1<<i)<=nB_;++i);
M=(ushort)1<<(i-1);
db=bsz*2*i;
return db*nb+65536;
}
void init(int *pb,int bsz_,int *wr,ushort *d_){
int i,j,*p;
ushort *q;
a=d_;
tbl=a+db*nb;
p=pb;
q=a;
for (i=0;i<nb;++i){
wm_build(p,wr,q,q+bsz);
p+=bsz_;
q+=db;
}
for (i=0;i<65536;++i){
for (j=i;j;j-=j&-j) ++tbl[i];
}
}
ll rangefreqs(int *pa_,int n_,int bsz_,ushort x,ushort y,int lside){
int i,s=0;
ushort ib,t,u,m,*pa,*ps,*tb,sz,l,r,l2,r2,nt=nbit;
sz=bsz;
tb=tbl;
for (i=0;i<n_;++i){
ib=pa_[i]/bsz_;
if (lside){
r2=r=pa_[i]-ib*bsz_;
l2=l=0;
}else{
l2=l=pa_[i]-ib*bsz_+1;
r2=r=bsz_;
}
pa=a+ib*db;
ps=pa+sz;
for (m=M;m;m>>=1){
if (l2<r2){
t=tb[pa[l2>>4]&~(USHRT_MAX<<(l2&15))]+ps[l2>>4];
u=tb[pa[r2>>4]&~(USHRT_MAX<<(r2&15))]+ps[r2>>4];
if (y&m){
s+=r2-l2-(u-t);
l2=nt-ps[sz-1]+t;
r2=l2+(u-t);
}else{
l2-=t;
r2-=u;
}
}
if (l<r){
t=tb[pa[l>>4]&~(USHRT_MAX<<(l&15))]+ps[l>>4];
u=tb[pa[r>>4]&~(USHRT_MAX<<(r&15))]+ps[r>>4];
if (x&m){
s-=r-l-(u-t);
l=nt-ps[sz-1]+t;
r=l+(u-t);
}else{
l-=t;
r-=u;
}
}
pa=ps+sz;
ps=pa+sz;
}
}
return s;
}
};
int *a,N,Bsz,bsz,nB,nb,*wi,*C;
ll *Bi;
ushort *pwm;
Bucket *Bk;
bucket *bk;
wm16 wm;
void BIT_init(int *BIT,int sz){
memset(BIT,0,(sz+1)*4);
}
void BIT_add1(int *BIT,int i,int sz){
for (++i;i<=sz;i+=i&-i) ++BIT[i];
}
ll BIT_sum(int *BIT,int i){
ll s=0LL;
for (++i;0<i;i-=i&-i) s+=(ll)BIT[i];
return s;
}
void initLR_lt(int *L,int *Li,int *R,int *w){
int i,j,sz=Bsz;
memset(w,0,nb*4);
for (i=j=0;i<sz;++i){
while (j<sz && L[i]>R[j]) ++j;
w[Li[i]]+=j;
}
}
void initLR_mt(int *L,int *R,int *Ri,int *w){
int i,j,sz=Bsz;
memset(w,0,nb*4);
for (i=j=0;i<sz;++i){
while (j<sz && L[j]<=R[i]) ++j;
w[Ri[i]]+=sz-j;
}
}
void C_build(int *pb){
int j,*p,*q,*r;
p=C;
q=pb;
r=C+N;
for (;p<r;p+=bsz,q+=bsz){
for (j=0;j<bsz;++j) ++p[q[j]];
}
for (p=C+bsz;p<r;p+=bsz){
for (j=0;j<bsz;++j) p[j]+=p[j-bsz];
}
for (p=C;p<r;p+=bsz){
for (j=1;j<bsz;++j) p[j]+=p[j-1];
}
}
void Bk_build_Bi(){
int i,j;
for (i=0;i<nB;++i){
ll t=Bk[i].inv[i];
for (j=i+1;j<nB;++j) Bk[i].inv[j]=t+Bk[i].bk[0].inv[j];
}
for (i=nB-1;0<=i;--i){
for (j=i+1;j<nB;++j) Bk[i].inv[j]+=Bk[i+1].inv[j];
}
}
void Bk_build_lr(int ib,int sz,int *pa,int *pBIT){
int i;
ll *q,*r;
BIT_init(pBIT,sz);
q=Bk[ib].li;
for (i=Bsz-1;0<=i;--i){
q[i]=BIT_sum(pBIT,pa[i]-1);
BIT_add1(pBIT,pa[i],sz);
}
for (i=Bsz-1;0<i;--i) q[i-1]+=q[i];
Bk[ib].inv[ib]=q[0];
BIT_init(pBIT,sz);
r=Bk[ib].ri;
for (i=0;i<Bsz;++i){
r[i]=i-BIT_sum(pBIT,pa[i]);
BIT_add1(pBIT,pa[i],sz);
}
for (i=1;i<Bsz;++i) r[i]+=r[i-1];
}
void bk_build_Bi(){
int i,j,k;
for (i=0;i<nB;++i){
for (j=0;j<i;++j){
initLR_lt(Bk[j].sa,Bk[j].ix,Bk[i].sa,Bk[j].w);
for (k=0;k<nb;++k) Bk[j].bk[k].inv[i]=Bk[j].w[k];
}
for (j=nB-1;i<j;--j){
initLR_mt(Bk[i].sa,Bk[j].sa,Bk[j].ix,Bk[j].w);
for (k=0;k<nb;++k) Bk[j].bk[k].inv[i]=Bk[j].w[k];
}
}
for (i=0;i<nB;++i){
for (j=0;j<nb;++j){
for (k=i-1;0<k;--k) Bk[i].bk[j].inv[k-1]+=Bk[i].bk[j].inv[k];
for (k=i+1;k+1<nB;++k) Bk[i].bk[j].inv[k+1]+=Bk[i].bk[j].inv[k];
}
}
for (i=0;i<nB;++i){
for (j=1;j<nb;++j){
for (k=0;k<i;++k) Bk[i].bk[j].inv[k]+=Bk[i].bk[j-1].inv[k];
}
for (j=nb-1;0<j;--j){
for (k=i+1;k<nB;++k) Bk[i].bk[j-1].inv[k]+=Bk[i].bk[j].inv[k];
}
}
}
ll C_query(int l,int l2,int lB,int r,int r2,int rB){
int i,ib,*pa,*pc,*p,sz=bsz;
ll s;
pa=a;
pc=C;
p=pc+(nb*nB-1)*sz;
s=(ll)((p[rB-1]-p[lB])*(r2-r));
for (i=l;i<l2;++i){
ib=pa[i]/sz;
if (ib){
p=pc+(ib-1)*sz;
s+=(ll)(p[rB-1]-p[lB]);
}
}
for (i=r;i<r2;++i){
ib=pa[i]/sz;
p=pc+ib*sz;
s-=(ll)(p[rB-1]-p[lB]);
}
return s;
}
static void qs2(int *a,int *a2,int argn){
int l,r,m,m2,t;
l=0,r=argn-1;
m=a[r>>1];
m2=a2[r>>1];
while (l<=r){
while ((a[l]<m) || (a[l]==m && a2[l]<m2)) ++l;
while ((a[r]>m) || (a[r]==m && a2[r]>m2)) --r;
if (l<=r){
t=a[l],a[l]=a[r],a[r]=t;
t=a2[l],a2[l]=a2[r],a2[r]=t;
++l,--r;
}
}
if (1<l) qs2(a,a2,l);
if (l+1<argn) qs2(a+l,a2+l,argn-l);
}
ll invLR(const int l,const int r,int lB,int rB){
int i,j,k,l2,r2,u,v,*pl,*pr,*il,*ir;
ll s=0LL;
pl=Bk[lB].sa;
il=Bk[lB].ix;
pr=Bk[rB].sa;
ir=Bk[rB].ix;
l2=l-lB*Bsz;
r2=r-rB*Bsz;
u=Bsz-l2;
v=r2;
for (i=j=0;u;++i){
if (l2<=il[i]){
for (k=pl[i];v && pr[j]<k;++j){
if (ir[j]<r2) --v;
}
s+=(ll)(r2-v);
--u;
}
}
return s;
}
ll invself(int ib,int l,int r){
int i,*ix,sz=Bsz;
ll nr,s=0LL;
nr=s;
ix=Bk[ib].ix;
l%=sz;
r=(r%sz)? r%sz:sz;
for (i=0;i<sz;++i){
if (l<=ix[i]){
if (r<=ix[i]) ++nr;
else s+=nr;
}
}
return s;
}
public:
// 数列A,数列長n(1<n)
RangeInversions2(int const *A,const int n){
int i,j,k,*p,*q;
ll *pl;
const double d=log((double)n)/log(2.0);
const int mBIT=1;
nB=bsz=sqrt((double)n/d);
for (nb=d;bsz*nb*nB<n;++nb);
Bsz=bsz*nb;
N=Bsz*nB;
i=N+(Bsz*2+nb)*nB+N+(N+Bsz+mBIT);
j=((nB+Bsz*2)+(nb*nB))*nB;
k=wm.calc_mem(nB,bsz,nb);
Bk=new Bucket[nB];
bk=new bucket[nb*nB];
a=new int[i];
Bi=new ll[j];
pwm=new ushort[k];
memset(a,0,i*4);
memset(Bi,0,j*8);
memset(pwm,0,k*2);
wi=a+N;
pl=Bi;
for (i=0;i<nB;++i){
Bk[i].a=a+i*Bsz;
Bk[i].sa=wi;wi+=Bsz;
Bk[i].ix=wi;wi+=Bsz;
Bk[i].w=wi;wi+=nb;
Bk[i].inv=pl;pl+=nB;
Bk[i].li=pl;pl+=Bsz;
Bk[i].ri=pl;pl+=Bsz;
Bk[i].bk=bk+i*nb;
for (j=0;j<nb;++j){
Bk[i].bk[j].a=Bk[i].a+j*bsz;
Bk[i].bk[j].inv=pl;pl+=nB;
}
}
C=wi;wi+=N;
memcpy(a,A,n*4);
for (i=0;i<N;++i) wi[i]=i;
qs2(a,wi,n);
for (i=0;i<n;++i) a[wi[i]]=i;
for (i=n;i<N;++i) a[i]=n;
for (i=0;i<N;++i) wi[i]/=Bsz;
C_build(wi);
wm.init(wi,bsz,wi+N,pwm);
p=a;
for (i=0;i<nB;++i,p+=Bsz){
memcpy(Bk[i].sa,p,Bsz*4);
q=Bk[i].ix;
for (j=0;j<Bsz;++j) q[j]=j;
qs2(Bk[i].sa,q,Bsz);
for (j=0;j<Bsz;++j) wi[q[j]]=j;
Bk_build_lr(i,Bsz,wi,wi+Bsz);
}
for (i=0;i<nB;++i){
memcpy(wi+i*Bsz,Bk[i].ix,Bsz*4);
for (j=0;j<Bsz;++j) Bk[i].ix[j]/=bsz;
}
bk_build_Bi();
for (i=0;i<nB;++i) memcpy(Bk[i].ix,wi+i*Bsz,Bsz*4);
Bk_build_Bi();
}
~RangeInversions2(){
delete [] Bk;delete [] bk;delete [] a;delete [] Bi;delete [] pwm;
}
// Aの[l,r)の転倒数
ll inversions_query(const int l,const int r){
if (l>=r) return 0LL;
int lb,rb,lB,rB;
ll s;
lB=l/Bsz;
rB=(r-1)/Bsz;
if (lB+1<rB){
lb=l/bsz;
rb=(r-1)/bsz;
s=Bk[lB].li[l%Bsz]+Bk[rB].ri[(r-1)%Bsz]+Bk[lB+1].inv[rB-1]+invLR(l,r,lB,rB)+C_query(l,(lb+1)*bsz,lB,rb*bsz,r,rB)+wm.rangefreqs(a+l,(lb+1)*bsz-l,bsz,lB+1,rB,1)+wm.rangefreqs(a+rb*bsz,r-rb*bsz,bsz,lB+1,rB,0);
if ((lb+1)%nb) s+=bk[lb+1].inv[rB-1];
if (rb%nb) s+=bk[rb-1].inv[lB+1];
}else if (lB<rB){
s=Bk[lB].li[l%Bsz]+Bk[rB].ri[(r-1)%Bsz]+invLR(l,r,lB,rB);
}else{
s=(r%Bsz)? Bk[lB].li[l%Bsz]-Bk[lB].li[r%Bsz]-invself(lB,l,r):Bk[lB].li[l%Bsz];
}
return s;
}
};
int main(void){
int i,t,l,r;
ll s;
const int n=1000000; // 数列長
int *a=new int[n];
// 数列を作成
for (i=0;i<n;++i){
a[i]=(double)rand() / (RAND_MAX + 1) * (INT_MAX);
}
// 構築
RangeInversions2 ri(a,n);
// クエリ1万回
for (i=0;i<10000;++i){
l=(double)rand() / (RAND_MAX + 1) * (n);
r=(double)rand() / (RAND_MAX + 1) * (n);
if (l>r) t=l,l=r,r=t;
s=ri.inversions_query(l,r);
}
delete [] a;
return 0;
}
動的WMのB木による高速化の可能性について
WM(ウェーブレット行列)は上位ビットからの基数ソートのようなものであること、基数ソートを1ビット単位で行うより8ビット単位で行うほうが速そうであること、以上のことから動的WMも8ビット単位で行うと速くなるかもしれないので試すことにした
次の処理ができればよい
insert(i,x) : 数列のi番目にxを追加(xは0以上255以下の整数)
remove(i) : 数列のi番目の要素を削除
frequency(l,r) : 数列の区間[l,r)に含まれる0から255までの値の出現回数
これらはnを数列長としたとき平衡二分木でO(log(n))でできる(各ノードにc[256]を持たせて、平衡処理においてそのノードを根とした木に値xがc[x]個含まれているという状態を保つようにする)
しかし、この方法だとfrequencyのO(log(n))は平衡二分木なので定数はlogの底2で256*log(n)なので遅い
多分木であるB木を使って木の高さを抑える。B木はM分木とする。(オリジナルのB木は順序タイプであり、ここで使うのは列タイプなのでやや異なるがここではこれもB木とする)
内部ノード
node{
n // 子ノードの個数
index[M] // 子ノードのインデックス
size[M] // size[i] : index[i]ノードを根とした木に含まれる要素数
c[M][256] // c[i][x] : 0-th から i-th children までに含まれるxの個数
}
葉ノード
leafNode{
n // 数列長
a[leafsize] // 数列
}B木のinsert処理の擬似コードを示す
数列のi番目にxを追加する
Iはノードのインデックスで初期値はB木の根
Btree_insert(I,i,x){
if (内部ノード){
for (j=0;node[I].size[j]<i;++j) i-=node[I].size[j]
node[I].size[j]+=1 // 追加先はindex[j]ノード
for (k=j;k<node[I].n;++k) node[I].c[k][x]+=1 // cの更新(注1)
ret=Btree_insert(node[I].index[j],i,x)
if (ret){
子index[j]が満杯なのでsplit処理(子ノードが1つ増える)
}
return (node[I].n==M)
}else{
// 葉ノード
leafNode[I].a[i]にxを追加
leafNode[I].n+=1
return (leafNode[I].n==leafsize)
}
} 注1の1行(メモリ上の等間隔に並んだ複数の場所を+1する処理)が並列処理などによって高速化できるとM分木のMを大きくできて木の高さが低くなるので全処理の高速化につながる(removeの場合は-1する処理)
下のコードだとinsert2、remove2内の次の部分
for (t=b[ib].n-pos;t;q+=256,--t) ++*q;
for (t=b[ib].n-pos;t;q+=256,--t) --*q;
上記よりは優先順位は低いと考えられるがvoid freq内の次の部分(複数の配列の足し算)も並列処理の対象となる
for (i=0;i<256;++i) a[i]=b[i]+c[i]+...
動的WMのinsertとquantileの擬似コードを示す
8ビット単位での処理なので、32ビットデータを扱うならばB木を4本使うことになる
上位ビットから下位ビットに進んでいく
BはB木
B.ranklt(x)は数列全体に含まれるx未満の要素の個数
B.rank(i,x)は数列の[0,i)に含まれるxの個数
sは注目する8ビットを(x>>s)&255で取り出すためのシフト量で、32ビットデータを扱うならばsの初期値は32-8=24になる
// 数列のi番目にxを追加
WM_insert(B,i,x,s){
y=(x>>s)&255
B.insert(i,y)
if (0<s){
j=B.ranklt(y)
k=B.rank(i,y)
WM_insert(Bの下位の木,j+k,s-8)
}
}
// [l,r)のk番目の要素
WM_quantile(B,l,r,k,s){
f=B.frequency(l,r) // fはサイズ256の配列
for (x=0;f[x]<=k;++x) k-=f[x] // k番目の要素はx
if (0<s){
i=B.ranklt(x)
L=B.rank(l,x)
R=B.rank(r,x)
ret=quantil(Bの下位の木,i+L,i+R,k,s-8)
return (x<<s)+ret
}else{
return x
}
}
結果
100万要素の32ビットデータをinsert後、各クエリ100万回、100万要素removeで計測(removeの時間を1とし右に行くほど遅い)
remove(1),insert(1.05),rank(1.31),ranklt(1.95),select(1.96),quantile(2.62)
ビットベクトル版WMの速度は処理間の定数倍の差が小さいので上の並びのどこか狭い範囲に対応するはず
下のB木のコードで動的WMを作るための部品は全てそろっているので書いて動かすとおおよその感じはつかめるかもしれない
(WMのコードは計算量の定数倍改善をどの程度行うかなどについて決められていないため未完成)
その他
B木を使うのでB木の持つ利点を生かせる(入力件数が大きい場合に有利、データベースへの応用など)
メモリ使用量が元データの10倍近くになる
8ビット単位での処理なので補助データを使う場合のメモリ使用量がビットベクトル版の1/8程になる
#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <algorithm>
#define uchar unsigned char
class Btree256{
private:
static const int M=32;
static const int leafsize=1024;
static const int crem=12;
static const int crem_l=409;
struct node{
int n,ip,pos,a[M],b[M],c[M*256];
};
struct leaf_node{
int n;
uchar a[leafsize];
};
node *b;
leaf_node *bl;
int memi[16],wi[256],hi,ui,mi,ui_l,mi_l,root,*ps[10],*qs[10];
uchar Tx;
int getindex(int internalnode){
int i;
if (internalnode){
if (ui!=-1) i=ui,ui=b[i].n;
else i=++mi;
b[i].n=0;
}else{
if (ui_l!=-1) i=ui_l,ui_l=bl[i].n;
else i=++mi_l;
bl[i].n=0;
}
return i;
}
void putindex(int ib,int internalnode){
if (internalnode){
b[ib].n=ui;ui=ib;
}else{
bl[ib].n=ui_l;ui_l=ib;
}
}
void node_popl_shiftl(int ib,int l,int r){
if (l<r){
memmove(b[ib].a+l,b[ib].a+(l+1),(r-l)*4);
memmove(b[ib].b+l,b[ib].b+(l+1),(r-l)*4);
}
}
void node_popr_shiftr(int ib,int l,int r){
if (l<r){
memmove(b[ib].a+(l+1),b[ib].a+l,(r-l)*4);
memmove(b[ib].b+(l+1),b[ib].b+l,(r-l)*4);
}
}
int finsert(int ib,int k,int xa,int xb){
node_popr_shiftr(ib,k,b[ib].n);
b[ib].a[k]=xa;
b[ib].b[k]=xb;
return (++b[ib].n==M);
}
int leaf_insert(int ib,int k,uchar x){
memmove(bl[ib].a+(k+1),bl[ib].a+k,bl[ib].n-k);
bl[ib].a[k]=x;
return (++bl[ib].n==leafsize);
}
int fremove(int ib,int k){
node_popl_shiftl(ib,k,b[ib].n);
return (--b[ib].n<crem);
}
int leaf_remove(int ib,int k,uchar *x){
*x=bl[ib].a[k];
memmove(bl[ib].a+k,bl[ib].a+(k+1),--bl[ib].n-k);
return (bl[ib].n<crem_l);
}
void underflow(int ib,int pos,int internalnode){
int ic,rc,nl,nr,t,*p,*q,*r;
uchar x;
ic=b[ib].a[pos];
rc=b[ib].a[pos+1];
if (internalnode-1){
nl=b[ic].n;
nr=b[rc].n;
t=b[rc].b[0];
finsert(ic,nl,b[rc].a[0],t);
b[ib].b[pos]+=t;
b[ib].b[pos+1]-=t;
q=b[ic].c+nl*256;
r=b[rc].c;
memcpy(q,r,256*4);
c_add_c(q,q-256);
p=b[ib].c+pos*256;
c_add_c(p,r);
cs_sub_c(rc,0,nr-1);
c_popl_shiftl(rc,0,nr-1);
fremove(rc,0);
}else{
nl=bl[ic].n;
leaf_remove(rc,0,&x);
leaf_insert(ic,nl,x);
++b[ib].b[pos];
--b[ib].b[pos+1];
p=b[ib].c+pos*256;
++p[x];
}
}
int array_sum(int *p,int n_){
int s=0;
for (int i=0;i<n_;++i) s+=p[i];
return s;
}
int node_split(int ib,int newi,int internalnode){
if (internalnode){
int t=M/2;
memcpy(b[newi].a,b[ib].a+t,t*4);
memcpy(b[newi].b,b[ib].b+t,t*4);
b[ib].n=b[newi].n=t;
split_c(ib,newi);
return array_sum(b[newi].b,t);
}else{
int t=leafsize/2;
memcpy(bl[newi].a,bl[ib].a+t,t);
bl[ib].n=bl[newi].n=t;
return t;
}
}
int node_merge(int lb,int ib,int internalnode){
if (internalnode){
int t=b[lb].n;
int u=b[ib].n;
memcpy(b[lb].a+t,b[ib].a,u*4);
memcpy(b[lb].b+t,b[ib].b,u*4);
merge_c(lb,ib);
b[lb].n+=u;
return array_sum(b[lb].b,t+u);
}else{
int t=bl[lb].n;
int u=bl[ib].n;
memcpy(bl[lb].a+t,bl[ib].a,u);
bl[lb].n+=u;
return t+u;
}
}
void c_popr_shiftr(int ib,int l,int r){
if (l<r){
int *p=b[ib].c+l*256;
memmove(p+256,p,(r-l)*256*4);
}
}
void c_popl_shiftl(int ib,int l,int r){
if (l<r){
int *p=b[ib].c+l*256;
memmove(p,p+256,(r-l)*256*4);
}
}
void c_count_leaf(int ib,int n_,int *ans,int init=0){
uchar *p=bl[ib].a;
if (init) memset(ans,0,256*4);
for (int i=0;i<n_;++i) ++ans[p[i]];
}
void c_add_c(int *c,int *c2){
for (int i=0;i<256;++i) c[i]+=c2[i];
}
void c_sub_c(int *c,int *c2){
for (int i=0;i<256;++i) c[i]-=c2[i];
}
void cs_add_c(int ib,int i,int n_){
int *p=b[ib].c+i*256;
for (int *q=p+256;0<n_;--n_,q+=256){
for (int k=0;k<256;++k) q[k]+=p[k];
}
}
void cs_sub_c(int ib,int i,int n_){
int *p=b[ib].c+i*256;
for (int *q=p+256;0<n_;--n_,q+=256){
for (int k=0;k<256;++k) q[k]-=p[k];
}
}
void split_c(int ib,int rb){
cs_sub_c(ib,M/2-1,M/2);
memcpy(b[rb].c,b[ib].c+(M/2)*256,(M/2)*256*4);
}
void merge_c(int ib,int rb){
int nl=b[ib].n;
int nr=b[rb].n;
memcpy(b[ib].c+nl*256,b[rb].c,nr*256*4);
cs_add_c(ib,nl-1,nr);
}
void insert_c(int ib,int pos,int internalnode){
int ic,*p,*q;
ic=b[ib].a[pos];
c_popr_shiftr(ib,pos,b[ib].n);
p=b[ib].c+pos*256;
if (internalnode-1){
q=b[ic].c+((M/2-1)*256);
memcpy(p,q,256*4);
if (pos) c_add_c(p,p-256);
}else{
if (pos) memcpy(p,p-256,256*4);
c_count_leaf(ic,leafsize/2,p,!pos);
}
}
void remove_c(int ib,int pos){
int *p=b[ib].c+pos*256;
memmove(p,p+256,(b[ib].n-(pos+1))*256*4);
}
int insert2(int ib,int k,uchar x,int ip){
int t,pos,ic,*p,*q,h,ret;
for (h=hi;1<h;--h){
b[ib].ip=ip;
p=b[ib].b;
for (pos=0;p[pos]<k;++pos) k-=p[pos];
++p[pos];
b[ib].pos=pos;
ic=b[ib].a[pos];
q=b[ib].c+(pos*256+x);
for (t=b[ib].n-pos;t;q+=256,--t) ++*q;
ip=ib;
ib=ic;
}
memmove(bl[ib].a+(k+1),bl[ib].a+k,bl[ib].n-k);
bl[ib].a[k]=x;
ret=(++bl[ib].n==leafsize);
for (ic=ib,ib=ip;ret && ib!=-1;++h){
int newi=getindex(h-1);
pos=b[ib].pos;
insert_c(ib,pos,h);
t=node_split(ic,newi,h-1);
b[ib].b[pos]-=t;
ret=finsert(ib,pos+1,newi,t);
ic=ib;
ib=b[ib].ip;
}
return ret;
}
int remove2(int ib,int k,int ip){
int t,pos,ic,*p,*q,rc,h,ret;
uchar x;
for (h=hi;1<h;--h){
b[ib].ip=ip;
p=b[ib].b;
for (pos=0;k>=p[pos];++pos) k-=p[pos];
--p[pos];
b[ib].pos=pos;
ip=ib;
ib=b[ib].a[pos];
}
Tx=x=bl[ib].a[k];
memmove(bl[ib].a+k,bl[ib].a+(k+1),--bl[ib].n-k);
ret=(bl[ib].n<crem_l);
for (ic=ib,ib=ip;ib!=-1;++h){
pos=b[ib].pos;
q=b[ib].c+(pos*256+x);
for (t=b[ib].n-pos;t;q+=256,--t) --*q;
if (ret){
if (pos+1<b[ib].n){
rc=b[ib].a[pos+1];
ret=(h-1)? b[rc].n<=crem:bl[rc].n<=crem_l;
if (!ret) underflow(ib,pos,h);
}else{
rc=ic;
ic=b[ib].a[--pos];
ret=(h-1)? (b[ic].n+b[rc].n<crem*2 || !b[rc].n):(bl[ic].n+bl[rc].n<crem_l*2 || !bl[rc].n);
}
if (ret){
remove_c(ib,pos);
b[ib].b[pos]=node_merge(ic,rc,h-1);
putindex(rc,h-1);
ret=fremove(ib,pos+1);
}
}
ic=ib;
ib=b[ib].ip;
}
return ret;
}
void freq(int **P,int **Q,int pq,int *ans){
int i,*p0,*p1,*p2,*p3,*p4,*p5,*p6,*q0,*q1,*q2,*q3,*q4,*q5,*q6;
switch (pq){
case 1:
q0=Q[0];
for (i=0;i<256;++i) ans[i]=-q0[i];break;
case 10:
p0=P[0];
for (i=0;i<256;++i) ans[i]=p0[i];break;
case 11:
p0=P[0],q0=Q[0];
for (i=0;i<256;++i) ans[i]=p0[i]-(q0[i]);break;
case 12:
p0=P[0],q0=Q[0],q1=Q[1];
for (i=0;i<256;++i) ans[i]=p0[i]-(q0[i]+q1[i]);break;
case 13:
p0=P[0],q0=Q[0],q1=Q[1],q2=Q[2];
for (i=0;i<256;++i) ans[i]=p0[i]-(q0[i]+q1[i]+q2[i]);break;
case 14:
p0=P[0],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3];
for (i=0;i<256;++i) ans[i]=p0[i]-(q0[i]+q1[i]+q2[i]+q3[i]);break;
case 15:
p0=P[0],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4];
for (i=0;i<256;++i) ans[i]=p0[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]);break;
case 16:
p0=P[0],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5];
for (i=0;i<256;++i) ans[i]=p0[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]);break;
case 17:
p0=P[0],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5],q6=Q[6];
for (i=0;i<256;++i) ans[i]=p0[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]+q6[i]);break;
case 20:
p0=P[0],p1=P[1];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i];break;
case 21:
p0=P[0],p1=P[1],q0=Q[0];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]-(q0[i]);break;
case 22:
p0=P[0],p1=P[1],q0=Q[0],q1=Q[1];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]-(q0[i]+q1[i]);break;
case 23:
p0=P[0],p1=P[1],q0=Q[0],q1=Q[1],q2=Q[2];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]-(q0[i]+q1[i]+q2[i]);break;
case 24:
p0=P[0],p1=P[1],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]-(q0[i]+q1[i]+q2[i]+q3[i]);break;
case 25:
p0=P[0],p1=P[1],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]);break;
case 26:
p0=P[0],p1=P[1],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]);break;
case 27:
p0=P[0],p1=P[1],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5],q6=Q[6];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]+q6[i]);break;
case 30:
p0=P[0],p1=P[1],p2=P[2];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i];break;
case 31:
p0=P[0],p1=P[1],p2=P[2],q0=Q[0];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]-(q0[i]);break;
case 32:
p0=P[0],p1=P[1],p2=P[2],q0=Q[0],q1=Q[1];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]-(q0[i]+q1[i]);break;
case 33:
p0=P[0],p1=P[1],p2=P[2],q0=Q[0],q1=Q[1],q2=Q[2];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]-(q0[i]+q1[i]+q2[i]);break;
case 34:
p0=P[0],p1=P[1],p2=P[2],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]-(q0[i]+q1[i]+q2[i]+q3[i]);break;
case 35:
p0=P[0],p1=P[1],p2=P[2],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]);break;
case 36:
p0=P[0],p1=P[1],p2=P[2],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]);break;
case 37:
p0=P[0],p1=P[1],p2=P[2],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5],q6=Q[6];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]+q6[i]);break;
case 40:
p0=P[0],p1=P[1],p2=P[2],p3=P[3];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i];break;
case 41:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],q0=Q[0];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]-(q0[i]);break;
case 42:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],q0=Q[0],q1=Q[1];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]-(q0[i]+q1[i]);break;
case 43:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],q0=Q[0],q1=Q[1],q2=Q[2];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]-(q0[i]+q1[i]+q2[i]);break;
case 44:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]-(q0[i]+q1[i]+q2[i]+q3[i]);break;
case 45:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]);break;
case 46:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]);break;
case 47:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5],q6=Q[6];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]+q6[i]);break;
case 50:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i];break;
case 51:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],q0=Q[0];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]-(q0[i]);break;
case 52:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],q0=Q[0],q1=Q[1];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]-(q0[i]+q1[i]);break;
case 53:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],q0=Q[0],q1=Q[1],q2=Q[2];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]-(q0[i]+q1[i]+q2[i]);break;
case 54:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]-(q0[i]+q1[i]+q2[i]+q3[i]);break;
case 55:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]);break;
case 56:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]);break;
case 57:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5],q6=Q[6];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]+q6[i]);break;
case 60:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i];break;
case 61:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],q0=Q[0];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]-(q0[i]);break;
case 62:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],q0=Q[0],q1=Q[1];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]-(q0[i]+q1[i]);break;
case 63:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],q0=Q[0],q1=Q[1],q2=Q[2];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]-(q0[i]+q1[i]+q2[i]);break;
case 64:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]-(q0[i]+q1[i]+q2[i]+q3[i]);break;
case 65:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]);break;
case 66:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]);break;
case 67:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5],q6=Q[6];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]+q6[i]);break;
case 70:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i];break;
case 71:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6],q0=Q[0];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i]-(q0[i]);break;
case 72:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6],q0=Q[0],q1=Q[1];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i]-(q0[i]+q1[i]);break;
case 73:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6],q0=Q[0],q1=Q[1],q2=Q[2];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i]-(q0[i]+q1[i]+q2[i]);break;
case 74:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i]-(q0[i]+q1[i]+q2[i]+q3[i]);break;
case 75:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]);break;
case 76:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]);break;
case 77:
p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6],q0=Q[0],q1=Q[1],q2=Q[2],q3=Q[3],q4=Q[4],q5=Q[5],q6=Q[6];
for (i=0;i<256;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i]-(q0[i]+q1[i]+q2[i]+q3[i]+q4[i]+q5[i]+q6[i]);break;
}
}
void ini2(int il,int n_,int h_){
int i,j,k,ib,ul=mi;
for (i=n_;i;i-=k){
ib=mi++;
k=(M/2<=i)? M/2:i;
b[ib].n=k;
for (j=0;j<k;++j,++il){
b[ib].a[j]=il;
b[ib].b[j]=(1<h_)? array_sum(b[il].b,b[il].n):bl[il].n;
}
}
if (ul+1<mi && b[ib].n<crem){
memcpy(b[ib-1].a+b[ib-1].n,b[ib].a,b[ib].n*4);
memcpy(b[ib-1].b+b[ib-1].n,b[ib].b,b[ib].n*4);
b[ib-1].n+=b[ib].n;
--mi;
}
if (ul+1==mi){
root=--mi;
hi=h_+1;
}else{
ini2(ul,mi-ul,h_+1);
}
}
void ini3(int ib,int *pc,int internalnode){
int i,t,*p,*q;
if (internalnode){
t=b[ib].n;
if (internalnode==1) memset(b[ib].c,0,t*256*4);
p=b[ib].a;
q=b[ib].c;
for (i=0;i<t;++i){
ini3(p[i],q+256*i,internalnode-1);
}
for (i=1;i<t;++i,q+=256){
c_add_c(q+256,q);
}
if (pc!=NULL){
memcpy(pc,q,256*4);
}
}else{
c_count_leaf(ib,bl[ib].n,pc);
}
}
public:
// n:最大入力件数
Btree256(int n){
int i,ni,nl;
double c=0.4;
// double h=(log((double)n)-log((double)leafsize)+log((double)M))/(log(c)+log((double)M));
nl=(double)n/(c*leafsize);
ni=((double)nl-1)/(c*M-1);
++nl,++ni;
// crem=c*M;
// crem_l=c*leafsize;
b=new node[ni];
bl=new leaf_node[nl];
for (i=0;i<16;++i) memi[i]=0;
for (i=0;i<nl;++i) bl[i].n=0;
ui=mi=ui_l=-1;
mi_l=root=0;
hi=1;
}
~Btree256(){
delete [] b;delete [] bl;
}
// 初期化時に数列の最後尾にxを追加する
void init(uchar x){
++memi[x>>4];
if (bl[mi_l].n==(leafsize/2)) ++mi_l;
bl[mi_l].a[bl[mi_l].n++]=x;
}
// init後にこれを実行すると初期化完了
void init_end(){
if (mi_l && bl[mi_l].n<crem_l){
memcpy(bl[mi_l-1].a+bl[mi_l-1].n,bl[mi_l].a,bl[mi_l].n);
bl[mi_l-1].n+=bl[mi_l].n;
--mi_l;
}
if (mi_l==0) return;
mi=0;
ini2(0,mi_l+1,1);
ini3(root,NULL,hi-1);
}
// i番目の要素を返す
uchar access(int i){
int *p,pos,ib=root;
for (int h=hi;1<h;--h){
p=b[ib].b;
for (pos=0;i>=p[pos];++pos) i-=p[pos];
ib=b[ib].a[pos];
}
return bl[ib].a[i];
}
// i番目にxを追加
void insert(int i,uchar x){
++memi[x>>4];
if (insert2(root,i,x,-1)){
int ro=getindex(1);
int r=getindex(hi-1);
b[ro].a[0]=root;
b[ro].a[1]=r;
if (1<hi){
memcpy(b[ro].c,b[root].c+((M/2-1)*256),256*4);
memcpy(b[ro].c+256,b[root].c+((M-1)*256),256*4);
b[ro].b[0]=array_sum(b[root].b,M/2);
b[ro].b[1]=node_split(root,r,1);
}else{
node_split(root,r,0);
c_count_leaf(root,leafsize/2,b[ro].c,1);
memcpy(b[ro].c+256,b[ro].c,256*4);
c_count_leaf(r,leafsize/2,b[ro].c+256);
b[ro].b[0]=b[ro].b[1]=leafsize/2;
}
b[ro].n=2;
root=ro;
++hi;
}
}
// i番目の要素を削除
void remove(int i){
if (remove2(root,i,-1) && 1<hi && b[root].n==1){
int i=b[root].a[0];
putindex(root,1);
root=i;
--hi;
}
--memi[Tx>>4];
}
// [l,r)に含まれる0から255までの値の出現回数をans配列に返す
void frequency(int l,int r,int *ans){
int i,pos,pos2,*p,*q,ib,ib2,np,nq,d,d2;
ib=ib2=root;
np=nq=d=d2=0;
for (i=hi;1<i;--i){
q=b[ib2].b;
for (pos2=0;q[pos2]<l;++pos2) l-=q[pos2];
if (pos2) qs[nq++]=b[ib2].c+(pos2-1)*256;
else qs[nq]=b[ib2].c;
ib2=b[ib2].a[pos2];
p=b[ib].b;
for (pos=0;p[pos]<r;++pos) r-=p[pos];
if (pos) ps[np++]=b[ib].c+(pos-1)*256;
ps[np]=b[ib].c;
ib=b[ib].a[pos];
if (ib==ib2 && pos) --np,--nq;
}
if (np){
if (q[pos2]<l+l){
if (pos2) qs[nq-1]+=256;
else ++nq;
d2=1;
}
if (p[pos]<r+r){
if (pos) ps[np-1]+=256;
else ++np;
d=1;
}
freq(ps,qs,np*10+nq,ans);
if (d==0){
for (uchar *u=bl[ib].a;r;++u,--r) ++ans[*u];
}else{
i=bl[ib].n-r;
for (uchar *u=bl[ib].a+r;i;++u,--i) --ans[*u];
}
if (d2==0){
for (uchar *u=bl[ib2].a;l;++u,--l) --ans[*u];
}else{
i=bl[ib2].n-l;
for (uchar *u=bl[ib2].a+l;i;++u,--i) ++ans[*u];
}
}else{
memset(ans,0,1024);
uchar *u=bl[ib].a;
for (i=l;i<r;++i) ++ans[u[i]];
}
}
// 数列全体に含まれるx未満の要素の個数
int ranklt(uchar x){
int i,s=0;
if (1<hi){
int *p=memi;
for (i=x>>4;i;s+=p[--i]);
p=b[root].c+((b[root].n-1)*256+(x&0xf0));
for (i=x&0x0f;i;s+=p[--i]);
}else{
uchar *q=bl[root].a;
for (i=bl[root].n;i;) if (q[--i]<x) ++s;
}
return s;
}
// [0,r)に含まれるxの個数
int rank(int r,uchar x){
int i,pos,*p,s,ib=root;
pos=-1;
s=0;
for (i=hi;1<i;--i){
p=b[ib].b;
for (pos=0;p[pos]<r;++pos) r-=p[pos];
p=b[ib].c+(pos*256+x);
if (pos) s+=p[-256];
ib=b[ib].a[pos];
}
if (r+r<=bl[ib].n || pos==-1){
for (uchar *q=bl[ib].a;r;++q,--r){
if (*q==x) ++s;
}
}else{
s+=pos? *p-p[-256]:*p;
i=bl[ib].n-r;
for (uchar *q=bl[ib].a+r;i;++q,--i){
if (*q==x) --s;
}
}
return s;
}
// i番目のxの位置(存在しない場合は-1を返す)
int select(int i,uchar x){
int j,pos,*p,r,s,ib=root;
uchar *q;
s=0;
++i;
for (int h=hi;1<h;--h){
p=b[ib].c+x;
r=b[ib].n;
for (pos=0;pos<r && *p<i;++pos) p+=256;
if (r<=pos) return -1;
if (pos){
i-=p[-256];
p=b[ib].b;
for (j=0;j<pos;++j) s+=p[j];
}
ib=b[ib].a[pos];
}
q=bl[ib].a;
for (j=0;1;++j){
if (q[j]==x && --i==0) return s+j;
}
}
};
int main(void){
int i,j,t,l,r,ans[256];
uchar x;
const int N=10000000; // 最大入力件数
const int n=1000000; // 削除及び入力件数
Btree256 b(N);
// N要素で初期化
for (i=0;i<N;++i){
x= (double)rand() / (RAND_MAX + 1) * (256);
b.init(x);
}
b.init_end();
// n要素削除
for (i=0;i<n;++i){
j= (double)rand() / (RAND_MAX + 1) * (N-i);
b.remove(j);
}
// n要素追加
for (i=n;0<i;--i){
j= (double)rand() / (RAND_MAX + 1) * (N-i);
x= (double)rand() / (RAND_MAX + 1) * (256);
b.insert(j,x);
}
// frequency10万回
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;
b.frequency(l,r,ans);
}
// 最後のクエリの結果を表示
std::cout <<"frequency(l,r)"<< std::endl;
std::cout <<"l=" <<l<< ",r="<<r << std::endl;
std::cout <<std::endl;
std::cout <<"i,frequency(i)"<< std::endl;
for (i=t=0;i<256;++i){
t+=ans[i];
std::cout <<i <<","<< ans[i] << std::endl;
}
std::cout <<std::endl;
std::cout <<"sigma frequency(i)=" << t << std::endl;
std::cout <<"r-l=" << (r-l) << std::endl;
return 0;
}
静的区間転倒数
https://suisen-kyopro.hatenablog.com/entry/2022/10/15/010252
が面白かったです。読んで考えたことについて書きます。
0<=x<=1,bs=pow(n,x)として転倒数の計算量について次を示す(pow(x,y)=xのy乗、bs=バケットサイズ、n=数列長)
領域 O(max(pow(n,2(1-x)),pow(n,1+x)))
構築 O(max(pow(n,2-x),pow(n,1+x)))
クエリ O(bs)
タイプ1(クエリ区間が2つ以上のブロック(バケット)に跨るケース)は次のように計算する
I(XYZ) = I(X) + I(Y) + I(Z) + I(X,Y) + I(Y,Z) + I(X,Z)
I(X),I(Y),I(Z)は前計算で求めておく
I(X,Z)は尺取り法(バケットの要素を添字と一緒にソートしておく)
I(X,Y),I(Y,Z)は説明する
要素列(1次元)とバケット法(2次元)の表記を小文字と大文字で区別する
aやAなどのサイズの宣言は、a[n]、A[n/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]] )
ここまでは前回とほぼ同じでさらに
SA[i]とB[i]バケットをB[i]の要素をキーにしてソートしたものを新たにSA[i]とB[i]とする
これに合わせてd列も変更する
x=a[i],x∈A[I],x∈SA[J]のとき、SA[J]内でxからxの左(右)側の任意要素までにxより大きい(小さい)要素がいくつあるかを累積和のテーブルを作りO(1)で取得できるようにする(A[I](xと同一バケット)に含まれる要素についてはカウントしないことに注意)
テーブルのxの行をtb[i]とする(iはa[i]のi)
列成分をtb[i][0],tb[i][1],...とする
下の表はn=16の平方分割の例(テーブルは1個しか書いてない)
I(X,Y)を求める
XはA[L-1]内の連続部分列、YはA[L],A[L+1],...,A[R-1]バケットからなる連続部分列とする
x∈X について I({x},Y) を1要素ずつ求めていく
x=a[i],x∈SA[J]とする
t = Yの要素のうちSA[0]からSA[J-1]までに含まれる要素の個数 = C[R-1][J-1] - C[L-1][J-1]
u = A[0]からA[R-1]までの要素のうちSA[J]に含まれる要素の個数 = C[R-1][J] - C[R-1][J-1]
I({x},Y) = t + tb[i][u-1]
I(Y,Z)も同様である
タイプ2(クエリ区間が1つのブロックに含まれるケース)は前計算したブロック内転倒数などを使うとできる
領域計算量のpow(n,2(1-x))はバケット間転倒数とCの計算量、pow(n,1+x)はテーブルの計算量
構築計算量のpow(n,2-x)はバケット間転倒数の計算量、pow(n,1+x)はテーブルの計算量
メモ
I(XYZ) = I(XY) + I(YZ) - I(Y) + I(X,Z)
という冒頭の記事にある等式はブロック境界と無関係に成り立つので、連続部分列の転倒数を求めるのと同じ計算量でI(X,Z)という2つの部分列の転倒数みたいなものが求められる
テーブルの累積和は+1ずつしか増えないのでビットベクトルで書くと省メモリになる
バケット間転倒数とCは昇順列の集まりなのでデータ圧縮できるかもしれない
領域O(n)、クエリO(sqrt(n))は領域とクエリのバランスがよいので作りたかったけれどできなかった
下のコードについて
bsがpow(n,1/3)のとき領域が最小になるのでそれで書いてある
#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <algorithm>
#define ll long long
template <typename T>
class RangeInversions{
private:
short **tbl;
int *A,bsz,nb,*b,**C,*d,N,**Bl,**Br;
ll **B,*pB;
template <typename U>
void qs2(U *a,int *a2,int argn){
U m,u;
int l,r,m2,t;
l=0,r=argn-1;
m=a[r>>1];
m2=a2[r>>1];
while (l<=r){
while ((a[l]<m) || (a[l]==m && a2[l]<m2)) ++l;
while ((a[r]>m) || (a[r]==m && a2[r]>m2)) --r;
if (l<=r){
u=a[l],a[l]=a[r],a[r]=u;
t=a2[l],a2[l]=a2[r],a2[r]=t;
++l,--r;
}
}
if (1<l) qs2(a,a2,l);
if (l+1<argn) qs2(a+l,a2+l,argn-l);
}
ll initLR(int *L,int *R,int sz){
int i,j,s;
for (i=j=s=0;i<sz;++i){
while (j<sz && L[i]>R[j]) ++j;
s+=j;
}
return (ll)s;
}
void C_build(int *pb){
int i,j;
for (i=0;i<N;++i) ++C[pb[i]/bsz][i/bsz];
for (int i=1;i<nb;++i){
for (j=0;j<nb;++j) C[i][j]+=C[i-1][j];
}
for (i=0;i<nb;++i){
for (j=1;j<nb;++j) C[i][j]+=C[i][j-1];
}
}
void B_build_self(int *p,int ib){
int i,j,s,*q;
q=Bl[ib];
for (i=bsz-2;0<=i;--i){
s=0;
for (j=i+1;j<bsz;++j){
if (p[i]>p[j]) ++s;
}
q[i]=q[i+1]+s;
}
B[ib][ib]=(ll)q[0];
q=Br[ib];
for (i=1;i<bsz;++i){
s=0;
for (j=0;j<i;++j){
if (p[j]>p[i]) ++s;
}
q[i]=q[i-1]+s;
}
}
void B_build(int *p){
int i,j,*q;
for (i=0;i<nb;++i,p+=bsz){
q=p+bsz;
for (j=i+1;j<nb;++j,q+=bsz) B[i][j]=initLR(p,q,bsz);
}
for (i=0;i<nb;++i){
for (j=i+1;j<nb;++j) B[i][j]+=B[i][j-1];
}
for (i=nb-1;0<=i;--i){
for (j=i+1;j<nb;++j) B[i][j]+=B[i+1][j];
}
}
void table_build2(short *tb,int *sa,int *ix,int x,int ib,const int i2){
int j,pos=ib*bsz;
short t=0;
for (j=i2%bsz;0<=j;--j) tb[j]=(x<sa[j] && ix[j]<pos)? ++t:t;
pos+=bsz;
t=0;
for (j=(i2%bsz)+1;j<bsz;++j) tb[j]=(x>sa[j] && ix[j]>=pos)? ++t:t;
}
void table_build(int *sa,int *ix,int *pd){
int i,j,k,x,t;
short *tb=tbl[0];
for (i=k=0;i<nb;++i){
for (j=0;j<bsz;++j,++k,tb+=bsz){
x=A[k];
t=(x/bsz)*bsz;
table_build2(tb,sa+t,ix+t,x,i,pd[k]);
}
}
}
ll invLM(int l,int lb,int rb){
int pos,*pl,*pr,*pd,sz=bsz;
short *tb=tbl[0]+l*sz;
ll s=0LL;
pl=C[lb-1];
pr=C[rb];
pd=d+l;
for (l=lb*sz-l;l;--l,tb+=sz,++pd){
pos=*pd;
s+=(ll)(pr[pos-1]-pl[pos-1]+(int)tb[pr[pos]-pr[pos-1]-1]);
}
return s;
}
ll invMR(int r,int lb,int rb){
int pos,*pl,*pr,*pd,sz=bsz;
short *tb=tbl[0]+r*sz;
ll s=(ll)((rb+1-lb)*sz*(r%sz+1));
pl=C[lb-1];
pr=C[rb];
pd=d+r;
for (r-=(rb+1)*sz-1;r;--r,tb-=sz,--pd){
pos=*pd;
s+=(ll)(-(pr[pos]-pl[pos])+(int)tb[pl[pos]-pl[pos-1]]);
}
return s;
}
ll invLR(int l,int r){
int i,j,k,s,*L,*R,*Li,*Ri,sz=bsz;
i=(l/sz)*sz;
j=(r/sz)*sz;
L=A+i;
R=A+j;
Li=b+i;
Ri=b+j;
l-=i;
r-=j;
for (i=j=k=s=0;i<sz;++i){
if (l<=Li[i]){
for (;j<sz && L[i]>R[j];++j){
if (Ri[j]<=r) ++k;
}
s+=k;
}
}
return (ll)s;
}
ll invself(int l,int r){
int i,s,nr,*ix,sz=bsz;
ix=b+(l/sz)*sz;
l%=sz;
r%=sz;
for (i=s=nr=0;i<sz;++i){
if (r<ix[i]) ++nr;
else if (l<=ix[i]) s+=nr;
}
return (ll)s;
}
public:
RangeInversions(T const *a,int n){
int i,j,tbsz,*p;
ll *q;
bsz=(int)pow((double)n,1.0/3.0)*2;
if (bsz&1) ++bsz;
nb=n/bsz;
if (bsz*nb<n) ++nb;
N=bsz*nb;
tbsz=bsz*bsz;
i=N*5+(nb+1)*(nb+1)+(tbsz/2)*nb;
j=((nb+1)*nb)/2;
B=new ll*[nb];
Bl=new int*[nb*4+1];
A=new int[i];
pB=new ll[j];
memset(A,0,i*4);
memset(pB,0,j*8);
Br=Bl+nb;
C=Br+(nb+1);
tbl=(short**)(C+nb);
b=A+N;
d=b+N;
p=d+(N+1);
C[-1]=p;p+=nb;
for (i=0;i<nb;++i,p+=nb+1) C[i]=p;
for (i=0;i<nb;++i,p+=tbsz/2) tbl[i]=(short*)p;
for (i=0;i<nb;++i,p+=bsz) Bl[i]=p;
for (i=0;i<nb;++i,p+=bsz) Br[i]=p;
for (q=pB,i=0;i<nb;q-=++i) B[i]=q,q+=nb;
for (i=0;i<N;++i) b[i]=i;
T *a_=new T[n];
memcpy(a_,a,sizeof(T)*n);
qs2(a_,b,n);
delete [] a_;
for (i=0;i<N;++i) A[b[i]]=i;
C_build(b);
p=Bl[0];
for (i=0;i<N;++i) p[i]=i;
for (i=j=0;i<nb;++i,j+=bsz) qs2(b+j,p+j,bsz);
for (i=0;i<N;++i) d[b[i]]=i;
table_build(p,b,d);
for (i=0;i<N;++i) d[i]/=bsz;
memset(p,0,N*4);
for (i=0;i<nb;++i) B_build_self(A+i*bsz,i);
for (i=0;i<N;++i) b[i]=i%bsz;
for (i=j=0;i<nb;++i,j+=bsz) qs2(A+j,b+j,bsz);
B_build(A);
}
~RangeInversions(){
delete [] B;delete [] Bl;delete [] A;delete [] pB;
}
ll inversions_query(int l,int r){
if (l>=r) return 0LL;
int lb=l/bsz;
int rb=--r/bsz;
ll s;
if (lb+1<rb){
s=invLM(l,lb+1,rb-1)+invMR(r,lb+1,rb-1)+invLR(l,r)+B[lb+1][rb-1]+(ll)(Bl[lb][l%bsz]+Br[rb][r%bsz]);
}else if (lb<rb){
s=invLR(l,r)+(ll)(Bl[lb][l%bsz]+Br[rb][r%bsz]);
}else{
s=((r+1)%bsz)? (ll)(Bl[lb][l%bsz]-Bl[lb][(r+1)%bsz])-invself(l,r):(ll)Bl[lb][l%bsz];
}
return s;
}
};
int main(void){
int i,t,l,r,*a;
ll s;
const int n=100000; // 数列長
a=new int[n];
// intの数列を作る
for (i=0;i<n;++i){
a[i]=(double)rand() / (RAND_MAX + 1) * (n);
}
// 構築
RangeInversions<int> ri(a,n);
// 転倒数クエリ10万回
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;
s=ri.inversions_query(l,r);
}
delete [] a;
return 0;
}
動的平方分割(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;
}
動的ウェーブレット行列
ウェーブレット行列を動的ビットベクトルを使って動的にしてある。
1件の大きさを最大入力値のビットサイズとしたときn件の入力に対してメモリ8n程で動く。
(1ノードに128ビット持たせてビット演算とかに負荷をかけるとメモリ5n程で動くみたいなトレードオフがある)
64ビットデータの入力に対応しているけどほとんどの処理がO(log入力件数*log最大入力値)なので最大入力値が大きくなると重くなる。
処理の大半がrankかselectの場合は前回のウェーブレット木、それ以外は今回のこれを使うというのが使い分けのひとつの目安になる。
動的ビットベクトルについて
平衡二分探索木でnode[i].rankをiノードの左部分木の「総ノード数」としてた定義を左部分木に含まれる「総ビット数」と読み替えたら大枠を変更せずに動的ビットベクトルのコードとして使うことができた。
(他のノードメンバーは、左部分木に含まれる1のビットの総数、部分ビットベクトル(64ビット)、部分ビットベクトルを何ビット使ってるかとそこに含まれる1のビットの総数)
1ビットを追加削除するときどのノードも32ビット(64ビットの半分)以上使用している状態を保つようにする(入力ビット数が32未満のときは例外)。
たとえば1ビット追加する場合
node[i].rankを頼りに追加するノード(iノードとする)を探索する。
iノードに空きがある(部分ビットベクトルの使用量が64ビット未満)ならば1ビット追加して終了。
空きがなければiノードの左部分木の最右ノードに空きがあればiノードから溢れた1ビットを移して終了。
空きがなければ最右ノードとiノードの間に新ノードを追加してそこにiノードの部分ビットベクトルの約半分を移してから平衡処理に移行して終了。
みたいな感じでやる。
#include <iostream>
#include <stdlib.h>
#include <time.h>
using namespace std;
#define uint unsigned int
#define ull unsigned long long
#define bitdel(x,k) (((((x)>>(k))>>1)<<(k))|(((1ULL<<(k))-1ULL)&(x)))
#define bitin(x,k,key) (((~((1ULL<<(k))-1ULL)&(x))<<1) | ((key)<<(k)) | ((x)&((1ULL<<(k))-1ULL)))
#define pc32a(x) (((((((x)&0x33333333)+0x11111111)>>1))&0x33333333)+(((((((x)>>2)&0x33333333)+0x11111111)>>1))&0x33333333))
#define pc32b(x) ((x)+((x)>>16))
#define pc32c(x) (((x)&0xf)+(((x)>>4)&0xf)+(((x)>>8)&0xf)+(((x)>>12)&0xf))
class WaveletMatrix2{
private:
struct node{
ull key;
uint rank,rnk1;
int l,r;
unsigned char ub,n1;
};
struct recstack{
int *p,d;
node *m;
};
int N,nptr,bitsize;
struct node *a0;
struct recstack *pstack;
inline uint pc64(ull x){
x=(x&0x5555555555555555ULL)+((x>>1)&0x5555555555555555ULL);
x=(x&0x3333333333333333ULL)+((x>>2)&0x3333333333333333ULL);
x=(x+(x>>4))&0x0f0f0f0f0f0f0f0fULL;
x=x+(x>>8);
x=x+(x>>16);
x=x+(x>>32);
return x&0x7f;
}
inline uint pc32(uint x){
x=(x&0x55555555)+((x>>1)&0x55555555);
x=(x&0x33333333)+((x>>2)&0x33333333);
x=(x+(x>>4))&0x0f0f0f0f;
x+=(x>>8);
return (x+(x>>16))&0x3f;
}
inline int ll_rot(node *a, node *i0,node *i1){
i0->rank-=i1->rank+i1->ub;
i0->rnk1-=i1->rnk1+i1->n1;
i0->l=i1->r;
i1->r=i0-a;
i1->l&=INT_MAX;
return i1-a;
}
inline int w_rot(node *a,node *i0,node *i1,node *i2,int lrrot){
if (lrrot){
i0->rank-=i1->rank+i1->ub+i2->rank+i2->ub;
i0->rnk1-=i1->rnk1+i1->n1+i2->rnk1+i2->n1;
}else{
i0->rank-=i2->rank+i2->ub;
i0->rnk1-=i2->rnk1+i2->n1;
}
i2->rank+=i1->rank+i1->ub;
i2->rnk1+=i1->rnk1+i1->n1;
i0->l=i2->r&INT_MAX;
i0->r|=i2->l&INT_MIN;
i1->l|=i2->r&INT_MIN;
i1->r=i2->l&INT_MAX;
i2->l=i1-a;
i2->r=i0-a;
return i2-a;
}
inline int rr_rot(node *a,node *i0,node *i1){
i1->rank+=i0->rank+i0->ub;
i1->rnk1+=i0->rnk1+i0->n1;
i0->r=i1->l;
i1->l=i0-a;
i1->r&=INT_MAX;
return i1-a;
}
inline int le_rot(node *a,node *i0,node *i1){
i0->rank-=i1->rank+i1->ub;
i0->rnk1-=i1->rnk1+i1->n1;
i0->l=INT_MIN|i1->r;
i1->r=INT_MIN|(i0-a);
return i1-a;
}
inline int re_rot(node *a,node *i0,node *i1){
i1->rank+=i0->rank+i0->ub;
i1->rnk1+=i0->rnk1+i0->n1;
i0->r=INT_MIN|i1->l;
i1->l=INT_MIN|(i0-a);
return i1-a;
}
void insert2(node *a,ull key,int k){
ull u,ky;
uint kysize,n1;
recstack *ps=pstack;
ps->p=&a->r;
if (key) a->rnk1++;
for (ps->m=a+a->r;ps->m!=a;++ps){
if (k<ps->m->rank){
ps[1].m=a+(*(ps[1].p=&ps->m->l)&INT_MAX);
ps->d=0;
ps->m->rank++;
ps->m->rnk1+=key;
}else if (ps->m->rank+ps->m->ub<k){
ps[1].m=a+(*(ps[1].p=&ps->m->r)&INT_MAX);
ps->d=1;
k-=ps->m->rank+ps->m->ub;
}else{
break;
}
}
if (ps->m!=a){
if (ps->m->ub<64){
k-=ps->m->rank;
ps->m->key=bitin(ps->m->key,k,key);
ps->m->ub++;
ps->m->n1+=key;
return;
}else if (ps->m->l){
node *q,*p=ps->m;
k-=p->rank;
ps->d=0;
ps[1].p=&ps->m->l;
ps[1].m=a+(ps->m->l&INT_MAX);
for (++ps;ps->m!=a;++ps){
ps->d=1;
ps[1].m=a+(*(ps[1].p=&ps->m->r)&INT_MAX);
}
q=ps[-1].m;
if (q->ub<64){
if (k){
ky=p->key&1;
p->key>>=1;
p->key=bitin(p->key,k-1,key);
p->n1+=key-ky;
}else{
ky=key;
}
p->rank++;
p->rnk1+=ky;
q->key|=ky<<q->ub;
q->ub++;
q->n1+=ky;
return;
}
if (k<32){
u=bitin(p->key,k,key);
ky=u&0xffffffff;
p->key>>=31;
p->rank+=32;
n1=pc32(ky);
p->rnk1+=n1;
p->n1-=n1-key;
}else{
ky=p->key&0xffffffff;
p->key>>=32;
p->key=bitin(p->key,k-32,key);
p->rank+=32;
n1=pc32(ky);
p->rnk1+=n1;
p->n1-=n1-key;
}
p->ub=33;
kysize=32;
}else{
ps[1].p=&ps->m->l;
if (k<32){
u=bitin(ps->m->key,k,key);
ky=u&0xffffffff;
ps->m->key>>=31;
}else{
ky=ps->m->key&0xffffffff;
ps->m->key>>=32;
ps->m->key=bitin(ps->m->key,k-32,key);
}
ps->m->ub=33;
kysize=32;
ps->m->rank=32;
n1=pc32(ky);
ps->m->rnk1=n1;
ps->m->n1-=n1-key;
ps->d=0;
++ps;
}
}else{
n1=ky=key;
kysize=1;
}
if (a->l){
int t=a->l;
a->l=(ps->m=a+(*ps->p=t))->l;
}else{
ps->m=a+(*ps->p=++(a->rank));
}
ps->m->l=ps->m->r=ps->m->rank=0;
ps->m->rnk1=0;
ps->m->key=ky;
ps->m->ub=kysize;
ps->m->n1=n1;
for (int ch=ps->m-a;pstack<=--ps;ch=ps->m-a){
if (ps->m->l&INT_MIN){
if (ps->d) ps->m->l&=INT_MAX;
else if (a[ch].l&INT_MIN) *ps->p=(*ps->p&INT_MIN)+ll_rot(a,ps->m,a+ch);
else *ps->p=(*ps->p&INT_MIN)+w_rot(a,ps->m,a+ch,a+(a[ch].r&INT_MAX),1);
break;
}else if (ps->m->r&INT_MIN){
if (!ps->d) ps->m->r&=INT_MAX;
else if (a[ch].r&INT_MIN) *ps->p=(*ps->p&INT_MIN)+rr_rot(a,ps->m,a+ch);
else *ps->p=(*ps->p&INT_MIN)+w_rot(a,a+ch,ps->m,a+(a[ch].l&INT_MAX),0);
break;
}else if (ps->d){
ps->m->r|=INT_MIN;
}else{
ps->m->l|=INT_MIN;
}
}
}
uint remove2(node *a,int k){
ull ky,u;
recstack *rempos,*ps=pstack;
node *p;
ps->p=&a->r;
for (ps->m=a+a->r;1;++ps){
if (ps->m->rank>k){
ps[1].m=a+(*(ps[1].p=&ps->m->l)&INT_MAX);
ps->m->rank--;
ps->d=0;
}else if (ps->m->rank+ps->m->ub<=k){
ps[1].m=a+(*(ps[1].p=&ps->m->r)&INT_MAX);
k-=ps->m->rank+ps->m->ub;
ps->d=1;
}else{
break;
}
}
ky=(ps->m->key>>(k-ps->m->rank))&1;
if (ky){
a->rnk1--;
for (rempos=pstack;rempos<ps;++rempos){
if (rempos->d==0) rempos->m->rnk1--;
}
}
if (ps->m->ub>32){
k-=ps->m->rank;
ps->m->key=bitdel(ps->m->key,k);
ps->m->ub--;
ps->m->n1-=ky;
return ky;
}else if (ps->m->l){
if (ps->m->r){
rempos=ps;
ps->d=0;
++ps;
for (ps->m=a+(*(ps->p=&ps[-1].m->l)&INT_MAX);ps->m->r;++ps){
ps->d=1;
ps[1].m=a+(*(ps[1].p=&ps->m->r)&INT_MAX);
}
if (ps->m->ub>32){
k-=rempos->m->rank;
rempos->m->key=bitdel(rempos->m->key,k);
rempos->m->n1-=ky;
ps->m->ub--;
u=ps->m->key>>ps->m->ub;
ps->m->key=bitdel(ps->m->key,ps->m->ub);
ps->m->n1-=u;
rempos->m->key=((rempos->m->key)<<1)|u;
rempos->m->n1+=u;
rempos->m->rank--;
rempos->m->rnk1-=u;
return ky;
}
k-=rempos->m->rank;
rempos->m->key=bitdel(rempos->m->key,k);
rempos->m->ub--;
rempos->m->n1-=ky;
ps->m->rank=rempos->m->rank-ps->m->ub;
ps->m->rnk1=rempos->m->rnk1-ps->m->n1;
ps->m->key|=rempos->m->key<<ps->m->ub;
ps->m->ub+=rempos->m->ub;
ps->m->n1+=rempos->m->n1;
*ps->p=(*ps->p&INT_MIN)+(ps->m->l&INT_MAX);
ps->m->l=rempos->m->l;
ps->m->r=rempos->m->r;
*rempos->p=(*rempos->p&INT_MIN)+(ps->m-a);
rempos->m->l=a->l;a->l=rempos->m-a;
rempos->m=ps->m;
(rempos+1)->p=&ps->m->l;
}else{
p=a+(ps->m->l&INT_MAX);
if (p->ub>32){
k-=ps->m->rank;
ps->m->key=bitdel(ps->m->key,k);
ps->m->n1-=ky;
p->ub--;
u=p->key>>p->ub;
p->n1-=u;
ps->m->key=(ps->m->key<<1)|u;
ps->m->rank--;
ps->m->rnk1-=u;
ps->m->n1+=u;
p->key=bitdel(p->key,p->ub);
return ky;
}
k-=ps->m->rank;
ps->m->key=bitdel(ps->m->key,k);
ps->m->n1-=ky;
ps->m->key=(ps->m->key<<p->ub)|p->key;
ps->m->ub+=p->ub-1;
ps->m->rank=ps->m->rnk1=0;
ps->m->n1+=p->n1;
ps->m->l=0;
p->l=a->l;a->l=p-a;
}
}else if (ps->m->r){
p=a+(ps->m->r&INT_MAX);
if (p->ub>32){
ps->m->key=bitdel(ps->m->key,k);
ps->m->ub--;
ps->m->n1-=ky;
u=p->key&1ULL;
ps->m->key|=u<<ps->m->ub;
ps->m->ub++;
ps->m->n1+=u;
p->key>>=1;
p->ub--;
p->n1-=u;
return ky;
}
ps->m->key=bitdel(ps->m->key,k);
ps->m->ub--;
ps->m->n1-=ky;
ps->m->key|=p->key<<ps->m->ub;
ps->m->ub+=p->ub;
ps->m->n1+=p->n1;
ps->m->r=0;
p->l=a->l;a->l=p-a;
}else{
if (ps==pstack){
ps->m->key=bitdel(ps->m->key,k);
ps->m->ub--;
ps->m->n1-=ky;
return ky;
}
ps->m->key=bitdel(ps->m->key,k);
ps->m->ub--;
ps->m->n1-=ky;
p=ps[-1].m;
if (p->ub>32){
if (ps[-1].d==0){
u=p->key&1;
p->key>>=1;
p->ub--;
p->rank++;
p->rnk1+=u;
p->n1-=u;
ps->m->key|=u<<ps->m->ub;
ps->m->ub++;
ps->m->n1+=u;
}else{
p->ub--;
u=p->key>>p->ub;
p->key=bitdel(p->key,p->ub);
p->n1-=u;
ps->m->key=((ps->m->key)<<1)|u;
ps->m->ub++;
ps->m->n1+=u;
}
return ky;
}else{
if (ps[-1].d==0){
p->key=(p->key<<ps->m->ub)|ps->m->key;
p->rank=p->rnk1=0;
p->ub+=ps->m->ub;
p->n1+=ps->m->n1;
}else{
p->key|=ps->m->key<<p->ub;
p->ub+=ps->m->ub;
p->n1+=ps->m->n1;
}
*ps->p&=INT_MIN;
ps->m->l=a->l; a->l=ps->m-a;
}
}
while (pstack<=--ps){
if (ps->m->l&INT_MIN){
if (ps->d){
node *m=a+(ps->m->l&INT_MAX);
if (m->l&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+ll_rot(a,ps->m,m);
}else if (m->r&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+w_rot(a,ps->m,m,a+((m->r)&INT_MAX),1);
}else{
*ps->p=(*ps->p&INT_MIN)+le_rot(a,ps->m,m);
break;
}
}else{
ps->m->l&=INT_MAX;
}
}else if (ps->m->r&INT_MIN){
if (!ps->d){
node *m=a+(ps->m->r&INT_MAX);
if (m->l&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+w_rot(a,m,ps->m,a+((m->l)&INT_MAX),0);
}else if (m->r&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+rr_rot(a,ps->m,m);
}else{
*ps->p=(*ps->p&INT_MIN)+re_rot(a,ps->m,m);
break;
}
}else{
ps->m->r&=INT_MAX;
}
}else{
if (ps->d) ps->m->l|=INT_MIN;
else ps->m->r|=INT_MIN;
break;
}
}
return ky;
}
int rank1(node *a, int k){
if (k<0) return 0;
int s=0;
node *m=a+a->r;
while (1){
if (k<m->rank){
m=a+(m->l&INT_MAX);
}else if (m->rank+m->ub<=k){
k-=m->rank+m->ub;
s+=m->rnk1+m->n1;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
k-=m->rank;
if (k<32){
uint t=((m->key<<(~k&63))>>(~k&63))&0xffffffff;
t=pc32a(t);t+=t>>16;
return s+m->rnk1+pc32c(t);
}else{
return s+m->rnk1+pc64((m->key)<<(~k&63));
}
}
int rank0(node *a,int k){
if (k<0) return 0;
int s=k+1;
node *m=a+a->r;
while (1){
if (k<m->rank){
m=a+(m->l&INT_MAX);
}else if (m->rank+m->ub<=k){
k-=m->rank+m->ub;
s-=m->rnk1+m->n1;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
k-=m->rank;
if (k<32){
uint t=((m->key<<(~k&63))>>(~k&63))&0xffffffff;
t=pc32a(t);t+=t>>16;
return s-m->rnk1-pc32c(t);
}else{
return s-m->rnk1-pc64((m->key)<<(~k&63));
}
}
void rank1_2(node *a,int k,int k2,int *ret,int *ret2){
if (k<0){
*ret=0;
if (k2<0) *ret2=0;
else *ret2=rank1(a,k2);
return;
}
int s2,s=0;
node *m2,*m=a+a->r;
while (1){
if (k2<m->rank){
m=a+(m->l&INT_MAX);
}else if (m->rank+m->ub<=k){
k-=m->rank+m->ub;
k2-=m->rank+m->ub;
s+=m->rnk1+m->n1;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
m2=m;
s2=s;
while (1){
if (k<m->rank){
m=a+(m->l&INT_MAX);
}else if (m->rank+m->ub<=k){
k-=m->rank+m->ub;
s+=m->rnk1+m->n1;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
k-=m->rank;
if (k<32){
uint u=((m->key)<<(~k&63))>>(~k&63);
u=pc32a(u);u+=u>>16;
*ret=s+m->rnk1+pc32c(u);
}else{
*ret=s+m->rnk1+pc64((m->key)<<(~k&63));
}
while (1){
if (k2<m2->rank){
m2=a+(m2->l&INT_MAX);
}else if (m2->rank+m2->ub<=k2){
k2-=m2->rank+m2->ub;
s2+=m2->rnk1+m2->n1;
m2=a+(m2->r&INT_MAX);
}else{
break;
}
}
k2-=m2->rank;
if (k2<32){
uint u=((m2->key)<<(~k2&63))>>(~k2&63);
u=pc32a(u);u+=u>>16;
*ret2=s2+m2->rnk1+pc32c(u);
}else{
*ret2=s2+m2->rnk1+pc64((m2->key)<<(~k2&63));
}
}
void rank0_2(node *a,int k,int k2,int *ret,int *ret2){
if (k<0){
*ret=0;
if (k2<0) *ret2=0;
else *ret2=rank0(a,k2);
return;
}
*ret=k+1;
*ret2=k2+1;
int s2,s=0;
node *m2,*m=a+a->r;
while (1){
if (k2<m->rank){
m=a+(m->l&INT_MAX);
}else if (m->rank+m->ub<=k){
k-=m->rank+m->ub;
k2-=m->rank+m->ub;
s+=m->rnk1+m->n1;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
m2=m;
s2=s;
while (1){
if (k<m->rank){
m=a+(m->l&INT_MAX);
}else if (m->rank+m->ub<=k){
k-=m->rank+m->ub;
s+=m->rnk1+m->n1;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
k-=m->rank;
if (k<32){
uint u=((m->key)<<(~k&63))>>(~k&63);
u=pc32a(u);u+=u>>16;
*ret-=s+m->rnk1+pc32c(u);
}else{
*ret-=s+m->rnk1+pc64((m->key)<<(~k&63));
}
while (1){
if (k2<m2->rank){
m2=a+(m2->l&INT_MAX);
}else if (m2->rank+m2->ub<=k2){
k2-=m2->rank+m2->ub;
s2+=m2->rnk1+m2->n1;
m2=a+(m2->r&INT_MAX);
}else{
break;
}
}
k2-=m2->rank;
if (k2<32){
uint u=((m2->key)<<(~k2&63))>>(~k2&63);
u=pc32a(u);u+=u>>16;
*ret2-=s2+m2->rnk1+pc32c(u);
}else{
*ret2-=s2+m2->rnk1+pc64((m2->key)<<(~k2&63));
}
}
int select1(node *a,int k){
int s=0;
node *m=a+a->r;
while (1){
if (k<m->rnk1){
m=a+(m->l&INT_MAX);
}else if (m->rnk1+m->n1<=k){
k-=m->rnk1+m->n1;
s+=m->rank+m->ub;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
k-=m->rnk1;
s+=m->rank;
uint u=m->key&0xffffffff;
u=pc32a(u);
u+=u>>16;
u=pc32c(u);
if (k<u){
u=m->key&0xffffffff;
for (++k;k;u>>=1,++s) k-=u&1;
return s-1;
}else{
k-=u;
u=m->key>>32;
for (++k;k;u>>=1,++s) k-=u&1;
return s+32-1;
}
}
int select0(node *a,int k){
int s=0;
node *m=a+a->r;
while (1){
if (k<m->rank-m->rnk1){
m=a+(m->l&INT_MAX);
}else{
k-=m->rank-m->rnk1;
s+=m->rank;
if (k<m->ub-m->n1) break;
k-=m->ub-m->n1;
s+=m->ub;
m=a+(m->r&INT_MAX);
}
}
uint u=~(m->key)&0xffffffff;
u=pc32a(u);u+=u>>16;
u=pc32c(u);
if (k<u){
u=~(m->key);
for (++k;k;u>>=1,++s) k-=u&1;
return s-1;
}else{
k-=u;
u=~(m->key)>>32;
for (++k;k;u>>=1,++s) k-=u&1;
return s+32-1;
}
}
uint bv_access(node *a,int k){
node *m=a+a->r;
while (1){
if (k<m->rank){
m=a+(m->l&INT_MAX);
}else if (m->rank+m->ub<=k){
k-=m->rank+m->ub;
m=a+(m->r&INT_MAX);
}else{
break;
}
}
return (m->key>>(k-m->rank))&1;
}
void rangemink2(node *a,int s,int e,int k,ull m, ull *ans){
int i,j,n0,n1;
rank0_2(a,s-1,e,&i,&j);
n0=j-i;
n1=e-s+1-n0;
if (k<n0){
if (m>1ULL) rangemink2(a+N,i,j-1,k,m>>1,ans);
}else{
if (m>1ULL) rangemink2(a+N,i,j-1,n0-1,m>>1,ans);
k-=n0;
for (int t=n0+k;n0<=t;--t) ans[t]|=m;
if (m>1ULL){
s+=nptr-a->rnk1-i;
rangemink2(a+N,s,s+n1-1,k,m>>1,ans+n0);
}
}
}
void rangemaxk2(node *a,int s,int e,int k,ull m, ull *ans){
int i,j,n0,n1;
rank1_2(a,s-1,e,&i,&j);
n1=j-i;
n0=e-s+1-n1;
if (k<n1){
for (int t=0;t<=k;++t) ans[t]|=m;
if (m>1ULL){
s=nptr-a->rnk1+i;
rangemaxk2(a+N,s,s+n1-1,k,m>>1,ans);
}
}else{
for (int t=0;t<n1;++t) ans[t]|=m;
if (m>1ULL){
int t=nptr-a->rnk1+i;
rangemaxk2(a+N,t,t+n1-1,n1-1,m>>1,ans);
}
k-=n1;
if (m>1ULL){
s-=i;
rangemaxk2(a+N,s,s+n0-1,k,m>>1,ans+n1);
}
}
}
public:
// n:最大入力件数 v:最大入力値
WaveletMatrix2(int n,ull v){
ull u=v;
N=n/32+4;
for (bitsize=0;u;u>>=1) ++bitsize;
a0=new struct node[N*bitsize];
pstack=new struct recstack[50];
node *p=a0;
for (int i=0;i<bitsize;++i,p+=N){
p->key=0ULL;
p->l=p->r=p->rank=p->rnk1=0;
p->ub=p->n1=0;
}
nptr=0;
}
~WaveletMatrix2(){
delete [] a0,delete [] pstack;
}
// i番目に位置する値を返す
ull access(int i){
int j,k;
node* p=a0;
ull u=0ULL;
for (j=bitsize;j--;p+=N){
if (bv_access(p,i)){
u|=1ULL<<j;
i=nptr-p->rnk1+rank1(p,i-1);
}else{
i=rank0(p,i-1);
}
}
return u;
}
// 値xをk番目の位置に追加
void insert(int k,ull key){
int i,j;
node *p=a0;
ull u;
for (i=bitsize;i--;p+=N){
if ((key>>i)&1){
j=nptr-p->rnk1+rank1(p,k-1);
u=1ULL;
}else{
j=rank0(p,k-1);
u=0ULL;
}
insert2(p,u,k);
k=j;
}
++nptr;
}
// k番目に位置する値を削除
void remove(int k){
int i,j;
node *p=a0;
for (i=bitsize;i--;p+=N){
if (remove2(p,k)){
j=nptr-(p->rnk1+1)+rank1(p,k-1);
}else{
j=rank0(p,k-1);
}
k=j;
}
--nptr;
}
// [0,i]における値xの出現回数
int rank(int i,ull x){
int l=0;
int r=i;
node *p=a0;
for (int j=bitsize;j--;p+=N){
if ((x>>j)&1){
rank1_2(p,l-1,r,&l,&r);
l+=nptr-p->rnk1;
r+=nptr-p->rnk1-1;
}else{
rank0_2(p,l-1,r,&l,&r);
--r;
}
}
return r-l+1;
}
// [s,e]におけるrank
int rangerank(int s,int e,ull x){
int l=s;
int r=e;
node *p=a0;
for (int j=bitsize;j--;p+=N){
if ((x>>j)&1){
rank1_2(p,l-1,r,&l,&r);
l+=nptr-p->rnk1;
r+=nptr-p->rnk1-1;
}else{
rank0_2(p,l-1,r,&l,&r);
--r;
}
}
return r-l+1;
}
// i番目に値xが出現する位置(存在しない場合は-1を返す)
int select(int i,ull x){
int j,l,r;
node *p=a0;
l=0;r=nptr-1;
for (j=bitsize;j--;p+=N){
if ((x>>j)&1){
rank1_2(p,l-1,r,&l,&r);
l+=nptr-p->rnk1;
r+=nptr-p->rnk1-1;
}else{
rank0_2(p,l-1,r,&l,&r);
--r;
}
}
if (l+i>r) return -1;
i+=l;
p-=N;
for (j=0;j<bitsize;++j,p-=N){
if ((x>>j)&1){
i-=nptr-p->rnk1;
i=select1(p,i);
}else{
i=select0(p,i);
}
}
return i;
}
// s番目の位置を先頭としたときのselect
int rangeselect(int s,int i,ull x){
int j,l,r;
node *p=a0;
l=s;r=nptr-1;
for (j=bitsize;j--;p+=N){
if ((x>>j)&1){
rank1_2(p,l-1,r,&l,&r);
l+=nptr-p->rnk1;
r+=nptr-p->rnk1-1;
}else{
rank0_2(p,l-1,r,&l,&r);
--r;
}
}
if (l+i>r) return -1;
i+=l;
p-=N;
for (j=0;j<bitsize;++j,p-=N){
if ((x>>j)&1){
i-=nptr-p->rnk1;
i=select1(p,i);
}else{
i=select0(p,i);
}
}
return i;
}
// [s,e]をソートしたときのk番目に位置する値
ull quantile(int s,int e,int k){
int i,t,t2,n0;
node *p=a0;
ull u=0ULL;
for (i=bitsize;i--;p+=N){
rank1_2(p,s-1,e,&t,&t2);
n0=e-s+1-(t2-t);
if (k<n0){
s-=t;
e=s+n0-1;
}else{
k-=n0;
u|=1ULL<<i;
s=nptr-p->rnk1+t;
e=s+t2-t-1;
}
}
return u;
}
// [s,e]におけるx未満の要素数
int rank_less_than(int s,int e,ull x){
int i,t,t2,j=0;
node *p=a0;
for (i=bitsize;i--;p+=N){
rank0_2(p,s-1,e,&t,&t2);
if ((x>>i)&1){
j+=t2-t;
s+=nptr-p->rnk1-t;
e+=nptr-p->rnk1-t2;
}else{
s=t;
e=t2-1;
}
}
return j;
}
// [s,e]をソートしたときの小さいほうからk+1個をans配列に返す
void rangemink(int s,int e,int k,ull *ans){
for (int i=0;i<=k;++i) ans[i]=0ULL;
rangemink2(a0,s,e,k,1ULL<<(bitsize-1),ans);
}
// [s,e]をソートしたときの大きいほうからk+1個をans配列に返す
void rangemaxk(int s,int e,int k,ull *ans){
for (int i=0;i<=k;++i) ans[i]=0ULL;
rangemaxk2(a0,s,e,k,1ULL<<(bitsize-1),ans);
}
};
int main(void){
int i,n,t;
uint u,x;
// 最大入力件数
n=100000;
// 最大入力値
const ull v=1000000000;
WaveletMatrix2 wm(n,v);
// m種類の入力値
const int m=1000;
ull *c=new ull[m];
for (i=0;i<m;++i) c[i]=(double)rand() / (RAND_MAX + 1) * v;
// 初期化
for (i=0;i<n;++i){
t=(double)rand() / (RAND_MAX + 1) * m;
x=c[t];
wm.insert(i,x);
}
wm.remove(10);
wm.insert(x,10);
u=wm.access(10);
u=wm.quantile(10,100,10);
t=wm.rank(1000,x);
t=wm.select(10,x);
t=wm.rangerank(10,1000,x);
t=wm.rangeselect(1000,0,x);
t=wm.rank_less_than(10,1000,x);
ull ans[10];
wm.rangemink(10,1000,9,ans);
//wm.rangemaxk(10,1000,9,ans);
for (i=0;i<10;++i) cout<<ans[i]<<endl;
delete [] c;
return 0;
}
動的ウェーブレット木
ウェーブレット木的データ構造を平衡二分探索木を使って動的にしてある。
まだ書いてる途中でたぶんバグも残ってるので注意。
#include <iostream>
#include <stdlib.h>
#include <time.h>
using namespace std;
#define uint unsigned int
class WaveletTree2{
private:
struct node{
uint key;
int l,r,p,rank;
};
struct node2{
int l,r,rank;
};
struct recstack{
int *p,d;
node *m;
node2 *m2;
};
int N,nptr;
uint *wrk,skey,mkey;
struct node *a0;
struct node2 *b0,*b3,*stack2;
struct recstack *pstack;
void qs(uint *argl,uint *argr){
if (argr-argl<16){
for (uint t,k,*r,*l=argl+1;l<=argr;++l){
for (k=*(r=l);argl<r && *(r-1)>k;--r) *r=*(r-1);
*r=k;
}
}else{
uint *l=argl;
uint *r=argr;
for (uint t,k=*(l+((r-l)>>1));l<=r;){
while (*l<k) ++l;
while (k<*r) --r;
if (l<=r){
t=*l,*l=*r,*r=t;
++l,--r;
}
}
if (argl<l-1) quick_sort(argl,l-1);
if (l<argr) quick_sort(l,argr);
}
}
int selectk(uint *p,int argr,int k){
int l,r;
uint u;
while (1){
if (argr<16){
for (l=1;l<=argr;++l){
u=p[l];
for (r=l;r && p[r-1]>u;--r) p[r]=p[r-1];
p[r]=u;
}
return p[k];
}else{
l=0;r=argr;
u=p[r>>1];
while (l<=r){
while (p[l]<u) ++l;
while (p[r]>u) --r;
if (l<=r){
uint t=p[l];p[l]=p[r];p[r]=t;
++l;--r;
}
}
if (k<l){
argr=l-1;
}else{
p+=l;
k-=l;
argr-=l;
}
}
}
}
inline int ll_rot(node *a, node *i0,node *i1,int ip){
i0->rank-=i1->rank+1;
if (i1->r) a[i1->r&INT_MAX].p=i0-a;
i0->l=i1->r;
i1->r=i0-a;
i1->l&=INT_MAX;
i0->p=i1-a;
i1->p=ip;
return i1-a;
}
inline int b_ll_rot(node2 *a,node2 *i0,node2 *i1){
i0->rank-=i1->rank+1;
i0->l=i1->r;
i1->r=i0-a;
i1->l&=INT_MAX;
return i1-a;
}
inline int w_rot(node *a,node *i0,node *i1,node *i2,int ip,int lrrot){
if (lrrot) i0->rank-=i1->rank+i2->rank+2;
else i0->rank-=i2->rank+1;
if (i2->l) a[i2->l&INT_MAX].p=i1-a;
if (i2->r) a[i2->r&INT_MAX].p=i0-a;
i2->rank+=i1->rank+1;
i0->l=i2->r&INT_MAX;
i0->r|=i2->l&INT_MIN;
i1->l|=i2->r&INT_MIN;
i1->r=i2->l&INT_MAX;
i2->l=i1-a;
i2->r=i0-a;
i0->p=i1->p=i2-a;
i2->p=ip;
return i2-a;
}
inline int b_w_rot(node2 *a,node2 *i0,node2 *i1,node2 *i2,int lrrot){
if (lrrot) i0->rank-=i1->rank+i2->rank+2;
else i0->rank-=i2->rank+1;
i2->rank+=i1->rank+1;
i0->l=i2->r&INT_MAX;
i0->r|=i2->l&INT_MIN;
i1->l|=i2->r&INT_MIN;
i1->r=i2->l&INT_MAX;
i2->l=i1-a;
i2->r=i0-a;
return i2-a;
}
inline int rr_rot(node *a,node *i0,node *i1,int ip){
i1->rank+=i0->rank+1;
if (i1->l) a[i1->l&INT_MAX].p=i0-a;
i0->r=i1->l;
i1->l=i0-a;
i1->r&=INT_MAX;
i0->p=i1-a;
i1->p=ip;
return i1-a;
}
inline int b_rr_rot(node2 *a,node2 *i0,node2 *i1){
i1->rank+=i0->rank+1;
i0->r=i1->l;
i1->l=i0-a;
i1->r&=INT_MAX;
return i1-a;
}
inline int le_rot(node *a,node *i0,node *i1,int ip){
i0->rank-=i1->rank+1;
a[i1->r].p=i0-a;
i0->p=i1-a;
i0->l=INT_MIN|i1->r;
i1->r=INT_MIN|(i0-a);
i1->p=ip;
return i1-a;
}
inline int b_le_rot(node2 *a,node2 *i0,node2 *i1){
i0->rank-=i1->rank+1;
i0->l=INT_MIN|i1->r;
i1->r=INT_MIN|(i0-a);
return i1-a;
}
inline int re_rot(node *a,node *i0,node *i1,int ip){
i1->rank+=i0->rank+1;
a[i1->l].p=i0-a;
i0->p=i1-a;
i0->r=INT_MIN|i1->l;
i1->l=INT_MIN|(i0-a);
i1->p=ip;
return i1-a;
}
inline int b_re_rot(node2 *a,node2 *i0,node2 *i1){
i1->rank+=i0->rank+1;
i0->r=INT_MIN|i1->l;
i1->l=INT_MIN|(i0-a);
return i1-a;
}
int b_lbsrc2(node2 *a,uint x,int m,int k){
int l=0;
for (int i=a->r;i;){
uint u=a0[i].key>>m;
if (u==x){
int t=n2i(a0,i);
if (t==k) return l+a[i].rank;
if (k<t){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}else if (x<u){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}
return l;
}
int b_rbsrc2(node2 *a,uint x,int m,int k,int i){
int l=0;
while (i){
uint u=a0[i].key>>m;
if (u==x){
int t=n2i(a0,i);
if (t==k) return l+a[i].rank;
if (k<t){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}else if (x<u){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}
return l-1;
}
int b_lrbsrc2(node2 *a,uint x,int m,int s,int e,int *lpos=NULL,int *rpos=NULL){
int i,i2,t,t2,l,l2;
i=a->r;
l=0;
while (i && (a0[i].key>>m)!=x){
if (x<(a0[i].key>>m)){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}
if (i==0) return 0;
t2=t=n2i(a0,i);
i2=i;
l2=l;
do{
if (s==t) {l+=a[i].rank;break;}
if (s<t){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
while (i && (a0[i].key>>m)!=x){
if (x<(a0[i].key>>m)){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}
if (!i) break;
t=n2i(a0,i);
}while (1);
do{
if (e==t2) {l2+=a[i2].rank+1;break;}
if (e<t2){
i2=a[i2].l&INT_MAX;
}else{
l2+=a[i2].rank+1;
i2=a[i2].r&INT_MAX;
}
while (i2 && (a0[i2].key>>m)!=x){
if (x<(a0[i2].key>>m)){
i2=a[i2].l&INT_MAX;
}else{
l2+=a[i2].rank+1;
i2=a[i2].r&INT_MAX;
}
}
if (!i2) break;
t2=n2i(a0,i2);
}while (1);
if (lpos!=NULL){
*lpos=l;
*rpos=l2-1;
}
return l2-l;
}
int b_lbsrc_p0(node2 *a,uint x,int i){
int l=0;
while (i){
if (x<=a0[i].key){
i=a[i].l&INT_MAX;
}else{
l+=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}
return l;
}
void get_sequence(node2 *a,int s,int e,int argr,uint *w,node2 *st){
int i,k,l,r;
node2 *sp=st+1;
i=st->rank=a->r;st->l=0;st->r=argr;
while (st<sp){
--sp;
i=sp->rank;
l=sp->l;r=sp->r;
k=l+a[i].rank;
if (s<=k && k<=e) *w++=i;
if (l<=e && s<=r){
if (a[i].r){
sp->rank=a[i].r&INT_MAX;
sp->l=k+1;sp->r=r;
++sp;
}
if (a[i].l){
sp->rank=a[i].l&INT_MAX;
sp->l=l;sp->r=k-1;
++sp;
}
}
}
}
inline int n2i(node *a,int i){
int k=a[i].rank;
for (int c=i;i=a[i].p;c=i){
if ((a[i].r&INT_MAX)==c) k+=a[i].rank+1;
}
return k;
}
int get_ca(node2 *a,uint x,int m){
int i=a->r;
while (i && (a0[i].key>>m)!=x){
if (x<(a0[i].key>>m)) i=a[i].l&INT_MAX;
else i=a[i].r&INT_MAX;
}
return i;
}
int get_ca(node2 *a,uint x){
int i=a->r;
while (i && a0[i].key!=x){
if (x<a0[i].key) i=a[i].l&INT_MAX;
else i=a[i].r&INT_MAX;
}
return i;
}
inline uint getseq_selectk(node2 *p,int l,int r,int k,uint *w){
get_sequence(p,l,r,nptr-1,w,stack2);
for (uint *q=w+r-l;w<=q;--q) *q=a0[*q].key;
return selectk(w,r-l,k)&mkey;
}
uint nxt(node2 *a,uint x,int m){
if ((x&255)==255) return 256;
uint u,v=256;
for (int i=a->r;i;){
u=a0[i].key>>m;
if (x<u){
v=u;
i=a[i].l&INT_MAX;
}else{
i=a[i].r&INT_MAX;
}
}
return v&255;
}
uint pre(node2 *a,uint x,int m){
if ((x&255)==0) return 256;
uint u,v=256;
for (int i=a->r;i;){
u=a0[i].key>>m;
if (x<=u){
i=a[i].l&INT_MAX;
}else{
v=u;
i=a[i].r&INT_MAX;
}
}
return v&255;
}
uint quantile_d(int s,int e,int k){
int l,r,t;
uint u=0;
node2 *p=b0;
for (int i=24;0<=i;i-=8,p+=N){
for (uint v=255;v<256;k-=t){
t=b_lrbsrc2(p,(u>>i)|v,i,s,e,&l,&r);
if (t>k){
u|=v<<i;
break;
}
v=pre(p,(u>>i)|v,i);
}
if (i && t<1024) return getseq_selectk(p,l,r,r-l-k,wrk);
}
return u&mkey;
}
void insert2(node *a,uint key,int k){
recstack *ps=pstack;
ps->p=&a->r;
for (ps->m=a+a->r;ps->m!=a;++ps){
if (k<=ps->m->rank){
ps[1].m=a+(*(ps[1].p=&ps->m->l)&INT_MAX);
ps->m->rank++;
ps->d=0;
}else{
ps[1].m=a+(*(ps[1].p=&ps->m->r)&INT_MAX);
k-=ps->m->rank+1;
ps->d=1;
}
}
if (a->p){
int t=a->p;
a->p=(ps->m=a+(*ps->p=t))->l;
}else{
ps->m=a+(*ps->p=++(a->rank));
}
ps->m->key=key;
ps->m->l=ps->m->r=ps->m->rank=0;
ps->m->p= pstack<ps ? ps[-1].m-a: 0;
*ps->p=ps->m-a;
for (int ch=ps->m-a;pstack<=--ps;ch=ps->m-a){
if (ps->m->l&INT_MIN){
if (ps->d) ps->m->l&=INT_MAX;
else if (a[ch].l&INT_MIN) *ps->p=(*ps->p&INT_MIN)+ll_rot(a,ps->m,a+ch,ps->m->p);
else *ps->p=(*ps->p&INT_MIN)+w_rot(a,ps->m,a+ch,a+(a[ch].r&INT_MAX),ps->m->p,1);
break;
}else if (ps->m->r&INT_MIN){
if (!ps->d) ps->m->r&=INT_MAX;
else if (a[ch].r&INT_MIN) *ps->p=(*ps->p&INT_MIN)+rr_rot(a,ps->m,a+ch,ps->m->p);
else *ps->p=(*ps->p&INT_MIN)+w_rot(a,a+ch,ps->m,a+(a[ch].l&INT_MAX),ps->m->p,0);
break;
}else if (ps->d){
ps->m->r|=INT_MIN;
}else{
ps->m->l|=INT_MIN;
}
}
}
void b_insert2(node2 *a,int k){
recstack *ps=pstack;
ps->p=&a->r;
for (ps->m2=a+a->r;ps->m2!=a;++ps){
if (k<=ps->m2->rank){
ps[1].m2=a+(*(ps[1].p=&ps->m2->l)&INT_MAX);
ps->m2->rank++;
ps->d=0;
}else{
ps[1].m2=a+(*(ps[1].p=&ps->m2->r)&INT_MAX);
k-=ps->m2->rank+1;
ps->d=1;
}
}
if (a->l){
int t=a->l;
a->l=(ps->m2=a+(*ps->p=t))->l;
}else{
ps->m2=a+(*ps->p=++(a->rank));
}
ps->m2->l=ps->m2->r=ps->m2->rank=0;
*ps->p=ps->m2-a;
for (int ch=ps->m2-a;pstack<=--ps;ch=ps->m2-a){
if (ps->m2->l&INT_MIN){
if (ps->d) ps->m2->l&=INT_MAX;
else if (a[ch].l&INT_MIN) *ps->p=(*ps->p&INT_MIN)+b_ll_rot(a,ps->m2,a+ch);
else *ps->p=(*ps->p&INT_MIN)+b_w_rot(a,ps->m2,a+ch,a+(a[ch].r&INT_MAX),1);
break;
}else if (ps->m2->r&INT_MIN){
if (!ps->d) ps->m2->r&=INT_MAX;
else if (a[ch].r&INT_MIN) *ps->p=(*ps->p&INT_MIN)+b_rr_rot(a,ps->m2,a+ch);
else *ps->p=(*ps->p&INT_MIN)+b_w_rot(a,a+ch,ps->m2,a+(a[ch].l&INT_MAX),0);
break;
}else if (ps->d){
ps->m2->r|=INT_MIN;
}else{
ps->m2->l|=INT_MIN;
}
}
}
void remove2(node *a,int k){
recstack *rempos,*ps=pstack;
ps->p=&a->r;
for (ps->m=a+a->r;ps->m->rank!=k;ps->m=a+(*ps->p&INT_MAX)){
if (k<ps->m->rank){
ps[1].p=&ps->m->l;
ps->m->rank--;
ps->d=0;
}else{
ps[1].p=&ps->m->r;
k-=ps->m->rank+1;
ps->d=1;
}
++ps;
}
if (ps->m->l){
if (ps->m->r){
rempos=ps;
ps->d=0;
ps->m->rank--;
++ps;
for (ps->m=a+(*(ps->p=&ps[-1].m->l)&INT_MAX);ps->m->r;++ps){
ps->d=1;
ps[1].m=a+(*(ps[1].p=&ps->m->r)&INT_MAX);
}
*ps->p=(*ps->p&INT_MIN)+(ps->m->l&INT_MAX);
if (ps->m->l) a[ps->m->l&INT_MAX].p=ps->m->p;
ps->m->l=rempos->m->l;
ps->m->r=rempos->m->r;
ps->m->p=rempos->m->p;
ps->m->rank=rempos->m->rank;
*rempos->p=(*rempos->p&INT_MIN)+(ps->m-a);
if (ps->m->l){
a[ps->m->l&INT_MAX].p=ps->m-a;
}
a[ps->m->r&INT_MAX].p=ps->m-a;
rempos->m->l=a->p; a->p=rempos->m-a;
rempos->m=ps->m;
(rempos+1)->p=&ps->m->l;
}else{
*ps->p=(*ps->p&INT_MIN)+(ps->m->l&INT_MAX);
a[ps->m->l&INT_MAX].p=ps->m->p;
ps->m->l=a->p; a->p=ps->m-a;
}
}else if (ps->m->r){
*ps->p=(*ps->p&INT_MIN)+(ps->m->r&INT_MAX);
a[ps->m->r&INT_MAX].p=ps->m->p;
ps->m->l=a->p; a->p=ps->m-a;
}else{
*ps->p&=INT_MIN;
ps->m->l=a->p; a->p=ps->m-a;
}
while (pstack<=--ps){
if (ps->m->l&INT_MIN){
if (ps->d){
node *m=a+(ps->m->l&INT_MAX);
if (m->l&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+ll_rot(a,ps->m,m,ps->m->p);
}else if (m->r&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+w_rot(a,ps->m,m,a+((m->r)&INT_MAX),ps->m->p,1);
}else{
*ps->p=(*ps->p&INT_MIN)+le_rot(a,ps->m,m,ps->m->p);
break;
}
}else{
ps->m->l&=INT_MAX;
}
}else if (ps->m->r&INT_MIN){
if (!ps->d){
node *m=a+(ps->m->r&INT_MAX);
if (m->l&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+w_rot(a,m,ps->m,a+((m->l)&INT_MAX),ps->m->p,0);
}else if (m->r&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+rr_rot(a,ps->m,m,ps->m->p);
}else{
*ps->p=(*ps->p&INT_MIN)+re_rot(a,ps->m,m,ps->m->p);
break;
}
}else{
ps->m->r&=INT_MAX;
}
}else{
//*((&ps->m->r)-(ps->d))|=INT_MIN;return;
if (ps->d) ps->m->l|=INT_MIN;
else ps->m->r|=INT_MIN;
break;
}
}
}
void b_remove2(node2 *a,int k){
recstack *rempos,*ps=pstack;
ps->p=&a->r;
for (ps->m2=a+a->r;ps->m2->rank!=k;ps->m2=a+(*ps->p&INT_MAX)){
if (k<ps->m2->rank){
ps[1].p=&ps->m2->l;
ps->m2->rank--;
ps->d=0;
}else{
ps[1].p=&ps->m2->r;
k-=ps->m2->rank+1;
ps->d=1;
}
++ps;
}
if (ps->m2->l){
if (ps->m2->r){
rempos=ps;
ps->d=0;
ps->m2->rank--;
++ps;
for (ps->m2=a+(*(ps->p=&ps[-1].m2->l)&INT_MAX);ps->m2->r;++ps){
ps->d=1;
ps[1].m2=a+(*(ps[1].p=&ps->m2->r)&INT_MAX);
}
*ps->p=(*ps->p&INT_MIN)+(ps->m2->l&INT_MAX);
ps->m2->l=rempos->m2->l;
ps->m2->r=rempos->m2->r;
ps->m2->rank=rempos->m2->rank;
*rempos->p=(*rempos->p&INT_MIN)+(ps->m2-a);
rempos->m2->l=a->l;a->l=rempos->m2-a;
rempos->m2=ps->m2;
(rempos+1)->p=&ps->m2->l;
}else{
*ps->p=(*ps->p&INT_MIN)+(ps->m2->l&INT_MAX);
ps->m2->l=a->l;a->l=ps->m2-a;
}
}else if (ps->m2->r){
*ps->p=(*ps->p&INT_MIN)+(ps->m2->r&INT_MAX);
ps->m2->l=a->l; a->l=ps->m2-a;
}else{
*ps->p&=INT_MIN;
ps->m2->l=a->l; a->l=ps->m2-a;
}
while (pstack<=--ps){
if (ps->m2->l&INT_MIN){
if (ps->d){
node2 *m=a+(ps->m2->l&INT_MAX);
if (m->l&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+b_ll_rot(a,ps->m2,m);
}else if (m->r&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+b_w_rot(a,ps->m2,m,a+((m->r)&INT_MAX),1);
}else{
*ps->p=(*ps->p&INT_MIN)+b_le_rot(a,ps->m2,m);
break;
}
}else{
ps->m2->l&=INT_MAX;
}
}else if (ps->m2->r&INT_MIN){
if (!ps->d){
node2 *m=a+(ps->m2->r&INT_MAX);
if (m->l&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+b_w_rot(a,m,ps->m2,a+((m->l)&INT_MAX),0);
}else if (m->r&INT_MIN){
*ps->p=(*ps->p&INT_MIN)+b_rr_rot(a,ps->m2,m);
}else{
*ps->p=(*ps->p&INT_MIN)+b_re_rot(a,ps->m2,m);
break;
}
}else{
ps->m2->r&=INT_MAX;
}
}else{
if (ps->d) ps->m2->l|=INT_MIN;
else ps->m2->r|=INT_MIN;
break;
}
}
}
inline uint search(node *a,int k){
node *m=a+a->r;
while (m->rank!=k){
if (k<m->rank){
m=a+(m->l&INT_MAX);
}else{
k-=m->rank+1;
m=a+(m->r&INT_MAX);
}
}
return m->key;
}
inline int b_search(node2 *a,int k){
int i=a->r;
while (a[i].rank!=k){
if (k<a[i].rank){
i=a[i].l&INT_MAX;
}else{
k-=a[i].rank+1;
i=a[i].r&INT_MAX;
}
}
return i;
}
int init2_a_h(int i,int l,int r,int ip,int *h){
int t0,t1,h0,h1;
if (l<i){
t0=(l+i-1)>>1;
a0[i].l=t0;
t0=init2_a_h(t0,l,i-1,i,&h0);
}else{
t0=a0[i].l=0;
h0=-1;
}
if (i<r){
t1=(i+1+r)>>1;
a0[i].r=t1;
t1=init2_a_h(t1,i+1,r,i,&h1);
}else{
t1=a0[i].r=0;
h1=-1;
}
a0[i].p=ip;
a0[i].rank=t0;
if (h0<h1){
a0[i].r|=INT_MIN;
*h=h1+1;
}else if (h0>h1){
a0[i].l|=INT_MIN;
*h=h0+1;
}else{
*h=h0+1;
}
return t0+t1+1;
}
void init2_b1b2(node2 *b,int i,int l,int r,int j,int *h,int *rank){
int h0,h1,rnk0,rnk1;
if (l<i){
int t=(l+i-1)>>1;
b[j].l=a0[t].rank;
init2_b1b2(b,t,l,i-1,b[j].l,&h0,&rnk0);
}else{
b[j].l=0;
h0=-1;rnk0=0;
}
if (i<r){
int t=(i+1+r)>>1;
b[j].r=a0[t].rank;
init2_b1b2(b,t,i+1,r,b[j].r,&h1,&rnk1);
}else{
b[j].r=0;
h1=-1;rnk1=0;
}
b[j].rank=rnk0;
*rank=rnk0+rnk1+1;
if (h0<h1){
b[j].r|=INT_MIN;
*h=h1+1;
}else if (h0>h1){
b[j].l|=INT_MIN;
*h=h0+1;
}else{
*h=h0+1;
}
}
void init2_b(node2 *a,int n,int s,uint *c){
int i,j,t;
uint u;
if (n<32){
for (i=1;i<n;++i){
t=a[i].l;
u=(a0[t].key>>s)&255;
for (j=i;0<j;--j){
if (((a0[a[j-1].l].key>>s)&255)<=u) break;
a[j].l=a[j-1].l;
}
a[j].l=t;
}
for (i=0;i<n;++i) a[i].rank=a[i].l;
if (s){
node2 *p=a+N;
for (i=0;i<n;++i) p[i].rank=p[i].l=a[i].l;
for (i=0;i<n;i=j){
u=(a0[a[i].l].key>>s)&255;
for (j=i+1;j<n;++j){
if (((a0[a[j].l].key>>s)&255)!=u) break;
}
init2_b(p+i,j-i,s-8,c+258);
}
}
}else{
for (i=0;i<=256;++i) c[i]=0;
for (i=0;i<n;++i) ++c[a[i].r=((a0[a[i].l].key)>>s)&255];
for (i=1;i<=256;++i) c[i]+=c[i-1];
for (i=n-1;0<=i;--i) a[--c[a[i].r]].rank=a[i].l;
if (s){
node2 *p=a+N;
for (i=0;i<n;++i) p[i].l=a[i].rank;
for (i=0;i<256;++i){
if (c[i]<c[i+1]){
init2_b(p+c[i],c[i+1]-c[i],s-8,c+258);
}
}
}
}
}
public:
// n:最大入力件数 v:最大入力値
WaveletTree2(int n, uint v=4294967295){
N=n+1;
a0=new struct node[N];
b0=new struct node2[N*4+50];
pstack=new struct recstack[50];
wrk=new uint[258+1026];
b3=b0+N*3;
stack2=b3+N;
a0->key=a0->l=a0->p=a0->r=a0->rank=0;
for (node2 *p=b0;p<stack2;p+=N) p->l=p->r=p->rank=0;
for (skey=0;(v<<skey)<0x80000000;++skey);
mkey=UINT_MAX>>skey;
nptr=0;
}
~WaveletTree2(){
delete [] a0,delete [] pstack,delete [] wrk,delete [] b0;
}
// [0,i]における値xの出現回数
int rank(int i,uint x){
x+=(x<<skey)&~mkey;
int j=get_ca(b3,x);
return b_rbsrc2(b3,x,0,i,j)+1-b_lbsrc_p0(b3,x,j);
}
// [s,e]におけるrank
int rangerank(int s,int e,uint x){
return b_lrbsrc2(b3,x+((x<<skey)&~mkey),0,s,e);
}
// i番目に値xが出現する位置(存在しない場合は-1を返す)
int select(int i,uint x){
x+=(x<<skey)&~mkey;
int l=b_lbsrc_p0(b3,x,b3->r);
if ((l+i>=nptr) || (a0[l=b_search(b3,l+i)].key!=x)) return -1;
return n2i(a0,l);
}
// s番目の位置を先頭としたときのselect
int rangeselect(int s,int i,uint x){
x+=(x<<skey)&~mkey;
int l=b_lbsrc2(b3,x,0,s);
if ((l+i>=nptr) || (a0[l=b_search(b3,l+i)].key!=x)) return -1;
return n2i(a0,l);
}
// [s,e]をソートしたときのk番目に位置する値
uint quantile(int s,int e,int k){
if (k>((e-s)>>1)) return quantile_d(s,e,e-s-k);
int l,r,t;
uint u=0;
node2 *p=b0;
for (int i=24;0<=i;i-=8,p+=N){
for (uint v=0;v<256;k-=t){
t=b_lrbsrc2(p,(u>>i)|v,i,s,e,&l,&r);
if (t>k){
u|=v<<i;
break;
}
v=nxt(p,(u>>i)|v,i);
}
if (i && t<1024) return getseq_selectk(p,l,r,k,wrk);
}
return u&mkey;
}
// [s,e]をソートしたときの小さいほうからk+1個をans配列に返す
void rangemink(int s,int e,int k,uint *ans){
int l,r,t,k2=k;
uint u=0;
node2 *p=b0;
uint *q=ans;
for (int i=24;0<=i;i-=8,p+=N){
for (uint v=0;v<256;++v,k-=t){
t=b_lrbsrc2(p,(u>>i)|v,i,s,e,&l,&r);
if (t>k){
if (i){
u|=v<<i;
}else{
for (u=(u|v)&mkey;0<=k;--k) *q++=u;
}
break;
}else if (t){
get_sequence(p,l,r,nptr-1,q,stack2);
for (int j=t;j--;++q) *q=a0[*q].key&mkey;
}
}
if (i && t<1024){
get_sequence(p,l,r,nptr-1,wrk,stack2);
for (int j=0;j<t;++j) wrk[j]=a0[wrk[j]].key&mkey;
qs(wrk,wrk+t-1);
for (int j=0;j<=k;++j) *q++=wrk[j];
break;
}
}
qs(ans,ans+k2);
}
// [s,e]をソートしたときの大きいほうからk+1個をans配列に返す
void rangemaxk(int s,int e,int k,uint *ans){
int l,r,t,k2=k;
uint u=0;
node2 *p=b0;
uint *q=ans;
for (int i=24;0<=i;i-=8,p+=N){
for (uint v=256;v--;k-=t){
t=b_lrbsrc2(p,(u>>i)|v,i,s,e,&l,&r);
if (t>k){
if (i){
u|=v<<i;
}else{
for (u=(u|v)&mkey;0<=k;--k) *q++=u;
}
break;
}else if (t){
get_sequence(p,l,r,nptr-1,q,stack2);
for (int j=t;j--;++q){
*q=a0[*q].key&mkey;
}
}
}
if (i && t<1024){
get_sequence(p,l,r,nptr-1,wrk,stack2);
for (int j=0;j<t;++j) wrk[j]=a0[wrk[j]].key&mkey;
qs(wrk,wrk+t-1);
for (;0<=k;--k) *q++=wrk[--t];
break;
}
}
qs(ans,ans+k2);
for (l=0,r=k2;l<r;++l,--r){
u=ans[l];ans[l]=ans[r];ans[r]=u;
}
}
// 値xをk番目の位置に追加
void insert(uint x,int k){
node2 *p=b0;
x+=(x<<skey)&~mkey;
for (int i=24;0<=i;i-=8,p+=N){
b_insert2(p,b_lbsrc2(p,x>>i,i,k));
}
insert2(a0,x,k);
++nptr;
}
// k番目に位置する値を削除
uint remove(int k){
node2 *p=b0;
uint v=search(a0,k);
for (int i=24;0<=i;i-=8,p+=N){
b_remove2(p,b_lbsrc2(p,v>>i,i,k));
}
remove2(a0,k);
--nptr;
return v&mkey;
}
// k番目に位置する値を返す
uint access(int k){
int i=a0->r;
while (a0[i].rank!=k){
if (k<a0[i].rank){
i=a0[i].l&INT_MAX;
}else{
k-=a0[i].rank+1;
i=a0[i].r&INT_MAX;
}
}
return a0[i].key&mkey;
}
// 初期化用のinsert
void init(uint x,int i){
a0[i+1].key=x+((x<<skey)&~mkey);
++nptr;
}
// init後にこれを実行すると初期化完了
void init2(){
node2 *p=b0;
int i,j,x,m=a0->r=(1+nptr)>>1;
for (i=1;i<=nptr;++i) b0[i].l=i;
init2_b(b0+1,nptr,24,wrk);
for (j=4;j--;p+=N){
p->r=p[m].rank;
for (i=1;i<=nptr;++i) a0[i].rank=p[i].rank;
init2_b1b2(p,m,1,nptr,p->r,&x,&x);
p->rank=nptr;
}
init2_a_h(m,1,nptr,0,&x);
a0->rank=nptr;
}
};
int main(void){
int i,n,t;
uint u,x,*c;
// 最大入力件数
n=100000;
// 最大入力値
const uint v=1000000000;
WaveletTree2 wt(n,v);
// m種類の入力値
const int m=1000;
c=new uint[m];
for (i=0;i<m;++i) c[i]=(double)rand() / (RAND_MAX + 1) * v;
// 初期化
for (i=0;i<n;++i){
t=(double)rand() / (RAND_MAX + 1) * m;
x=c[t];
wt.init(x,i);
}
wt.init2();
wt.remove(10);
wt.insert(x,10);
u=wt.access(10);
u=wt.quantile(10,100,10);
t=wt.rank(1000,x);
t=wt.select(10,x);
t=wt.rangerank(10,1000,x);
t=wt.rangeselect(1000,0,x);
uint ans[10];
wt.rangemink(10,1000,9,ans);
//wt.rangemaxk(10,1000,9,ans);
for (i=0;i<10;++i) cout<<ans[i]<<endl;
delete [] c;
return 0;
}
テスト
テスト
int i;
// テスト int i;
