静的区間転倒数2
領域O(n)、構築O(n*sqrt(n/log(n)))、クエリO(sqrt(n*log(n)))、を示す
n : 数列長
B : 大バケットのサイズ
b : 小バケットのサイズ
n=Bb
大バケットがb個ある、小バケットがB個ある、b≦sqrt(n)≦B、となる
(Bとbは整数で、Bはbの整数倍とする)
(必要に応じて番兵を入れてnを調整する)
WM(ウェーブレット行列)を使うので数列は0以上n未満の整数値に変換する
順列に変換しても問題ないので変換後の数列は順列とする
I(X) : 数列の区間Xの転倒数
I(X,Y) : 集合{(x,y)|x>y,x∈X,y∈Y}の要素数
図のX' : クエリ区間で大バケットの左にはみ出した小バケットちょうど何個分の区間
図のY : クエリ区間で大バケットちょうど何個分の区間
I(X',Y)を前処理で全て計算してもn=Bbなので領域O(n)で保存できる(I(Y,Z')も同様)
図のX : クエリ区間で小バケットの左にはみ出した区間
I(X,Y)はWMで1要素ずつ計算する(I(Y,Z)も同様)
クエリ区間[L,R)が大バケットを1つ以上含む場合の転倒数
=I(XX')+I(Z'Z)+I(Y)+I(X',Y)+I(Y,Z')+I(X,Y)+I(Y,Z)+I(XX',Z'Z)
各項の時間計算量
I(XX') : 前処理によりO(1)
I(Z'Z) : 同上
I(Y) : 前処理によりO(1)
I(X',Y) : 前処理によりO(1)
I(Y,Z') : 同上
I(X,Y) : WMによりO(b*log(n))
I(Y,Z) : 同上
I(XX',Z'Z) : 大バケットを(値,添字)の値をキーにしてソートした列を使う尺取り法によりO(B)
クエリ区間が大バケットを含まず大バケット間の境界を跨ぐ場合の区間をXX'Z'Zとする(図からYを取り除いたもの)
I(XX'Z'Z)=I(XX')+I(Z'Z)+I(XX',Z'Z)
計算量はO(B)
クエリ区間が1つの大バケットに含まれる場合の区間をX(図のX)とする
I(X)=I(XX')-I(X')-I(X,X')
計算量はO(B)
以上によりクエリの計算量はO(max(b*log(n),B))となる
B=sqrt(n*log(n))、b=sqrt(n/log(n))、とすると冒頭の計算量を達成できる
Xをある大バケット、YをX以外の大(小)バケットとするとき、I(X,Y)やI(Y,X)の構築は数列の要素がn未満の整数値なので次のようにXに関する累積和みたいなものを作り、Xに含まれるx未満(以上)の要素数をO(1)で求められるようにするのがおそらく一番簡単である
サイズnの作業領域wを0で初期化
for each x in X{
w[x]+=1
}
for (i=1;i<n;++i){
w[i]+=w[i-1]
}このI(X,Y)やI(Y,X)をもとに、図のI(Y)やI(X',Y)やI(Y,Z')を構築する
変種1
n=Bbの制約を外して、B=sqrt(n)、b=sqrt(n)/log(n)、とすると領域O(n*log(n))、クエリO(sqrt(n))になる
変種2
前に「平方分割でrankltはO(sqrt(n))でできる」ことを書いたが、この平方分割をバケット法として一般化して考えると「バケット法でrankltはクエリ区間が大バケットちょうど何個分ならばO(b)でできる(bはn=Bbのb)」ということがいえる
さらに、このO(b)は小バケット内の全探索なのでWMを使うとここがO(log(b))になるので「バケット法とWMでrankltはクエリ区間が大バケットちょうど何個分ならばO(log(b))でできる」になる
これを使うと図のI(X,Y)とI(Y,Z)はO(b*log(n))からO(b*log(b))になるので、n=b*b*log(b)のとき、クエリがO(b*log(b))になる
(I(X,Y)とI(XX',Z'Z)の計算量が等しい ⇔ b*log(b)=B (B=n/b))
O(b*log(b))がnの式でどう書けるのか分からないので計算量は分からない
(下のコードでは、B=sqrt(n*log(n))、b=sqrt(n/log(n))、であるがこのバケット法を使っているのでWMの処理の部分はO(b*log(n))ではなくO(b*log(b))であり定数倍の改善になっている)
変種3
WMを使わなくてもバケット法でI(X,Y)とI(Y,Z)はO(b*b)なので、B=pow(n,2/3)、b=pow(n,1/3)、とすると領域O(n)、クエリO(pow(n,2/3))になる
#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <algorithm>
#define ll long long
#define ushort unsigned short
class RangeInversions2{
private:
struct bucket{
int *a;
ll *inv;
};
struct Bucket{
int *a,*sa,*ix,*w;
ll *inv,*li,*ri;
bucket *bk;
};
class wm16{
private:
ushort *a,bsz,nb,*tbl,db,M,nbit;
void add1(ushort *pa,ushort *ps,ushort i){
pa[i>>4]|=(ushort)1<<(i&15);
++ps[(i>>4)+1];
}
void wm_build(int *pb,int *wr,ushort *pa,ushort *ps){
int i,n0,n1,*tp,m,nt=nbit;
for (m=M;m;m>>=1){
n0=n1=nt;
for (i=0;i<nt;++i){
if (pb[i]&m){
add1(pa,ps,i);
--n0;
}
}
for (i=1;i<bsz;++i) ps[i]+=ps[i-1];
for (i=nt-1;0<=i;--i){
if (pb[i]&m) wr[--n1]=pb[i];
else wr[--n0]=pb[i];
}
tp=pb,pb=wr,wr=tp;
pa=ps+bsz;
ps=pa+bsz;
}
}
public:
int calc_mem(int nB_,int bsz_,int nb_){
int i;
nbit=bsz_;
bsz=(bsz_%16)? bsz_/16+1:bsz_/16;
++bsz;
nb=nb_*nB_;
for (i=0;(1<<i)<=nB_;++i);
M=(ushort)1<<(i-1);
db=bsz*2*i;
return db*nb+65536;
}
void init(int *pb,int bsz_,int *wr,ushort *d_){
int i,j,*p;
ushort *q;
a=d_;
tbl=a+db*nb;
p=pb;
q=a;
for (i=0;i<nb;++i){
wm_build(p,wr,q,q+bsz);
p+=bsz_;
q+=db;
}
for (i=0;i<65536;++i){
for (j=i;j;j-=j&-j) ++tbl[i];
}
}
ll rangefreqs(int *pa_,int n_,int bsz_,ushort x,ushort y,int lside){
int i,s=0;
ushort ib,t,u,m,*pa,*ps,*tb,sz,l,r,l2,r2,nt=nbit;
sz=bsz;
tb=tbl;
for (i=0;i<n_;++i){
ib=pa_[i]/bsz_;
if (lside){
r2=r=pa_[i]-ib*bsz_;
l2=l=0;
}else{
l2=l=pa_[i]-ib*bsz_+1;
r2=r=bsz_;
}
pa=a+ib*db;
ps=pa+sz;
for (m=M;m;m>>=1){
if (l2<r2){
t=tb[pa[l2>>4]&~(USHRT_MAX<<(l2&15))]+ps[l2>>4];
u=tb[pa[r2>>4]&~(USHRT_MAX<<(r2&15))]+ps[r2>>4];
if (y&m){
s+=r2-l2-(u-t);
l2=nt-ps[sz-1]+t;
r2=l2+(u-t);
}else{
l2-=t;
r2-=u;
}
}
if (l<r){
t=tb[pa[l>>4]&~(USHRT_MAX<<(l&15))]+ps[l>>4];
u=tb[pa[r>>4]&~(USHRT_MAX<<(r&15))]+ps[r>>4];
if (x&m){
s-=r-l-(u-t);
l=nt-ps[sz-1]+t;
r=l+(u-t);
}else{
l-=t;
r-=u;
}
}
pa=ps+sz;
ps=pa+sz;
}
}
return s;
}
};
int *a,N,Bsz,bsz,nB,nb,*wi,*C;
ll *Bi;
ushort *pwm;
Bucket *Bk;
bucket *bk;
wm16 wm;
void BIT_init(int *BIT,int sz){
memset(BIT,0,(sz+1)*4);
}
void BIT_add1(int *BIT,int i,int sz){
for (++i;i<=sz;i+=i&-i) ++BIT[i];
}
ll BIT_sum(int *BIT,int i){
ll s=0LL;
for (++i;0<i;i-=i&-i) s+=(ll)BIT[i];
return s;
}
void initLR_lt(int *L,int *Li,int *R,int *w){
int i,j,sz=Bsz;
memset(w,0,nb*4);
for (i=j=0;i<sz;++i){
while (j<sz && L[i]>R[j]) ++j;
w[Li[i]]+=j;
}
}
void initLR_mt(int *L,int *R,int *Ri,int *w){
int i,j,sz=Bsz;
memset(w,0,nb*4);
for (i=j=0;i<sz;++i){
while (j<sz && L[j]<=R[i]) ++j;
w[Ri[i]]+=sz-j;
}
}
void C_build(int *pb){
int j,*p,*q,*r;
p=C;
q=pb;
r=C+N;
for (;p<r;p+=bsz,q+=bsz){
for (j=0;j<bsz;++j) ++p[q[j]];
}
for (p=C+bsz;p<r;p+=bsz){
for (j=0;j<bsz;++j) p[j]+=p[j-bsz];
}
for (p=C;p<r;p+=bsz){
for (j=1;j<bsz;++j) p[j]+=p[j-1];
}
}
void Bk_build_Bi(){
int i,j;
for (i=0;i<nB;++i){
ll t=Bk[i].inv[i];
for (j=i+1;j<nB;++j) Bk[i].inv[j]=t+Bk[i].bk[0].inv[j];
}
for (i=nB-1;0<=i;--i){
for (j=i+1;j<nB;++j) Bk[i].inv[j]+=Bk[i+1].inv[j];
}
}
void Bk_build_lr(int ib,int sz,int *pa,int *pBIT){
int i;
ll *q,*r;
BIT_init(pBIT,sz);
q=Bk[ib].li;
for (i=Bsz-1;0<=i;--i){
q[i]=BIT_sum(pBIT,pa[i]-1);
BIT_add1(pBIT,pa[i],sz);
}
for (i=Bsz-1;0<i;--i) q[i-1]+=q[i];
Bk[ib].inv[ib]=q[0];
BIT_init(pBIT,sz);
r=Bk[ib].ri;
for (i=0;i<Bsz;++i){
r[i]=i-BIT_sum(pBIT,pa[i]);
BIT_add1(pBIT,pa[i],sz);
}
for (i=1;i<Bsz;++i) r[i]+=r[i-1];
}
void bk_build_Bi(){
int i,j,k;
for (i=0;i<nB;++i){
for (j=0;j<i;++j){
initLR_lt(Bk[j].sa,Bk[j].ix,Bk[i].sa,Bk[j].w);
for (k=0;k<nb;++k) Bk[j].bk[k].inv[i]=Bk[j].w[k];
}
for (j=nB-1;i<j;--j){
initLR_mt(Bk[i].sa,Bk[j].sa,Bk[j].ix,Bk[j].w);
for (k=0;k<nb;++k) Bk[j].bk[k].inv[i]=Bk[j].w[k];
}
}
for (i=0;i<nB;++i){
for (j=0;j<nb;++j){
for (k=i-1;0<k;--k) Bk[i].bk[j].inv[k-1]+=Bk[i].bk[j].inv[k];
for (k=i+1;k+1<nB;++k) Bk[i].bk[j].inv[k+1]+=Bk[i].bk[j].inv[k];
}
}
for (i=0;i<nB;++i){
for (j=1;j<nb;++j){
for (k=0;k<i;++k) Bk[i].bk[j].inv[k]+=Bk[i].bk[j-1].inv[k];
}
for (j=nb-1;0<j;--j){
for (k=i+1;k<nB;++k) Bk[i].bk[j-1].inv[k]+=Bk[i].bk[j].inv[k];
}
}
}
ll C_query(int l,int l2,int lB,int r,int r2,int rB){
int i,ib,*pa,*pc,*p,sz=bsz;
ll s;
pa=a;
pc=C;
p=pc+(nb*nB-1)*sz;
s=(ll)((p[rB-1]-p[lB])*(r2-r));
for (i=l;i<l2;++i){
ib=pa[i]/sz;
if (ib){
p=pc+(ib-1)*sz;
s+=(ll)(p[rB-1]-p[lB]);
}
}
for (i=r;i<r2;++i){
ib=pa[i]/sz;
p=pc+ib*sz;
s-=(ll)(p[rB-1]-p[lB]);
}
return s;
}
static void qs2(int *a,int *a2,int argn){
int l,r,m,m2,t;
l=0,r=argn-1;
m=a[r>>1];
m2=a2[r>>1];
while (l<=r){
while ((a[l]<m) || (a[l]==m && a2[l]<m2)) ++l;
while ((a[r]>m) || (a[r]==m && a2[r]>m2)) --r;
if (l<=r){
t=a[l],a[l]=a[r],a[r]=t;
t=a2[l],a2[l]=a2[r],a2[r]=t;
++l,--r;
}
}
if (1<l) qs2(a,a2,l);
if (l+1<argn) qs2(a+l,a2+l,argn-l);
}
ll invLR(const int l,const int r,int lB,int rB){
int i,j,k,l2,r2,u,v,*pl,*pr,*il,*ir;
ll s=0LL;
pl=Bk[lB].sa;
il=Bk[lB].ix;
pr=Bk[rB].sa;
ir=Bk[rB].ix;
l2=l-lB*Bsz;
r2=r-rB*Bsz;
u=Bsz-l2;
v=r2;
for (i=j=0;u;++i){
if (l2<=il[i]){
for (k=pl[i];v && pr[j]<k;++j){
if (ir[j]<r2) --v;
}
s+=(ll)(r2-v);
--u;
}
}
return s;
}
ll invself(int ib,int l,int r){
int i,*ix,sz=Bsz;
ll nr,s=0LL;
nr=s;
ix=Bk[ib].ix;
l%=sz;
r=(r%sz)? r%sz:sz;
for (i=0;i<sz;++i){
if (l<=ix[i]){
if (r<=ix[i]) ++nr;
else s+=nr;
}
}
return s;
}
public:
// 数列A,数列長n(1<n)
RangeInversions2(int const *A,const int n){
int i,j,k,*p,*q;
ll *pl;
const double d=log((double)n)/log(2.0);
const int mBIT=1;
nB=bsz=sqrt((double)n/d);
for (nb=d;bsz*nb*nB<n;++nb);
Bsz=bsz*nb;
N=Bsz*nB;
i=N+(Bsz*2+nb)*nB+N+(N+Bsz+mBIT);
j=((nB+Bsz*2)+(nb*nB))*nB;
k=wm.calc_mem(nB,bsz,nb);
Bk=new Bucket[nB];
bk=new bucket[nb*nB];
a=new int[i];
Bi=new ll[j];
pwm=new ushort[k];
memset(a,0,i*4);
memset(Bi,0,j*8);
memset(pwm,0,k*2);
wi=a+N;
pl=Bi;
for (i=0;i<nB;++i){
Bk[i].a=a+i*Bsz;
Bk[i].sa=wi;wi+=Bsz;
Bk[i].ix=wi;wi+=Bsz;
Bk[i].w=wi;wi+=nb;
Bk[i].inv=pl;pl+=nB;
Bk[i].li=pl;pl+=Bsz;
Bk[i].ri=pl;pl+=Bsz;
Bk[i].bk=bk+i*nb;
for (j=0;j<nb;++j){
Bk[i].bk[j].a=Bk[i].a+j*bsz;
Bk[i].bk[j].inv=pl;pl+=nB;
}
}
C=wi;wi+=N;
memcpy(a,A,n*4);
for (i=0;i<N;++i) wi[i]=i;
qs2(a,wi,n);
for (i=0;i<n;++i) a[wi[i]]=i;
for (i=n;i<N;++i) a[i]=n;
for (i=0;i<N;++i) wi[i]/=Bsz;
C_build(wi);
wm.init(wi,bsz,wi+N,pwm);
p=a;
for (i=0;i<nB;++i,p+=Bsz){
memcpy(Bk[i].sa,p,Bsz*4);
q=Bk[i].ix;
for (j=0;j<Bsz;++j) q[j]=j;
qs2(Bk[i].sa,q,Bsz);
for (j=0;j<Bsz;++j) wi[q[j]]=j;
Bk_build_lr(i,Bsz,wi,wi+Bsz);
}
for (i=0;i<nB;++i){
memcpy(wi+i*Bsz,Bk[i].ix,Bsz*4);
for (j=0;j<Bsz;++j) Bk[i].ix[j]/=bsz;
}
bk_build_Bi();
for (i=0;i<nB;++i) memcpy(Bk[i].ix,wi+i*Bsz,Bsz*4);
Bk_build_Bi();
}
~RangeInversions2(){
delete [] Bk;delete [] bk;delete [] a;delete [] Bi;delete [] pwm;
}
// Aの[l,r)の転倒数
ll inversions_query(const int l,const int r){
if (l>=r) return 0LL;
int lb,rb,lB,rB;
ll s;
lB=l/Bsz;
rB=(r-1)/Bsz;
if (lB+1<rB){
lb=l/bsz;
rb=(r-1)/bsz;
s=Bk[lB].li[l%Bsz]+Bk[rB].ri[(r-1)%Bsz]+Bk[lB+1].inv[rB-1]+invLR(l,r,lB,rB)+C_query(l,(lb+1)*bsz,lB,rb*bsz,r,rB)+wm.rangefreqs(a+l,(lb+1)*bsz-l,bsz,lB+1,rB,1)+wm.rangefreqs(a+rb*bsz,r-rb*bsz,bsz,lB+1,rB,0);
if ((lb+1)%nb) s+=bk[lb+1].inv[rB-1];
if (rb%nb) s+=bk[rb-1].inv[lB+1];
}else if (lB<rB){
s=Bk[lB].li[l%Bsz]+Bk[rB].ri[(r-1)%Bsz]+invLR(l,r,lB,rB);
}else{
s=(r%Bsz)? Bk[lB].li[l%Bsz]-Bk[lB].li[r%Bsz]-invself(lB,l,r):Bk[lB].li[l%Bsz];
}
return s;
}
};
int main(void){
int i,t,l,r;
ll s;
const int n=1000000; // 数列長
int *a=new int[n];
// 数列を作成
for (i=0;i<n;++i){
a[i]=(double)rand() / (RAND_MAX + 1) * (INT_MAX);
}
// 構築
RangeInversions2 ri(a,n);
// クエリ1万回
for (i=0;i<10000;++i){
l=(double)rand() / (RAND_MAX + 1) * (n);
r=(double)rand() / (RAND_MAX + 1) * (n);
if (l>r) t=l,l=r,r=t;
s=ri.inversions_query(l,r);
}
delete [] a;
return 0;
}
