動的WM(補助データを使う)

rangereport(l,r,x,y) : 数列の区間[l,r)に含まれるx以上y未満の要素を列挙する
rangemink(l,r,k) : 数列の区間[l,r)に含まれる要素のうち 小さい方 からk個を列挙する
rangemaxk(l,r,k) : 数列の区間[l,r)に含まれる要素のうち 大きい方 からk個を列挙する

静的WM(ウェーブレット行列)によるrangereportの処理の手順はおおよそ次のようになる

report_me(i,l,r,x) : WMのi行の区間[l,r)に含まれる x以上 の要素を列挙する
report_lt(i,l,r,x) : WMのi行の区間[l,r)に含まれる x未満 の要素を列挙する
report_all(i,l,r)  : WMのi行の区間[l,r)に含まれる 全 て の要素を列挙する

x<yとする
WMの最初の行から開始し、xとyのビットが等しい間はWM上をその方向に進む
xが0、yが1になる最初のビットの行の次の行の0側に対してreport_meを呼びx以上を列挙し、1側に対してreport_ltを呼びy未満を列挙する
(最初の行が最後の行のとき((x+1=y かつ xは偶数)のとき)はxを最後の行のビットベクトルの[l,r)に含まれる0の数だけ出力して終了)
report_me、report_lt内で必要があればreport_allを呼ぶ

WMは上位ビットからの基数ソートによって構築されるが、基数ソートを行うときに使用する作業用配列を書き換えずに全て保存したものがこの場合の補助データになる

msb : 要素の最大値の最上位の1のビットの位置(1桁目を0-thビットとする)
w(msb) : 対象となる数列そのもの
w(i) : w(i+1)の要素を(i+1)-thビットをキーにしてソートした数列
のとき0≦i≦msbのw(i)を保存するということ

静的な場合、mを要素の最大値、kをrangereportの解の個数とするとrangereportをWMだけで行うとO(log(m)*k)であるが、補助データを使うとi行にいるときのreport_all(i,l,r)がO(r-l)になるためO(log(m)+k)になる

動的な場合、nを入力件数とするとrangereportをWMだけで行うとO(log(m)*log(n)*k)であるが、静的な場合の補助データを平衡木によって実現するとreport_allがO(log(n)+r-l)になるためO(log(m)*log(n)+k)になる

rangeminkとrangemaxkは探索過程で解として確定した要素はreport_allで列挙し確定していない要素は次の行以降で探索するというようにするとrangereportと同じ計算量になる(こちらのkは引数で指定される定数であるという違いはある)

rangepreport_posのように処理名に_posがついているものは値を列挙するだけでなく値の数列上の位置を返す処理でありO(log(m)*log(n)*k)である
(この計算量は補助データを使わないWMと等しいがメモリ使用量の定数倍の増大によってO(log(m)*log(n)+log(n)*k)に改善できる)

その他

上の説明は1ビット単位での基数ソートによるWMなので32ビットデータを扱う場合補助データとして静的な場合は32本の配列、動的な場合は32本の平衡木を使用する必要があるがこれがメモリ的に問題になる場合は4ビットや8ビット単位での基数ソートによるWMに変更すると使用する配列や平衡木の本数を減らすことができる

補助データをWM(補助データを使う)で構築するという高次元化が理屈上は可能であるが実用上はメモリの制約のため低次元でしか動かせない(Range Treeなどの高次元化に伴う問題と同様である)

WMは何がうれしいのか
静的な構造を自然な拡張によって動的な構造に変換してもならし計算量にならないところ

重み付きインターバルデータに対するTop-k範囲探索
これは少し前に何かで検索してたら見つけた問題でrangereportが使えないかなどと考えていたらクエリO(log(N)+k)になる気がしたのでいつかコードとかを書くかもしれない

#include 
#include 
#include 
#include 
#include 

#define uchar unsigned char
#define uint unsigned int



void qs2(uint *a,uint *b,uint n_,uint *w){
	if (!n_) return;
	uint l,r,l_,r_,m,t,sp=0;
	l_=0;r_=n_-1;

	while (1){
		l=l_;r=r_;
		m=a[(l+r)>>1];
		while (l<=r){
			while (a[l]<m) ++l;
			while (m<a[r]) --r;
			if (l<=r){
				t=a[l],a[l]=a[r],a[r]=t;
				t=b[l],b[l]=b[r],b[r]=t;
				++l;
				if (r) --r;		}
		}

		if (l_+1<l){
			if (l<r_) w[sp]=l,w[sp+1]=r_,sp+=2;
			r_=l-1;
		}else if (l<r_){
			l_=l;
		}else if (sp){
			r_=w[--sp];l_=w[--sp];
		}else{
			break;
		}
	}
}

void qs_k(uint *a,uint k,uint n_){
	if (k>=n_) return;
	uint l,r,l_,r_,m,t;
	l_=0;r_=n_-1;

	while (1){
		l=l_;
		r=r_;
		m=a[(l+r)>>1];

		while (l<=r){
			while (a[l]<m) ++l;
			while (m<a[r]) --r;
			if (l<=r){
				t=a[l],a[l]=a[r],a[r]=t;
				++l;
				if (r) --r;
			}
		}
		if (k>1;
		m=a[t];
		i=ix[t];
		while (l<=r){
			while (a[l]>SH)&255];
		}
	}

	uint array_sum(uint *p,uint n_){
		uint s=0;
		for (uint i=0;i<n_;++i) s+=p[i];
		return s;
	}

	uint node_split(uint ib,uint newi,uint internalnode){
		if (internalnode){
			uint 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{
			uint t=leafsize/2;
			memcpy(bl[newi].a,bl[ib].a+t,t*szx);
			bl[ib].n=bl[newi].n=t;
			bl[newi].pleaf=bl[ib].pleaf;bl[ib].pleaf=newi;
			return t;
		}
	}

	uint node_merge(uint lb,uint ib,uint internalnode){
		if (internalnode){
			uint t=b[lb].n;
			uint 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{
			uint t=bl[lb].n;
			uint u=bl[ib].n;
			memcpy(bl[lb].a+t,bl[ib].a,u*szx);
			bl[lb].n+=u;
			bl[lb].pleaf=bl[ib].pleaf;
			return t+u;
		}
	}

	void c_popr_shiftr(uint ib,uint l,uint r){
		if (l<r){
			uint *p=b[ib].c+l*256;
			memmove(p+256,p,(r-l)*256*4);
		}
	}

	void c_popl_shiftl(uint ib,uint l,uint r){
		if (l<r){
			uint *p=b[ib].c+l*256;
			memmove(p,p+256,(r-l)*256*4);
		}
	}

	void c_count_leaf(uint ib,uint n_,uint *ans,uint init=0){
		uint *p=bl[ib].a;
		const uint s=SH;
		if (init) memset(ans,0,256*4);
		for (uint i=0;i<n_;++i) ++ans[(p[i]>>s)&255];
	}

	uint leaf_rank(uint *p,uint n_,uint x,uint ms){
		uint s=0;
		for (;n_;--n_,++p){
			if ((*p&ms)==x) ++s;
		}
		return s;
	}

	void c_add_c(uint *c,uint *c2){
		for (uint i=0;i<256;++i) c[i]+=c2[i];
	}

	void c_sub_c(uint *c,uint *c2){
		for (uint i=0;i<256;++i) c[i]-=c2[i];
	}

	void cs_add_c(uint ib,uint i,uint n_){
		uint *p=b[ib].c+i*256;
		for (uint *q=p+256;0<n_;--n_,q+=256){
			for (uint k=0;k<256;++k) q[k]+=p[k];
		}
	}

	void cs_sub_c(uint ib,uint i,uint n_){
		uint *p=b[ib].c+i*256;
		for (uint *q=p+256;0<n_;--n_,q+=256){
			for (uint k=0;k<256;++k) q[k]-=p[k];
		}
	}

	void split_c(uint ib,uint 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(uint ib,uint rb){
		uint nl=b[ib].n;
		uint 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(uint ib,uint pos,uint internalnode){
		uint 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(uint ib,uint pos){
		uint *p=b[ib].c+pos*256;
		memmove(p,p+256,(b[ib].n-(pos+1))*256*4);
	}

	uint insert2_rnk(uint ib,uint k,uint x,uint ip,uint retrank){
		uint t,pos,ic,*p,*q,h,ret,s=0;
		const uint y=(x>>SH)&255;

		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];
			p=q=b[ib].c+(pos*256+y);
			if (pos) s+=q[-256];
			for (t=b[ib].n-pos;t;q+=256,--t) ++*q;
			ip=ib;
			ib=ic;
		}

		q=bl[ib].a;
		if (retrank){
			if (k+k<=bl[ib].n || ip==UINT_MAX){
				s+=leaf_rank(q,k,x&MS,MS);
			}else{
				s+=pos? (*p-p[-256])-1:*p-1;
				s-=leaf_rank(q+k,bl[ib].n-k,x&MS,MS);
			}
			Trank=s;
		}

		memmove(q+(k+1),q+k,(bl[ib].n-k)*szx);
		bl[ib].a[k]=x;
		ret=(++bl[ib].n==leafsize);

		for (ic=ib,ib=ip;ret && ib!=UINT_MAX;++h){
			uint 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;
	}

	uint remove2_rnk(uint ib,uint k,uint ip,uint retrank){
		uint t,pos,ic,*p,*q,rc,h,ret,x,s=0;

		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;
			p=b[ib].c+pos*256;
			ip=ib;
			ib=b[ib].a[pos];
		}

		q=bl[ib].a;
		t=Tx=q[k];
		x=(t>>SH)&255;

		if (retrank){
			if (k+k<=bl[ib].n || ip==UINT_MAX){
				s=leaf_rank(q,k,t&MS,MS);
			}else{
				p+=x;
				s=pos? (*p-p[-256])-1:*p-1;
				s-=leaf_rank(q+(k+1),bl[ib].n-(k+1),t&MS,MS);
			}
		}

		t=--bl[ib].n-k;
		memmove(q+k,q+(k+1),t*szx);
		ret=(bl[ib].n<crem_l);

		for (ic=ib,ib=ip;ib!=UINT_MAX;++h){
			pos=b[ib].pos;
			q=b[ib].c+(pos*256+x);
			if (pos) s+=q[-256];
			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;
		}
		Trank=s;
		return ret;
	}

	void freq(uint **P,uint **Q,uint pq,uint *ans){
		uint 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 freq1(uint **P,uint p,uint *ans,const uint i0,const uint i1){
		uint i,*p0,*p1,*p2,*p3,*p4,*p5,*p6;

		switch (p){
			case 1:
				p0=P[0];
				for (i=i0;i<i1;++i) ans[i]=p0[i];break;
			case 2:
				p0=P[0],p1=P[1];
				for (i=i0;i<i1;++i) ans[i]=p0[i]+p1[i];break;
			case 3:
				p0=P[0],p1=P[1],p2=P[2];
				for (i=i0;i<i1;++i) ans[i]=p0[i]+p1[i]+p2[i];break;
			case 4:
				p0=P[0],p1=P[1],p2=P[2],p3=P[3];
				for (i=i0;i<i1;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i];break;
			case 5:
				p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4];
				for (i=i0;i<i1;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i];break;
			case 6:
				p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5];
				for (i=i0;i<i1;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i];break;
			case 7:
				p0=P[0],p1=P[1],p2=P[2],p3=P[3],p4=P[4],p5=P[5],p6=P[6];
				for (i=i0;i<i1;++i) ans[i]=p0[i]+p1[i]+p2[i]+p3[i]+p4[i]+p5[i]+p6[i];break;
		}
	}

	void ini2(uint il,uint n_,uint h_){
		uint 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>SH)&255))>>4];
		if (bl[root].n==(leafsize/2)) ++root;
		bl[root].a[bl[root].n]=x;
		++bl[root].n;
	}

	void init2(uint i,uint x){
		uint j=i>>8;
		++memi[((x>>SH)&255)>>4];
		bl[j].a[i&255]=x;
		++bl[j].n;
		if (root<j) root=j;
	}

	uint init_access_x(uint i){
		return bl[i>>8].a[i&255];
	}

	void init_end(){
		if (root && bl[root].n<crem_l){
			memcpy(bl[root-1].a+bl[root-1].n,bl[root].a,bl[root].n*szx);
			bl[root-1].n+=bl[root].n;
			--root;
		}

		mi_l=root+1;
		if (root){
			for (uint i=0;i<root;++i) bl[i].pleaf=i+1;
			ini2(0,root+1,1);
			ini3(root,NULL,hi);
		}
	}

	uint access(uint k){
		uint *p,pos,ib=root;

		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}
		return bl[ib].a[k];
	}

	uint insert_rank(uint k,uint x,uint retrank){
		++memi[(x>>(SH+4))&15];

		if (insert2_rnk(root,k,x,UINT_MAX,retrank)){
			uint ro=getindex(1);
			uint 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;
		}
		return Trank;
	}

	uint remove_rank(uint k,uint retrank){
		if (remove2_rnk(root,k,UINT_MAX,retrank) && 1>(SH+4))&15];
		return Trank;
	}

	uint removed_x(){
		return Tx;
	}

	void frequency(uint l,uint r,uint *ans){
		uint i,pos,pos2,*p,*q,ib,ib2,np,nq,d,d2;
		uint **P=ps;
		uint **Q=qs;
		const uint s=SH;
		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) Q[nq++]=b[ib2].c+(pos2-1)*256;
			else Q[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) P[np++]=b[ib].c+(pos-1)*256;
			else P[np]=b[ib].c;
			ib=b[ib].a[pos];

			if (ib==ib2 && pos) --np,--nq;
		}

		if (np){
			if (q[pos2]<l+l){
				if (pos2) Q[nq-1]+=256;
				else ++nq;
				d2=1;
			}

			if (p[pos]<r+r){
				if (pos) P[np-1]+=256;
				else ++np;
				d=1;
			}

			freq(P,Q,np*10+nq,ans);

			if (d==0){
				for (uint *u=bl[ib].a;r;++u,--r) ++ans[(*u>>s)&255];
			}else{
				i=bl[ib].n-r;
				for (uint *u=bl[ib].a+r;i;++u,--i) --ans[(*u>>s)&255];
			}

			if (d2==0){
				for (uint *u=bl[ib2].a;l;++u,--l) --ans[(*u>>s)&255];
			}else{
				i=bl[ib2].n-l;
				for (uint *u=bl[ib2].a+l;i;++u,--i) ++ans[(*u>>s)&255];
			}
		}else{
			memset(ans,0,1024);
			i=r-l;
			for (uint *u=bl[ib2].a+l;i;++u,--i) ++ans[(*u>>s)&255];

		}
	}

	void frequency1(uint r,uint *ans,uint i0=0,uint i1=256){
		uint i,pos,*p,ib,np,d;
		uint **P=ps;
		const uint s=SH;
		ib=root;
		np=d=0;

		for (i=hi;1<i;--i){
			p=b[ib].b;
			for (pos=0;p[pos]<r;++pos) r-=p[pos];
			if (pos) P[np++]=b[ib].c+(pos-1)*256;
			else P[np]=b[ib].c;
			ib=b[ib].a[pos];
		}

		if (np){
			if (p[pos]<r+r){
				if (pos) P[np-1]+=256;
				else ++np;
				d=1;
			}

			freq1(P,np,ans,i0,i1);

			if (d==0){
				for (uint *u=bl[ib].a;r;++u,--r) ++ans[(*u>>s)&255];
			}else{
				i=bl[ib].n-r;
				for (uint *u=bl[ib].a+r;i;++u,--i) --ans[(*u>>s)&255];
			}
		}else{
			memset(ans,0,1024);
			for (uint *u=bl[ib].a;r;++u,--r) ++ans[(*u>>s)&255];

		}
	}

	uint ranklt(uchar x){
		uint i,s=0;
		
		if (1<hi){
			uint *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;
			uint *q=bl[root].a;
			for (i=0;i<bl[root].n;++i){
				if (((q[i]>>SH)&255)<x) ++s;
			}
		}
		return s;
	}

	uint rank(uint r,uchar x){
		uint i,*p,ib=root;
		uint pos=UINT_MAX;
		uint s=0;
		const uint y=(uint)x<<SH;
		const uint m=MS;

		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 (pos!=UINT_MAX && !*p) return s;

		if (r+r<=bl[ib].n || pos==UINT_MAX){
			for (uint *q=bl[ib].a;r;++q,--r){
				if ((*q&m)==y) ++s;
			}
		}else{
			s+=pos? *p-p[-256]:*p;
			i=bl[ib].n-r;
			for (uint *q=bl[ib].a+r;i;++q,--i){
				if ((*q&m)==y) --s;
			}
		}
		return s;
	}

	uint slct_s3(uint k,uint n_,uchar x,uint *ansp,uint *ansv){
		uint i,ib,lk,uk,pos,h,*pb,s,t,r,*p,k_,*w=W;
		const uint n_0=n_;
		const uint u=(uint)x<<SH;
		const uint m=MS;

		h=hi;
		w[0]=ib=root;
		w[1]=lk=0;
		w[2]=uk=UINT_MAX;
		s=pos=0;

		for (uint pop=1;n_;pop=1){
			if (1<h){
				if (k<uk){
					p=b[ib].c+(pos*256+x);
					pb=b[ib].b;
					for (r=b[ib].n;pos=uk) break;
						}
					}
				}
			}

			if (pop){
				++h;
				w-=5;
				ib=w[0];lk=w[1];uk=w[2];s=w[3];pos=w[4];
			}
		}
		return n_0-n_;
	}

	uint c_include(uint *pc,uint pos,uint xs,uint xr){
		if (pos){
			for (uint *p=pc-256;xs<xr;++xs){
				if (p[xs]<pc[xs]) return 1;
			}
		}else{
			for (;xs<xr;++xs){
				if (pc[xs]) return 1;
			}
		}
		return 0;
	}

	uint rangereport4(const uint ql,const uint qr,const uint xs,const uint xr,uint *ansp,uint *ansv,const uint x_0,const uint x_1,uint cmp){
		if (ql>=qr || xs>=xr) return 0;
		uint i,ib,pos,h,*pb,*pc,r,*p,nl,nr,*w=W;
		const uint h_=h=hi;
		const uint *ansp_=ansp;

		w[0]=ib=root;
		pos=nl=0;

		for (uint pop=1;1;pop=1){
			if (1<h){
				pb=b[ib].b;
				pc=b[ib].c+pos*256;
				for (r=b[ib].n;pos=qr || ql>=(nl+pb[pos]) || !c_include(pc,pos,xs,xr));pc+=256,++pos) nl+=pb[pos];

				if (pos<r){
					w[1]=nl+pb[pos]; // =nr
					w[2]=pos+1;
					w+=3;
					w[0]=ib=b[ib].a[pos];
					pos=0;
					pop=0;
					--h;
				}
			}else{
				i=(ql<=nl)? 0:ql-nl;
				nr=nl+bl[ib].n;
				r=(nr<qr)? bl[ib].n:qr-nl;
				p=bl[ib].a;

				if (cmp==1){
					for (;i<r;++i){
						if (x_0<=p[i]){
							*ansv=p[i];++ansv;
							*ansp=nl+i;++ansp;
						}
					}
				}else if (cmp==2){
					for (;i<r;++i){
						if (p[i]<x_1){
							*ansv=p[i];++ansv;
							*ansp=nl+i;++ansp;
						}
					}
				}else{
					for (;i<r;++i){
						if (x_0<=p[i] && p[i]<x_1){
							*ansv=p[i];++ansv;
							*ansp=nl+i;++ansp;
						}
					}
				}
			}

			if (pop){
				if (h==h_) break;
				w-=3;
				ib=w[0];nl=w[1];pos=w[2];
				++h;
			}
		}

		return ansp-ansp_;
	}

	void sequence(uint l,uint r,uint *ans){
		uint *p,pos,ib=root;
		uint t,u,k=l;

		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
		memcpy(ans,bl[ib].a+k,t*szx);

		for (u-=t;u;u-=t){
			ib=bl[ib].pleaf;
			ans+=t;
			t=(u<=bl[ib].n)? u:bl[ib].n;
			memcpy(ans,bl[ib].a,t*szx);
		}
	}

	void getcsum(uint *ans){
		if (1<hi){
			uint i,s,*p=b[root].c+(b[root].n-1)*256;
			for (i=s=0;i<256;++i){
				ans[i]=s;
				s+=p[i];
			}
		}else{
			c_count_leaf(root,bl[root].n,ans,1);
			memmove(ans+1,ans,255*4);
			ans[0]=0;
			for (uint i=2;i<256;++i) ans[i]+=ans[i-1];
		}
	}

	void leaf_sequence(uint ib,uint k,uint n_,uint *ans){
		uint t=bl[ib].n-k;
	
		t=(n_<=t)? n_:t;
		memcpy(ans,bl[ib].a+k,t*szx);

		for (n_-=t;n_;n_-=t){
			ib=bl[ib].pleaf;
			ans+=t;
			t=(n_<=bl[ib].n)? n_:bl[ib].n;
			memcpy(ans,bl[ib].a,t*szx);
		}
	}

	void seqs3(uint *ls,uint *ns,uint na,uint *ans){
		uint ib,d,dnext,pos,h,*pb,k,t,*w=W;
		const uint h_=h=hi;

		w[0]=ib=root;
		w[1]=d=0;
		w[2]=dnext=UINT_MAX;
		pos=0;

		for (uint pop=1;na;pop=1){
			if (1<h){
				if (*ls<dnext){
					pb=b[ib].b;
					for (k=*ls;d+pb[pos]<=k;++pos) d+=pb[pos];
					w[1]=d;
					w[3]=pos;
					
					w+=4;
					w[0]=ib=b[ib].a[pos];
					w[2]=dnext=d+pb[pos];
					pos=0;
					pop=0;
					--h;
				}
			}else{
				k=*ls-d;
				t=*ns;
				leaf_sequence(ib,k,t,ans);
				ans+=t;
				++ls;++ns;
				--na;
			}

			if (pop && h<h_){
				++h;
				w-=4;
				ib=w[0];d=w[1];dnext=w[2];pos=w[3];
			}
		}
	}

	uint rep_ls(uint l,uint r,const uint x_0,const uint x_1,uint *ans){
		uint *p,pos,ib=root;
		uint s,t,u,k=l;

		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		for (s=0;u;ib=bl[ib].pleaf){
			t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
			p=bl[ib].a+k;k=0;
			for (u-=t;t;--t,++p){
				if (x_0<=*p && *p<x_1) ans[s++]=*p;
			}
		}
		return s;
	}

	uint rep_ls(uint l,uint r,const uint x_0,const uint x_1,uint *ansp,uint *ansv){
		uint *p,*q,pos,ib=root;
		uint s,t,u,k=l;

		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		for (s=0;u;ib=bl[ib].pleaf){
			t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
			p=bl[ib].a+k;k=0;
			u-=t;
			for (q=p+t;p<q;++p,++l){
				if (x_0<=*p && *p<x_1){
					ansv[s]=*p;
					ansp[s]=l;
					++s;
				}
			}
		}
		return s;
	}

	uint repme_ls(uint l,uint r,const uint x,uint *ans){
		uint *p,pos,ib=root;
		uint s,t,u,k=l;
	
		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		for (s=0;u;ib=bl[ib].pleaf){
			t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
			p=bl[ib].a+k;k=0;
			for (u-=t;t;--t,++p){
				if (x<=*p) ans[s++]=*p;
			}
		}
		return s;
	}

	uint repme_ls(uint l,uint r,const uint x,uint *ansp,uint *ansv){
		uint *p,*q,pos,ib=root;
		uint s,t,u,k=l;

		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		for (s=0;u;ib=bl[ib].pleaf){
			t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
			p=bl[ib].a+k;k=0;
			u-=t;
			for (q=p+t;p<q;++p,++l){
				if (x<=*p){
					ansv[s]=*p;
					ansp[s]=l;
					++s;
				}
			}
		}
		return s;
	}

	uint replt_ls(uint l,uint r,const uint x,uint *ans){
		uint *p,pos,ib=root;
		uint s,t,u,k=l;
	
		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		for (s=0;u;ib=bl[ib].pleaf){
			t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
			p=bl[ib].a+k;k=0;
			for (u-=t;t;--t,++p){
				if (*p<x) ans[s++]=*p;
			}
		}
		return s;
	}

	uint replt_ls(uint l,uint r,const uint x,uint *ansp,uint *ansv){
		uint *p,*q,pos,ib=root;
		uint s,t,u,k=l;

		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		for (s=0;u;ib=bl[ib].pleaf){
			t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
			p=bl[ib].a+k;k=0;
			u-=t;
			for (q=p+t;p<q;++p,++l){
				if (*p<x){
					ansv[s]=*p;
					ansp[s]=l;
					++s;
				}
			}
		}
		return s;
	}

	uint rank_ls(uint l,uint r,const uint x){
		uint *p,pos,ib=root;
		uint s,t,u,k=l;
	
		for (uint h=hi;1<h;--h){
			p=b[ib].b;
			for (pos=0;k>=p[pos];++pos) k-=p[pos];
			ib=b[ib].a[pos];
		}

		u=r-l;
		for (s=0;u;ib=bl[ib].pleaf){
			t=(u<=bl[ib].n-k)? u:bl[ib].n-k;
			p=bl[ib].a+k;k=0;
			for (u-=t;t;--t,++p){
				if (*p==x) ++s;
			}
		}
		return s;
	}
};



class wm_report_minkmaxk{
private:
	uint ntree,SH,W[256*3*3+1024*2];
	Btree256 *B;
	node *b;
	leaf_node *bl;

	uint repme(uint it,uint l,uint r,uint x,uint sh,uint *w0,uint *ans){
		uint i,j,l2,r2,t,*w1,*wc,ret=0;
		const uint d=(x>>sh)&255;

		if (r-l<=2048){
			ret=B[it].repme_ls(l,r,x,ans);

		}else if (sh){
			w1=w0+256;wc=w1+256;
			B[it].frequency1(l,w0,d,256);
			B[it].frequency1(r,w1,d,256);
			B[it].getcsum(wc);

			i=wc[d];
			l2=i+w0[d];
			r2=i+w1[d];

			j=0;			
			for (i=d+1;i<256;++i){
				if (w0[i]<w1[i]){
					w1[j]=w1[i]-w0[i];
					w0[j]=w0[i]+wc[i];
					ret+=w1[j];
					++j;
				}
			}

			t=(l2<r2)? repme(it-1,l2,r2,x,sh-8,wc,ans):0;
			if (j) B[it-1].seqs3(w0,w1,j,ans+t);
			ret+=t;

		}else{
			B[it].frequency(l,r,w0);
			const uint u=x&0xffffff00;
			for (i=d;i<256;++i){
				for (t=w0[i];t;--t)	ans[ret++]=u|i;
			}
		}
		return ret;
	}

	uint replt(uint it,uint l,uint r,uint x,uint sh,uint *w0,uint *ans){
		uint i,j,l2,r2,t,*w1,*wc,ret=0;
		const uint d=(x>>sh)&255;

		if (r-l<=2048){
			ret=B[it].replt_ls(l,r,x,ans);

		}else if (sh){
			w1=w0+256;wc=w1+256;
			B[it].frequency1(l,w0,0,d+1);
			B[it].frequency1(r,w1,0,d+1);
			B[it].getcsum(wc);

			i=wc[d];
			l2=i+w0[d];
			r2=i+w1[d];

			j=0;			
			for (i=0;i<d;++i){
				if (w0[i]<w1[i]){
					w1[j]=w1[i]-w0[i];
					w0[j]=w0[i]+wc[i];
					ret+=w1[j];
					++j;
				}
			}
			if (j) B[it-1].seqs3(w0,w1,j,ans);
			if (l2<r2) ret+=replt(it-1,l2,r2,x,sh-8,w0,ans+ret);

		}else{
			B[it].frequency(l,r,w0);
			const uint u=x&0xffffff00;
			for (i=0;i<d;++i){
				for (t=w0[i];t;--t)	ans[ret++]=u|i;
			}
		}
		return ret;
	}

	uint rngrep2(uint it,uint l,uint r,uint x_0,uint x_1,uint sh,uint *w0,uint *ans){
		uint i,j,l2,r2,l3,r3,t,*w1,*wc,ret=0;
		const uint xs=(x_0>>sh)&255;
		const uint xt=(x_1>>sh)&255;

		if (r-l<=2048){
			ret=B[it].rep_ls(l,r,x_0,x_1,ans);

		}else if (sh==0){
			B[it].frequency(l,r,w0);
			const uint u=x_0&0xffffff00;
			for (i=xs;i<xt;++i){
				for (t=w0[i];t;--t)	ans[ret++]=u|i;
			}

		}else if (xs==xt){
			l2=B[it].rank(l,xs);
			r2=B[it].rank(r,xs);
			if (l2<r2){
				i=B[it].ranklt(xs);
				ret=rngrep2(it-1,i+l2,i+r2,x_0,x_1,sh-8,w0,ans);
			}

		}else{
			w1=w0+256;wc=w1+256;
			B[it].frequency1(l,w0,xs,xt+1);
			B[it].frequency1(r,w1,xs,xt+1);
			B[it].getcsum(wc);

			i=wc[xs];
			l2=i+w0[xs];
			r2=i+w1[xs];

			i=wc[xt];
			l3=i+w0[xt];
			r3=i+w1[xt];
		
			j=0;
			for (i=xs+1;i<xt;++i){
				if (w0[i]<w1[i]){
					w1[j]=w1[i]-w0[i];
					w0[j]=w0[i]+wc[i];
					ret+=w1[j];
					++j;
				}
			}
			t=(l2<r2)? repme(it-1,l2,r2,x_0,sh-8,wc,ans):0;
			if (j) B[it-1].seqs3(w0,w1,j,ans+t);
			ret+=t;
			if (l3<r3) ret+=replt(it-1,l3,r3,x_1,sh-8,w0,ans+ret);
		}
		return ret;
	}

	void rngmnk2(uint it,uint l,uint r,uint k,const uint u,uint sh,uint *ans,uint *w0){
		uint i,j,t,*w1,*wc,k_=k;

		if (r-l<=1024){
			B[it].sequence(l,r,w0);
			qs_k(w0,k-1,r-l);
			memcpy(ans,w0,k*4);

		}else if (sh){
			w1=w0+256;wc=w1+256;
			B[it].frequency1(l,w0);
			B[it].frequency1(r,w1);
			B[it].getcsum(wc);

			for (i=j=0;k;++i){
				if (w0[i]<w1[i]){
					w1[j]=w1[i]-w0[i];
					w0[j]=w0[i]+wc[i];
					if (k<w1[j]) break;
					k-=w1[j];
					++j;
				}
			}
			if (j) B[it-1].seqs3(w0,w1,j,ans);
			if (k) rngmnk2(it-1,w0[j],w0[j]+w1[j],k,u|(i<<sh),sh-8,ans+(k_-k),w0);

		}else{
			B[it].frequency(l,r,w0);
			for (i=j=0;k;++i){
				for (t=w0[i];t && k;--t,--k) ans[j++]=u|i;
			}
		}
	}

	void select_nj_s(uint it,uint n_,uint x,uint sh,uint *ansp){
		uint i,t;
		uchar y;

		for (;it<ntree;++it,sh+=8){
			y=(x>>sh)&255;
			t=B[it].ranklt(y);
			if (t) for (i=0;i<n_;++i) ansp[i]-=t;
			B[it].slct_s3(ansp,n_,y);
		}
	}

	void rngmnk_pos2(uint it,uint l,uint r,uint k,const uint u,uint sh,uint *w0,uint *ansp,uint *ansv){
		uint i,t,*w1,*ansp_=ansp;

		if (r-l<=1024){
			t=r-l;
			w1=w0+t;
			B[it].sequence(l,r,w0);
			for (i=0;i<t;++i) w1[i]=l+i;
			qs2_i_k(w0,w1,k-1,t);
			qs2(w1,w0,k,w1+t);
			memcpy(ansv,w0,k*4);
			memcpy(ansp,w1,k*4);
			select_nj_s(it+1,k,u,sh+8,ansp);

		}else{
			w1=w0+256;
			B[it].frequency1(l,w0);
			B[it].frequency1(r,w1);

			for (i=0;k;++i){
				if (w0[i]<w1[i]){
					t=w1[i]-w0[i];
					if (k<t) break;
					k-=t;
					ansp+=t;
				}
			}

			if (ansp!=ansp_){
				t=ansp-ansp_;
				B[it].rangereport4(l,r,0,i,ansp_,ansv,0,u|(i<<sh),2);
				ansv+=t;
				select_nj_s(it+1,t,u,sh+8,ansp_);
			}

			if (k){
				if (sh){
					t=B[it].ranklt(i);
					uint l2=t+w0[i];
					uint r2=t+w1[i];
					rngmnk_pos2(it-1,l2,r2,k,u|(i<<sh),sh-8,w0,ansp,ansv);
				}else{
					B[it].slct_s3(w0[i],k,i,ansp,ansv);
					select_nj_s(it+1,k,u,sh+8,ansp);
				}
			}
		}
	}

	void rngmxk2(uint it,uint l,uint r,uint k,const uint u,uint sh,uint *ans,uint *w0){
		uint i,j,t,*w1,*wc;

		if (r-l<=1024){
			B[it].sequence(l,r,w0);
			qs_k(w0,(r-l)-k,r-l);
			memcpy(ans,w0+((r-l)-k),k*4);
			
		}else if (sh){
			w1=w0+256;wc=w1+256;
			B[it].frequency1(l,w0);
			B[it].frequency1(r,w1);
			B[it].getcsum(wc);

			for (i=j=255;k;--i){
				if (w0[i]<w1[i]){
					w1[j]=w1[i]-w0[i];
					w0[j]=w0[i]+wc[i];
					if (k<w1[j]) break;
					k-=w1[j];
					--j;
				}
			}
			if (j!=255) B[it-1].seqs3(w0+(j+1),w1+(j+1),255-j,ans+k);
			if (k) rngmxk2(it-1,w0[j],w0[j]+w1[j],k,u|(i<<sh),sh-8,ans,w0);

		}else{
			B[it].frequency(l,r,w0);
			for (i=255;k;--i){
				for (t=w0[i];t && k;--t) ans[--k]=u|i;
			}
		}
	}

	void rngmxk_pos2(uint it,uint l,uint r,uint k,const uint u,uint sh,uint *w0,uint *ansp,uint *ansv){
		uint i,t,*w1,*ansp_=ansp;

		if (r-l<=1024){
			t=r-l;
			w1=w0+t;
			B[it].sequence(l,r,w0);
			for (i=0;i<t;++i) w1[i]=l+i;
			qs2_i_k(w0,w1,t-k,t);
			w0+=t-k;w1+=t-k;
			qs2(w1,w0,k,w1+t);
			memcpy(ansv,w0,k*4);
			memcpy(ansp,w1,k*4);
			select_nj_s(it+1,k,u,sh+8,ansp);

		}else{
			w1=w0+256;
			B[it].frequency1(l,w0);
			B[it].frequency1(r,w1);

			for (i=255;k;--i){
				if (w0[i]<w1[i]){
					t=w1[i]-w0[i];
					if (k<t) break;
					k-=t;
					ansp+=t;
				}
			}

			if (ansp!=ansp_){
				t=ansp-ansp_;
				B[it].rangereport4(l,r,i+1,256,ansp_,ansv,u|((i+1)<<sh),0,1);
				ansv+=t;
				select_nj_s(it+1,t,u,sh+8,ansp_);
			}

			if (k){
				if (sh){
					t=B[it].ranklt(i);
					uint l2=t+w0[i];
					uint r2=t+w1[i];
					rngmxk_pos2(it-1,l2,r2,k,u|(i<<sh),sh-8,w0,ansp,ansv);
				}else{
					B[it].slct_s3(w0[i],k,i,ansp,ansv);
					select_nj_s(it+1,k,u,sh+8,ansp);
				}
			}
		}
	}

	uint repme_pos(uint it,uint l,uint r,uint x,uint sh,uint *ansp,uint *ansv){
		uint i,l2,r2,t,ret=0;
		const uint d=(x>>sh)&255;

		if (r-l<=1024){
			ret=B[it].repme_ls(l,r,x,ansp,ansv);
			if (ret) select_nj_s(it+1,ret,x,sh+8,ansp);

		}else if (sh){
			l2=B[it].rank(l,d);
			r2=B[it].rank(r,d);
			if (l2<r2){
				i=B[it].ranklt(d);
				ret=repme_pos(it-1,i+l2,i+r2,x,sh-8,ansp,ansv);
				ansp+=ret;ansv+=ret;
			}
			t=B[it].rangereport4(l,r,d+1,256,ansp,ansv,((x>>sh)+1)<<sh,0,1);
			if (t) select_nj_s(it+1,t,x,sh+8,ansp);
			ret+=t;

		}else{
			t=B[it].rangereport4(l,r,d,256,ansp,ansv,x,0,1);
			if (t) select_nj_s(it+1,t,x,sh+8,ansp);
			ret+=t;
		}
		return ret;
	}

	uint replt_pos(uint it,uint l,uint r,uint x,uint sh,uint *ansp,uint *ansv){
		uint i,l2,r2,ret=0;
		const uint d=(x>>sh)&255;

		if (r-l<=1024){
			ret=B[it].replt_ls(l,r,x,ansp,ansv);
			if (ret) select_nj_s(it+1,ret,x,sh+8,ansp);

		}else{
			ret=B[it].rangereport4(l,r,0,d,ansp,ansv,0,(x>>sh)<<sh,2);
			if (ret) select_nj_s(it+1,ret,x,sh+8,ansp);
			if (sh){
				l2=B[it].rank(l,d);
				r2=B[it].rank(r,d);
				if (l2<r2){
					i=B[it].ranklt(d);
					ret+=replt_pos(it-1,i+l2,i+r2,x,sh-8,ansp+ret,ansv+ret);
				}
			}
		}
		return ret;
	}

	uint rngrep_pos2(uint it,uint l,uint r,uint x_0,uint x_1,uint sh,uint *ansp,uint *ansv){
		uint i,l2,r2,l3,r3,t,ret=0;
		const uint xs=(x_0>>sh)&255;
		const uint xt=(x_1>>sh)&255;

		if (r-l<=1024){
			ret=B[it].rep_ls(l,r,x_0,x_1,ansp,ansv);
			if (ret) select_nj_s(it+1,ret,x_0,sh+8,ansp);

		}else if (sh==0){
			ret=B[it].rangereport4(l,r,xs,xt,ansp,ansv,x_0,x_1,3);
			if (ret) select_nj_s(it+1,ret,x_0,sh+8,ansp);

		}else if (xs==xt){
			l2=B[it].rank(l,xs);
			r2=B[it].rank(r,xs);
			if (l2<r2){
				i=B[it].ranklt(xs);
				ret=rngrep_pos2(it-1,i+l2,i+r2,x_0,x_1,sh-8,ansp,ansv);
			}
		}else{

			l2=B[it].rank(l,xs);
			r2=B[it].rank(r,xs);
			if (l2<r2){
				i=B[it].ranklt(xs);
				ret=repme_pos(it-1,i+l2,i+r2,x_0,sh-8,ansp,ansv);
				ansp+=ret;ansv+=ret;
			}

			if (xs+1<xt){
				t=B[it].rangereport4(l,r,xs+1,xt,ansp,ansv,((x_0>>sh)+1)<<sh,(x_1>>sh)<<sh,3);
				if (t) select_nj_s(it+1,t,x_0,sh+8,ansp);
				ansp+=t;ansv+=t;ret+=t;
			}

			l3=B[it].rank(l,xt);
			r3=B[it].rank(r,xt);
			if (l3<r3){
				i=B[it].ranklt(xt);
				ret+=replt_pos(it-1,i+l3,i+r3,x_1,sh-8,ansp,ansv);
			}
		}
		return ret;
	}

public:

	// n:最大入力件数
	// m:最大入力値
	wm_report_minkmaxk(uint n,uint m){
		uint i,j,ni,nl;
		double c=0.4;

		ntree=0;
		for (i=m;i;i>>=8) ++ntree;
		SH=(ntree-1)*8;

		nl=(double)n/(c*binfo.leafsize);
		ni=((double)nl-1)/(c*binfo.M-1);
		++nl,++ni;

		B=new Btree256[ntree];
		b=new node[ni*ntree];
		bl=new leaf_node[nl*ntree];

		for (i=0,j=ntree-1;i<ntree;++i,--j){
			B[j].init_mem(b+i*ni,bl+i*nl,nl,j*8);
		}
		memset(W,0,256*4);	
	}

	~wm_report_minkmaxk(){
		delete [] B;delete [] b;delete [] bl;
	}

	// 初期化時に数列の最後尾にxを追加する
	void init(uint x){
		B[ntree-1].init(x);
		++W[x>>SH];
	}

	// init後にこれを実行すると初期化完了
	void init_end(){
		uint i,j,k,n_,x,ix,*wc,*w1,s=SH;
		wc=W;w1=W+256;ix=0;
		for (i=n_=0;i<256;++i) n_+=wc[i];

		for (i=ntree-1;0<i;--i,s-=8){
			for (j=1;j<256;++j) wc[j]+=wc[j-1];
			memset(w1,0,256*4);
			for (j=n_;j--;){
				x=B[i].init_access_x(j);
				k=--wc[(x>>s)&255];
				B[i-1].init2(k,x);
				++w1[(x>>(s-8))&255];
			}
			memcpy(wc,w1,256*4);
		}

		for (i=ntree;i--;) B[i].init_end();
	}

	// i番目の要素を返す
	uint access(uint i){
		return B[ntree-1].access(i);
	}

	// i番目にxを追加
	void insert(uint i,uint x){
		uint j,k,s=SH;
		uchar y;
	
		for (j=ntree;j--;s-=8){
			y=(x>>s)&255;
			k=B[j].insert_rank(i,x,j);
			if (j) i=B[j].ranklt(y)+k;
		}
	}

	// i番目の要素を削除
	uint remove(uint i){
		uint j,k,x,s=SH;
		uint y=256;

		for (j=ntree;j--;s-=8){
			k=B[j].remove_rank(i,j);
			if (y==256) x=B[j].removed_x();
			y=(x>>s)&255;
			if (j) i=B[j].ranklt(y)+k;
		}
		return x;
	}

	// [l,r)に含まれるxの個数
	uint rank(uint l,uint r,uint x){
		uint i,j,y,s=SH;
		if (l>=r) return 0;

		for (i=ntree;i-- && l<r;s-=8){
			if (255<r-l){
				y=(x>>s)&255;
				j=i? B[i].ranklt(y):0;
				l=j+B[i].rank(l,y);
				r=j+B[i].rank(r,y);
			}else{
				return B[i].rank_ls(l,r,x);
			}
		}
		return r-l;
	}

	// [l,r)に含まれるx以上y未満の要素をans配列に返す
	// rangereportは解の個数を返す
	uint rangereport(uint l,uint r,uint x,uint y,uint *ans){
		if (l>=r || x>=y) return 0;
		return rngrep2(ntree-1,l,r,x,y,SH,W,ans);
	}

	// [l,r)に含まれる要素のうち小さい方からk個をans配列に返す
	void rangemink(uint l,uint r,uint k,uint *ans){
		if (l>=r || k>r-l) return;
		rngmnk2(ntree-1,l,r,k,0,SH,ans,W);
	}

	// [l,r)に含まれる要素のうち小さい方からk個をans_val配列に返す
	// ans_pos[i]=ans_val[i]の位置
	void rangemink_pos(uint l,uint r,uint k,uint *ans_pos,uint *ans_val){
		if (l>=r || k>r-l || !k) return;
		rngmnk_pos2(ntree-1,l,r,k,0,SH,W,ans_pos,ans_val);
	}

	// [l,r)に含まれる要素のうち大きい方からk個をans配列に返す
	void rangemaxk(uint l,uint r,uint k,uint *ans){
		if (l>=r || k>r-l || !k) return;
		rngmxk2(ntree-1,l,r,k,0,SH,ans,W);
	}

	// [l,r)に含まれる要素のうち大きい方からk個をans_val配列に返す
	// ans_pos[i]=ans_val[i]の位置
	void rangemaxk_pos(uint l,uint r,uint k,uint *ans_pos,uint *ans_val){
		if (l>=r || k>r-l || !k) return;
		rngmxk_pos2(ntree-1,l,r,k,0,SH,W,ans_pos,ans_val);
	}

	// [l,r)に含まれるx以上y未満の要素をans_val配列に返す
	// ans_pos[i]=ans_val[i]の位置
	// rangereport_posは解の個数を返す
	uint rangereport_pos(uint l,uint r,uint x,uint y,uint *ans_pos,uint *ans_val){
		if (l>=r || x>=y) return 0;
		return rngrep_pos2(ntree-1,l,r,x,y,SH,ans_pos,ans_val);
	}
};



int main(void){
	const uint N=1000000;      // 最大入力件数
	const uint M=UINT_MAX;     // 最大入力値
	const uint n=100000;       // クエリの回数(大)
	const uint n_small=10000;  // クエリの回数(小)
	const uint k=100;          // rangeminkのk
	
	uint i,j,t,l,r,x,y,ret;
	uint *anspos=new uint[N];
	uint *ansval=new uint[N];

	wm_report_minkmaxk wm(N,M);


	// N要素で初期化
	std::cout<<"init(input size="<<N<<")"<<std::endl;
	for (i=0;i<N;++i){
		x= (double)rand() / (RAND_MAX + 1) * (M);
		wm.init(x);
	}
	wm.init_end();
	std::cout<<std::endl;std::cout<<std::endl;std::cout<<std::endl;
	

	// n要素削除
	std::cout<<"remove *"<<n<<std::endl;
	for (i=0;i<n;++i){
		j= (double)rand() / (RAND_MAX + 1) * (N-i);
		wm.remove(j);
	}
	std::cout<<std::endl;std::cout<<std::endl;std::cout<<std::endl;


	// n要素追加
	std::cout<<"insert *"<<n<<std::endl;
	for (i=n;0<i;--i){
		j= (double)rand() / (RAND_MAX + 1) * (N-i);
		x= (double)rand() / (RAND_MAX + 1) * (M);
		wm.insert(j,x);
	}
	std::cout<<std::endl;std::cout<<std::endl;std::cout<<std::endl;


	// rangemink n回
	std::cout<<"rangemink *"<<n<<std::endl;
	for (i=0;i<n;++i){
		do{
			l= (double)rand() / (RAND_MAX + 1) * (N);
			r= (double)rand() / (RAND_MAX + 1) * (N);
			if (l>r) t=l,l=r,r=t;
		}while (k>r-l);
		wm.rangemink(l,r,k,ansval);	
//		wm.rangemaxk(l,r,k,ansval);	
	}

	// 最後のクエリの結果を表示
	std::cout<<"rangemink(l,r,k)"<<std::endl;
	std::cout<<"l="<<l<<",r="<<r<<",k="<<k<<std::endl;
	std::cout<<std::endl;
	for (i=0;i<k;++i){
		std::cout<<ansval[i]<<",";
	}
	std::cout<<std::endl;std::cout<<std::endl;std::cout<<std::endl;std::cout<<std::endl;


	// rangereport n回
	// l,r,xを乱数で決めてそれに応じて解がk個程になるようにyを決める
	// 最後に解の個数の平均を表示する
	std::cout<<"rangereport *"<<n<<std::endl;
	const double Mk=(double)M*(double)k;
	unsigned long long avg=0ULL;

	for (i=0;i<n;++i){
		do{
			l= (double)rand() / (RAND_MAX + 1) * (N);
			r= (double)rand() / (RAND_MAX + 1) * (N);
			if (l>r) t=l,l=r,r=t;
		}while (k>r-l);
		do{
			x= (double)rand() / (RAND_MAX + 1) * (M);
			y=x+Mk/(double)(r-l);
		}while (x>y || y>M);
		ret=wm.rangereport(l,r,x,y,ansval);
		avg+=(unsigned long long)ret;
	}

	// 最後のクエリの結果を表示
	std::cout<<"rangereport(l,r,x,y)"<<std::endl;
	std::cout<<"l="<<l<<",r="<<r<<std::endl;
	std::cout<<"x="<<x<<",y="<<y<<std::endl;
	std::cout<<"output size="<<ret<<std::endl;
	std::cout<<std::endl;
	for (i=0;i<ret;++i){
		std::cout<<ansval[i]<<",";
	}
	std::cout<<std::endl;std::cout<<std::endl;
	std::cout<<"average output size="<<(double)avg/(double)n<<std::endl;
	std::cout<<std::endl;std::cout<<std::endl;std::cout<<std::endl;


	// 位置の取得は遅いのでクエリを1万回にする

	// rangemink_pos n_small回
	std::cout<<"rangemink_pos *"<<n_small<<std::endl;
	for (i=0;i<n_small;++i){
		do{
			l= (double)rand() / (RAND_MAX + 1) * (N);
			r= (double)rand() / (RAND_MAX + 1) * (N);
			if (l>r) t=l,l=r,r=t;
		}while (k>r-l);
		wm.rangemink_pos(l,r,k,anspos,ansval);
//		wm.rangemaxk_pos(l,r,k,anspos,ansval);
	}

	// 最後のクエリの結果を表示
	// val(i)=access(pos(i)) でなくてはいけない
	std::cout<<"rangemink_pos(l,r,k,pos,val)"<<std::endl;
	std::cout<<"l="<<l<<",r="<<r<<",k="<<k<<std::endl;
	std::cout<<std::endl;
	std::cout<<"i,pos(i),val(i),access(pos(i))"<<std::endl;
	for (i=0;i<k;++i){
		t=wm.access(anspos[i]);
		std::cout<<i<<","<<anspos[i]<<","<<ansval[i]<<","<<t<<std::endl;
	}
	std::cout<<std::endl;std::cout<<std::endl;std::cout<<std::endl;


	// rangereport_pos n_small回
	std::cout<<"rangereport_pos *"<<n_small<<std::endl;
	avg=0ULL;
	for (i=0;i<n_small;++i){
		do{
			l= (double)rand() / (RAND_MAX + 1) * (N);
			r= (double)rand() / (RAND_MAX + 1) * (N);
			if (l>r) t=l,l=r,r=t;
		}while (k>r-l);
		do{
			x= (double)rand() / (RAND_MAX + 1) * (M);
			y=x+Mk/(double)(r-l);
		}while (x>y || y>M);
		ret=wm.rangereport_pos(l,r,x,y,anspos,ansval);
		avg+=(unsigned long long)ret;
	}

	// 最後のクエリの結果を表示
	// val(i)=access(pos(i)) でなくてはいけない
	std::cout<<"rangereport_pos(l,r,x,y,pos,val)"<<std::endl;
	std::cout<<"l="<<l<<",r="<<r<<std::endl;
	std::cout<<"x="<<x<<",y="<<y<<std::endl;
	std::cout<<"output size="<<ret<<std::endl;
	std::cout<<std::endl;
	std::cout<<"i,pos(i),val(i),access(pos(i))"<<std::endl;
	for (i=0;i<ret;++i){
		t=wm.access(anspos[i]);
		std::cout<<i<<","<<anspos[i]<<","<<ansval[i]<<","<<t<<std::endl;
	}
	std::cout<<std::endl;

	std::cout<<"average output size="<<(double)avg/(double)n_small<<std::endl;
	std::cout<<std::endl;std::cout<<std::endl;


	delete [] anspos;
	delete [] ansval;
	return 0;
}

 







  

静的区間転倒数2

領域O(n)、構築O(n*sqrt(n/log(n)))、クエリO(sqrt(n*log(n)))、を示す

数列を大バケットと小バケットの2種類のバケットで分割する

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;
}