静的区間転倒数

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