動的ウェーブレット行列

ウェーブレット行列を動的ビットベクトルを使って動的にしてある。

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