JavaScript: BWTとか圧縮とか

block sort(burrows wheeler transform, BWT)を利用した圧縮programを紹介します。ロードバイク並の加速力で、gzipを遥かに超える程度の圧縮力を誇るかどうか定かではないほどの性能です。

原理

処理の流れは BWT -> Move To Front -> Zero Length Encoding -> Range Coder となります。
日本語訳すると、整列により圧縮しやすくし、小さい数値に偏らせ、特に0は連長圧縮し、Range符号でとどめをさす、と言った具合です。

実装編

BWT

function compare(R,a,b,i,n,r){
	for(n-=a<b?b:a;i<n&&!(r=R[a+i]-R[b+i]);)i+=2;return r||b-a
}
function subsort(R,A,l,h,d,n){
	if(10<h-l)
		for(var i=l+2,j,m=i,t;i<h;A[j]=t)
			for(t=A[j=++i];m<j&&compare(R,t,A[j-3],d,n)<0;)A[j]=A[j-=3];
	for(i=l;i<h;A[j]=t)
		for(t=A[j=++i];l<j&&compare(R,t,A[j-1],d,n)<0;)A[j--]=A[j]
}
function median(R,A,a,b,c,d,z){
	var r=A[a]+d,s=A[b]+d;
	r=r<z?R[r]:-1,s=s<z?R[s]:-1,d+=A[c],d=d<z?R[d]:-1;
	return r<s?s<d?b:r<d?c:a:r<d?a:s<d?c:b
}
function getPivot(R,A,r,l,h,d,z){
	if(r<65)var m=l+(r>>1),p;
	else r>>=3,
		l=median(R,A,l,p=l+r,p+=r,d,z),
		m=median(R,A,p+=r,p+=r,p+=r,d,z),
		h=median(R,A,p+=r,p+=r,h,d,z);
	p=A[m=median(R,A,l,m,h,d,z)];A[m]=A[l];A[l]=p;
	return R[p+d]
}
function mqsort(R,I,l,h,d,n,L,H,D){
	for(var a,b,i,j,o,p,r,s=0,t;;){
		if(h-l<20){
			for(l<h&&subsort(R,I,l,h,d,n);l<=h;)R[I[l]]=l++;
			if(!s)break;
			l=L[--s];h=H[s];d=D[s];
			continue}
		p=getPivot(R,I,h-l+1,l,h,d,n);
		for(i=a=l,j=b=h;;I[j--]=t){
			for(;i<=j&&(r=((t=I[i]+d)<n?R[t]:-1))<=p;i++)
				r^p||(I[i]=I[a],I[a++]=t-d);
			for(;i<=j&&(r=((t=I[j]+d)<n?R[t]:-1))>=p;j--)
				r^p||(I[j]=I[b],I[b--]=t-d);
			if(i>j)break;
			t=I[i],I[i++]=I[j]}
		if((p=a-l)<(r=i-a))r=p;o=j;
		for(p=l;r--;I[o--]=t)t=I[p],I[p++]=I[o];
		if((p=h-b)<(r=b-j))r=p;o=h;
		for(p=i;r--;I[o--]=t)t=I[p],I[p++]=I[o];
		i+=l-a;j-=~h+b;
		if(j<=h)L[s]=j,H[s]=h,D[s++]=d;
		if(i<j)L[s]=i,H[s]=j-1,D[s++]=d+2;
		h=i-1}
}
//sort type B*
function sortTypeBstr(B,I,C,R,n){
	for(var A=C.slice(),p=1,r,i=n-1,j=0|Math.sqrt(n*2)+1,L=new Int32Array(j*3+3),H=L.subarray(j+1),D=H.subarray(j+1);i--;)
		if(r=B[i]-B[i+1])
			r<0&&p>0&&(I[A[B[i]<<8|B[i+1]]++]=i),p=r;
	for(;i<255;)
		for(j=++i;j<255;p-r>1&&mqsort(R,I,r,--p,2,n,L,H,D))
			p=A[r=i<<8|++j],r=C[r]
}
//set type A,B
function setTypeAB(B,I,C,n){
	var i=255,j,p,A=new Uint32Array(256),E=A.slice();
	for(A[i]=n;j=i;A[i]=E[i]){
		for(;j;)E[--j]=C[j<<8|i];
		for(j=C[i--<<8];j>E[i];)
			if((p=I[--j])>0&&B[p]>=B[--p])
				I[--E[B[p]]]=p}
	for(I[C[B[n-1]<<8]++]=n-1;j<n;)
		if((i=I[j++])>0&&B[i]<=B[--i]&&C[p=B[i]<<8|B[i+1]]<A[B[i]])
			I[C[p]++]=i
}
/* suffix sort
@T: input (Uint8Array / Array)
@A: suffix array (Int32Array / Array / null)
@n: input size (Number / null)
*/
function I2SRsort(T,A,n){n=n||T.length;
	var I=A||new Int32Array(n);
	if(n<2)return n?I[0]=0:0,I;
	if(n<3)return I[n=0|T[0]<T[1]]=1,I[1-n]=0,I;
	if(n<512&&!A)return I.map((a,b)=>b).sort((a,b,c,d)=>{for(c=a<b?n-b:n-a;c--;)if(d=T[a++]-T[b++])return d;return b-a});
	for(var i=n,c,C=new Uint32Array(65537),R=I.slice();i;)C[c|(c=T[--i])<<8]++;
	for(n=0;i<65537;n+=c)c=C[i],C[i++]=n;
	for(i=n,c=0;i;)R[--i]=C[(c|(c=T[i])<<8)+1]-1;
	sortTypeBstr(T,I,C,R,n);setTypeAB(T,I,C,n);
	return I
}
/* BWT
@T: input (Uint8Array / Array)
@U: output (Uint8Array / Array / null)
@n: input size (Number / null)
*/
function BWTe(T,U,n){n=n||T.length;
	var i=0,a=n<257?1:n<65537?2:n<16777217?3:4,p,SA=I2SRsort(T,0,n),B=U||new Uint8Array(n+a);
	for(B[0]=T[n-1];p=SA[i++];)B[i]=T[--p];
	for(p=i;i<n;)B[i]=T[SA[i++]-1];
	if(U)return p;--p;
	for(;a--;p>>=8)B[i++]=p&255;
	return B
}

main

//setImmediateもどき
(function(a){
	var F=[],hit="\0";
	a.wait0=function(f){
		F[F.length]=f;a.postMessage(hit,"*")
	};
	a.onmessage=function(e){
		if(e.source===a&&e.data===hit)
			e.stopPropagation(),F[0]&&F.shift()()
	}
})(this);
/* compressor
@A	:input(Array/Uint8Array)
@bs	:block size. 32KB-2GB
@done: call back of last process.
	done(A,a,z)
	@A: compressed/decompressed array of input.
	@a: input size
	@z: compressed size
@rate: call back of progress
	rate(a,z)
	@a: current position
	@z: last position
*/
function BMZenc(In,bs,done=a=>a,rate){
	// binary range coder
	function encBit(P,c,b){
		var p=P[c],f=x1+(x2-x1>>>12)*p;
		P[c]-=(p-(b<<12)>>5)+b;
		for(b?x2=f:x1=f+1;!((x1^x2)>>24);x2=x2<<8|255)
			Out[o++]=x2>>>24,x1<<=8
	}
	function body(){
		var b,c,step=goal,h,i=0,n,s,d=Date.now();
		for(rate(a,size);;){
			if(!--step&&Date.now(step=goal)-d>200)return wait(body);
			for(s=a<wall?In[a++]:256;MTF[i]^s;)i++;
			if(!i)run++;
			else{// encode 0 run length
				for(h=i;run;run>>=1)
					for(c=1,i=8,n=1&--run;i;c+=c+b)encBit(P254,c,b=n>>--i&1);
				// encode rank
				if(i=8,++h<255)
					for(c=1;i;c+=c+b)encBit(P254,c,b=h>>--i&1);
				else{// rank is lager than 254
					for(i=9,c=1;--i;c+=c+1)encBit(P254,c,1);
					encBit(P255,0,b=h>255);
					if(b){encBit(P255,1,b=h>256);if(b)break}
				}
				// update MTF table
				for(--h;h;)MTF[h]=MTF[--h];MTF[h]=s
			}
		}
		// write pidx
		for(i=0;width>>++i;);
		if(width>1)for(;i;encBit(P255,0,bp>>--i&1))P255[0]=2048;
		return wait(head)
	}
	function head(){
		if(a>=size){// write EOB
			for(width=a=1;a++<9;width+=width+1)encBit(P254,width,1);
			encBit(P255,0,1);encBit(P255,1,1);
			Out[o++]=x2>>>24;done(Out,size,o);
			return Out
		}
		var U=In.subarray(a,a+bs);wall+=width=U.length;
		bp=BWTe(U.slice(),U);
		return wait(body)
	}
	In instanceof Uint8Array||(In=new Uint8Array(In));
	var a=257,o=0,x1=0,x2=-1>>>0,
		run=0,width,wall=1<<15,bp;
	const size=In.length,Out=[],goal="function"!=typeof rate?0:8192,wait=wait0,
		MTF=new Uint16Array(a),
		P254=new Uint16Array(a),P255=new Uint16Array([2048,2048]);

	if(!goal)wait=a=>a(),rate=a=>a;// not async
	for(;MTF[--a]=a;)P254[a]=2048;
	bs&=-1;if(bs<wall)bs=wall;if(bs>size)bs=size;wall=0;
	return head()
}
/* decompressor
@A	:input(Array/Uint8Array)
@done: call back of last process.
	done(A,a,z)
	@A: compressed/decompressed array of input.
	@a: input size
	@z: compressed size
@rate: call back of progress
	rate(a,z)
	@a: current position
	@z: last position
*/
function BMZdec(In,done=a=>a,rate){
	function decBit(c){
		var p=P255[c],f=x1+(x2-x1>>>12)*p,b=v>>>0<=f;
		P255[c]-=(p-(b<<12)>>5)+b;
		for(b?x2=f:x1=f+1;!((x1^x2)>>24);v=v<<8|In[a++])
			x1=x1<<8>>>0,x2=x2<<8|255;
		return b
	}
	function body(){
		var b,c,i,step=goal,p,d=Date.now();
		for(rate(a,size);;){
			if(!--step&&Date.now(step=goal)-d>200)return wait(body);
			for(c=1;c<256;c+=c+b)
				for(p=P254[c],b=v>>>0<=(i=x1+(x2-x1>>>12)*p),P254[c]-=(p-(b<<12)>>5)+b,b?x2=i:x1=i+1;!((x1^x2)>>24);v=v<<8|In[a++])
					x1=x1<<8>>>0,x2=x2<<8|255;
			if((c&=255)<2)run+=++c<<n++;
			else{
				if(run)for(n=0;I[u]=Hit[U[u++]=r0]++,--run;);
				if(c-->254&&(c+=decBit(0))>254&&decBit(1))break;
				for(I[u]=Hit[U[u++]=r0=MTF[c]]++;c;)MTF[c]=MTF[--c];MTF[c]=r0
			}
		}
		if(!u){done(Out,size,o);return Out}
		for(c=b=0;u>>++c;);
		// decode index
		if(u>1)for(;c--;n|=decBit(0)<<c)P255[0]=2048;
		// unBWT
		for(c=0;u>b;b+=p)p=Hit[c],Hit[c++]=b;
		for(c=0;u;c<n&&c++)c=I[c]+Hit[Out[o+--u]=U[c]];
		o+=b;Hit.fill(r0=n=0);
		return wait(body)
	}
	var a=256,o=0,n=0,r0=0,run=0,u=0,v=2048,x1=0,x2=-1;
	const size=In.length,goal="function"!=typeof rate?0:8192,wait=wait0,
		Hit=new Uint32Array(a),I=[],Out=[],U=[],
		MTF=new Uint16Array(a),P254=new Uint16Array(a),P255=new Uint16Array([v,v]);

	for(;MTF[--a]=a;)P254[a]=v;
	for(v=0;a<4;)v=v<<8|In[a++];
	if(!goal)wait=a=>a(),rate=a=>a;// not async
	return body()
}

BMZencの第2引数には圧縮元配列の分割区間幅を指定。通常大きい方が圧縮率が高くなり、空間消費量が増え、速度が犠牲になります。
doneとrateはcallback関数ですが、指定するほどのものではありません。面倒臭いだけです。非同期関数のように振る舞い、CPU負荷が若干減るかもしれない程度のもの。

使用例

let A=Array.from("That that is is that that is not is not is that it it is",a=>a.charCodeAt());
let e=BMZenc(A,A.length), d=BMZdec(e);
console.log(A.length,"->",e.length, String.fromCharCode(...d))

file圧縮伸長実験

See the Pen BWT+MTF+ZLE compressor by xezz (@xezz) on CodePen.

余談

bzip2とほぼ同じ事をしています。大きな違いはhuffman符号ではなくRange符号を使っている点。
本作は符号化の詰めが甘ったるく、手抜き力を遺憾なく発揮しています