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