JavaScript: 圧縮とかPPMd

今回は例によって圧縮の話題となります。高圧縮率を誇るPPMの中でも特に高速な部類に属するPPMdについて紹介します。開発者はDmitry Shkarin氏。初版から20年以上経過していますが、いまだに現役バリバリの骨董programです。

## PPMDに搭載されている技術

## 接尾辞木
文脈情報を接尾辞木のようなもので管理する事により高速化。多分。

### Secondary Escape Estimation
これは記号数や頻度等の色々な要素から作るhash値によって、似たような連中を寄せ集めてescape確率推定を行います。

### 遺伝法
今までより長い次数の文脈が出てきた場合に、文脈を木として持っている場合ならば、親の情報を子に遺伝させることにより、始めからかなり高い予想性能を持たせます。

## 実装編
```js:sub
const N1=4,N2=4,N3=4,N4=128+3-1*N1-2*N2-3*N3>>2,
	N_INDEXES=N1+N2+N3+N4,UP_FREQ=5,MAX_FREQ=124,
	ExpEscape=new Uint8Array([51,43,18,12,11,9,8,7,6,5,4,3,3,2,2,2]), //Tabulated escapes for exponential symbol distribution
	NS2BSIndx=new Uint8Array(256),QTable=new Uint8Array(260);

function CLAMP(x,l,h){return x<l?l:h<x?h:x}
//constants initialization
function PPMD_STARTUP(){
	var i=256,k=1,m=UP_FREQ,s=1;
	for(NS2BSIndx[0]=0;i>29;)NS2BSIndx[--i]=6;
	for(NS2BSIndx[1]=NS2BSIndx[2]=2;i>3;)NS2BSIndx[--i]=4;
	for(i=0;i<m;)QTable[i]=i++;
	for(;i<260;--k||(k=++s,m++))QTable[i++]=m
}
function StartSubAllocator(as,z){
	as=++as*1048576;
	var HeapStart=0,HiUnit=HeapStart+as,LoUnit=0|HiUnit-as/8*7|0,UnitsStart=LoUnit<z?LoUnit:LoUnit=z,T=new Uint8Array(LoUnit);
	T.as=as;T.p=T.hs=HeapStart;T.hu=HiUnit;T.lu=T.us=LoUnit;return T
}
function StartModelRare(T,mo,binSum,SEE2Cont){
	var i=256,k=0,s,i2f=[],EscCoef=new Int8Array([16,-10,1,51,14,89,23,35,64,26,-42,43]),maxC={ns:255,sum:257,son:i2f,flag:0};
	T.o=mo;
	//we are in solid mode
//	if(mo<2){for(s=maxC;s.next;s=s.next)OrderFall--;return}
	// InitSubAllocator
	T.hu=(T.p=T.hs)+T.as;
	T.lu=T.us=0|T.hu-T.as/8*7;
	//alloc and init order0 context
	for(;i;)i2f[--i]={s:i,c:1};
	for(i2f=[];i<25;i2f[i++]=k+1)for(;QTable[k]==i;)k++;
	//bin SEE init
	for(k=64;k--;){
		for(s=i=0;i<6;i++)s+=EscCoef[2*i|k>>i&1];
		s=CLAMP(s,32,224)<<7;
		for(i=25;i--;)binSum[i<<6|k]=BIN_SCALE-(0|s/i2f[i])<<18
	}
	//masked SEE init
	for(;++i<23;)for(k=32;k;)SEE2Cont[i<<5|--k]=new SEE2Context(8*i+5);
	return maxC
}
var INT_BITS=7,PERIOD_BITS=7,TOT_BITS=INT_BITS+PERIOD_BITS,BIN_SCALE=1<<TOT_BITS;
//SEE-contexts for PPM-contexts with masked symbols
function SEE2Context(s){this.s=s<<this.lr}
SEE2Context.prototype={lr:PERIOD_BITS-4,c:7,
init:function(s){this.s=s<<(this.lr=PERIOD_BITS-4);this.c=7},

update:function(){if(!--this.c){
	var i=this.s>>this.lr;
	i=PERIOD_BITS-(i>40)-(i>280)-(i>1020);
	if(i<this.lr)this.s>>=1,this.lr--;
	else if(i>this.lr)this.s<<=1,this.lr++;
	this.c=5<<this.lr
}}};
function getLinks(T,hitS,p,pc,skip){
	var s=1,b=hitS.s,cf,s0,up=hitS.link,a=0,A=[];
	if(!skip)A[a++]=hitS,s=pc.next;
	if(s){
		if(p)pc=pc.next,s=0;
		do{
			if(s){
				pc=s;s=pc.son;
				//increment current symbol's freq in lower order contexts
				//more partial updates?
				if(pc.ns){
					//find b node
					for(p=0;s[p].s^b;p++);
					//increment freq if limit allows
					p=s[p];p.c<MAX_FREQ-1&&(p.c+=2,pc.sum+=2)
				}else p=s[0],p.c>15||pc.next.ns||p.c++
			}
			if(p.link!==up){pc=p.link;if(!a)return pc;break}
			A[a++]=p
		}while(s=pc.next)
	}
	skip=16*(b>63)|8*((b=T[up++])>63);
	s=pc.son;//pc is minC,the context used for encoding
	if(pc.ns){
		for(p=0;s[p].s^b;p++);
		cf=s[p].c-1;s0=pc.sum-pc.ns-cf;
		cf=2*cf<s0?1+(12*cf>s0):(cf/=s0)<5?3+cf|0:7
	}else cf=s[0].c;
	//attach the new node to all orders
	for(;A[--a].link=pc={ns:0,flag:skip,son:[{s:b,c:cf,link:up}],next:pc},a;);
	return pc
}
function ReduceOrder(T,hitS,p,pc,maxC){
	var s,b=hitS.s,i=0,q=pc,up=hitS.link=T.p;
	if(p)pc=pc.next,i--;
	for(T.o++;;T.o++){
		if(++i){
			if(!pc.next)return pc;
			pc=pc.next;s=pc.son;
			if(pc.ns){
				for(p=0;s[p].s^b;p++);
				p=s[p];p.c<MAX_FREQ-3&&(p.c+=2,pc.sum+=2)
			}else p=s[0],p.c+=p.c<11
		}
		if(p.link)break;p.link=up
	}
	if(p.link<=up)p.link=getLinks(T,p,0,pc);
	if(T.o==1&&q==maxC)
		hitS.link=p.link,T.p--;
	return p.link
}
function rescale(p,q,o,i){
	var s=q.son,of=o>0,a,n=q.ns,f0=p.c,sf=q.sum,e=sf-f0;
	for(q.flag&=20;i;)s[i]=s[--i];s[i]=p;
	for(q.sum=p.c=f0+of>>1;i<n;){
		p=s[++i];e-=a=p.c;
		q.sum+=p.c=a=a+of>>1;
		if(a&&p.s>63)q.flag|=8;
		if(a>s[i-1].c){
			for(o=i;a>(s[o--]=s[o]).c;);s[o+1]=p
		}
	}
	p=s[0];//remove the zero freq nodes
	if(!s[i].c){
		for(;e++,!s[--i].c;);s.length=i+1;
		if(!(q.ns=i)){
			q.flag&=24;
			p.c=0|Math.min(MAX_FREQ/3,(2*p.c+e-1)/e);
			return p
		}
	}
	q.sum+=e+1>>1;
	if(!of&&4&q.flag)a=2;
	else a=(sf-=e)-f0,a=CLAMP(0|(f0*q.sum-sf*p.c+a-1)/a,2,0|MAX_FREQ/2-18);
	p.c+=a;q.sum+=a;q.flag|=4;
	return p
}
function UpdateModel(T,hitS,minC,maxC,bsum){
	var fs=hitS.s,ns,cf,sf,s0,fc=hitS.c,link,flink=hitS.link,pc=minC.next,p,s;
	//partial update for the suffix context
	if(pc){
		s=pc.son;p=s[0];
		if(pc.ns){
			if(p.s^fs){
				for(p=0;s[++p].s^fs;);
				if(s[p].c>=s[p-1].c)s0=s[p],s[p--]=s[p],s[p]=s0;p=s[p]
			}
			if(p.c<MAX_FREQ-3)p.c+=cf=2+(fc<28),pc.sum+=cf
		}else p.c+=p.c<14}
	//try increasing the order
	if(!T.o&&flink){
		if(p=getLinks(T,hitS,p,minC,1))return hitS.link=p;
		return 0
	}
	T[link=T.p++]=fs;
	if(++link>=T.length)return 0;
	if(flink){
		if(flink<T.length)flink=getLinks(T,hitS,p,minC)
	}else flink=ReduceOrder(T,hitS,p,minC,maxC);
	if(!flink)return 0;
	if(!--T.o)link=flink,T.p-=maxC!==minC;
	s0=minC.sum-fc--;ns=minC.ns;
	for(pc=maxC;pc!==minC;pc=pc.next){
		s=pc.son;
		if(pc.ns)sf=pc.sum+=QTable[ns+4]>>3;//non-binary context
		else p=s[0],//escaped binary context
			p.c=p.c<=MAX_FREQ/3?2*p.c-1:MAX_FREQ-15,
			sf=pc.sum=p.c+(ns>1)+ExpEscape[QTable[bsum>>6]];//update escape
		//inheritance
		cf=fc*(5+sf);sf+=s0;
		//this is a weighted rescaling of symbol's freq into a new context(cf/sf)
		if(cf>3*sf)
			pc.sum+=cf=5+(cf>5*sf)+(cf>6*sf)+(cf>8*sf)+(cf>10*sf)+(cf>12*sf);
		else //if the new freq is too small the we increase the escape freq too
			pc.sum+=4,cf=1+(2*cf>sf)+(2*cf>3*sf);
		s[++pc.ns]={s:fs,c:cf,link:link};
		pc.flag|=8*(fs>63) //flag if last added symbol was>63
	}return flink
}
//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);
```
```js:main
/* compressor
@A: input(Array / Uint8Array)
@mo: context order(0-255). mo+1
@mm: memory usage(0-4095). mm MB
@done: call back of last process.
	done(A,a,z)
	@A: compressed array of input.
	@a: input size
	@z: compressed size
@rate: call back of progress
	rate(a,z)
	@a: current position
	@z: last position
*/
function PPMDenc(A,mo,mm,done=a=>a,rate){
	function enc(c,f,t){
		R=R/t>>>0;L+=R*c;
		for(R*=f;R<16777216;R*=256,L=L<<8>>>0)
			if((c=L>0xffffffff)||255>L>>>24)
				for(O[o++]=255&c--+B,B=L>>>24;C;C--)O[o++]=255&c;
			else++C
	}
	PPMD_STARTUP();
	var a=A.length,o=2,h=0,now=1,z=a,
		L=256,R=-1>>>0,B,C=0,// range coder
		T=StartSubAllocator(mm&=4095,z),O=[mo&=255,mm&255,mm>>8],
		D=new Uint8Array(L),X=new Uint32Array(L),binSum=new Uint32Array(1600),
		DummySEE2Cont=new SEE2Context,SEE2Cont=[],
		maxC=StartModelRare(T,++mo,binSum,SEE2Cont),
		run=mo<13?-mo:-13,run0=run,wait=wait0;

	if(typeof rate!="function")rate=a=>a,wait=a=>a();

	for(;O[++o]=B=a&255,a>>>=8;)O[2]+=16;// write size
	for(;L;)D[--L]=512/(L+8);

	function fn(){
		var d=+new Date,b,c,f,i,j,l,p,q,s,t,ns,sf,see2,masked=0,bsum;
		for(rate(a,z);a<z&&(a&8191||Date.now()-d<200);){
			q=maxC;b=A[a++];
			if(ns=q.ns){ //encode in unmasked(max order)context
				s=q.son,p=s[i=l=0],f=p.c,t=q.sum;
				if(p.s===b)p.c+=4,q.sum+=4,h=0;//2*f>t
				else s:{
					for(l=f,h=0;p=s[++i];l+=f){
						f=p.c;
						if(p.s===b){
							q.sum+=4;
							if((p.c+=4)>s[i-1].c)s[i--]=s[i],s[i]=p;break s
						}
					}
					for(masked=ns;X[s[--i].s]=now,i;);f=t-l
				}
				enc(l,f,t);
				p&&p.c>MAX_FREQ&&rescale(p,q,T.o,i)
			}
			else{ //encode in binary(deterministic)context
				p=q.son[l=0],t=binSum[i=QTable[p.c-1]<<6|NS2BSIndx[q.next.ns]+h+q.flag+(run>>26&32)]>>>9;
				h=p.s===b;binSum[i]+=((h<<23)-t)*D[i=binSum[i]&255]&~255|i<255;
				bsum=f=t>>11||1;
				if(h)p.c+=p.c<196,run++;
				else X[p.s]=now,masked=p=0,l=f,f=4096-f;
				enc(l,f,4096)
			}
			//encode in masked context
			for(;!p;enc(l,f,t)){
				for(;T.o++,q=q.next,q.ns===masked;);
				s=q.son,ns=q.ns;
				if(ns^255)
					see2=SEE2Cont[QTable[ns+3]-4<<5|(q.sum>10*(ns+1))+2*(2*ns<q.next.ns+masked)+q.flag],
					sf=1+(see2.s>>see2.lr);
				else see2=DummySEE2Cont,sf=0;
				for(i=j=f=l=0;p=s[i];i++)
					if(X[c=p.s]^now)
						X[c]=now,f+=p.c,c^b||(j=i,l=f);
				t=f+sf;
				if(l){
					p=s[j];l-=f=p.c;
					if(sf>2)see2.s-=sf;
					see2.update();q.sum+=4;
					(p.c+=4)>MAX_FREQ&&rescale(p,q,T.o,j);
					run=run0;now++
				}else masked=ns,see2.s+=l=f,f=sf
			}
			if(T.o||p.link<T.length){
				if(p=UpdateModel(T,p,q,maxC,bsum))maxC=p
			}else p=maxC=p.link;
			p||(maxC=StartModelRare(T,mo,binSum,SEE2Cont),X.fill(h=0),run=run0,now=1)
		}
		if(a<z)return wait(fn);
		for(i=5;i--;L=L<<8>>>0)
			if((c=L>0xffffffff)||255>L>>>24)
				for(O[o++]=255&c--+B,B=L>>>24;C;C--)O[o++]=255&c;
			else++C
		done(O,z,o);return O
	}
	return fn()
}
/* decompressor
@In: input(Array / Uint8Array)
@done: call back of last process.
	done(A,a,z)
	@A: decompressed array of input.
	@a: last position
	@z: decompressed size
@rate: call back of progress
	rate(a,z)
	@a: current position
	@z: last position
*/
function PPMDdec(A,done=a=>a,rate){
	function dec(c,f){for(L-=R*c,R*=f;R<16777216;R*=256)L=(L<<8|A[a++])>>>0}
	PPMD_STARTUP();

	var a=3,h=A[2]>>4,mo=A[0]+1,o=1,now=1,z=0,y=A.length,
		L=256,R=-1>>>0,I=[],O=[],
		D=new Uint8Array(L),X=new Uint32Array(L),binSum=new Uint32Array(1600),
		SEE2Cont=[],DummySEE2Cont=new SEE2Context,
		run=mo<13?-mo:-13,run0=run,wait=wait0;

	if(typeof rate!="function")rate=a=>a,wait=a=>a();

	for(;z+=A[a++]*o,h--;)o*=256;
	var T=StartSubAllocator(A[1]|A[2]<<8&4095,z),maxC=StartModelRare(T,mo,binSum,SEE2Cont);

	for(o=5;L;)D[--L]=512/(L+8);
	for(h=0;--o;)L=(L<<8|A[a++])>>>0;

	function fn(){
		var c,f,i,d=+new Date,j,l,p,q,s,t,ns,sf,see2,masked=0,bsum;
		for(rate(a,y);o<z&&(o&8191||Date.now()-d<200);){
			q=maxC;
			if(ns=q.ns){
				s=q.son,p=s[i=l=0],f=p.c,t=q.sum;c=L/(R=R/t>>>0)>>>0;
				if(c<f)p.c+=4,q.sum+=4,h=0;//2*f>t
				else s:{
					for(l=f,h=0;p=s[++i];l+=f){
						f=p.c;
						if(l+f>c){
							q.sum+=4;
							if((p.c+=4)>s[i-1].c)s[i--]=s[i],s[i]=p;break s
						}
					}
					for(masked=ns;X[s[--i].s]=now,i;);f=t-l
				}
				dec(l,f);
				p&&p.c>MAX_FREQ&&rescale(p,q,T.o,i)
			}
			else{
				p=q.son[l=0],t=binSum[i=QTable[p.c-1]<<6|NS2BSIndx[q.next.ns]+h+q.flag+(run>>26&32)]>>>9;
				bsum=f=t>>11||1;h=L/(R>>>=12)>>>0<f;
				binSum[i]+=((h<<23)-t)*D[i=binSum[i]&255]&~255|i<255;
				if(h)p.c+=p.c<196,run++;
				else X[p.s]=now,masked=p=0,l=f,f=4096-f;
				dec(l,f)
			}
			for(;!p;dec(l,f)){
				for(;T.o++,q=q.next,q.ns===masked;);
				s=q.son,ns=q.ns;
				if(ns^255)
					see2=SEE2Cont[QTable[ns+3]-4<<5|(q.sum>10*(ns+1))+2*(2*ns<q.next.ns+masked)+q.flag],
					sf=1+(see2.s>>see2.lr);
				else see2=DummySEE2Cont,sf=0;
				for(i=j=f=0;p=s[i];i++)
					if(X[c=p.s]^now)
						X[c]=now,f+=p.c,I[j++]=i;
				t=f+sf;c=L/(R=R/t>>>0)>>>0;
				if(c<f){
					for(i=l=0;(l+=s[j=I[i]].c)<=c;i++);
					p=s[j];l-=f=p.c;
					if(sf>2)see2.s-=sf;
					see2.update();q.sum+=4;
					(p.c+=4)>MAX_FREQ&&rescale(p,q,T.o,j);
					run=run0;now++
				}else masked=ns,see2.s+=l=f,f=sf
			}
			O[o++]=p.s;
			if(T.o||p.link<T.length){
				if(p=UpdateModel(T,p,q,maxC,bsum))maxC=p
			}else p=maxC=p.link;
			p||(maxC=StartModelRare(T,mo,binSum,SEE2Cont),X.fill(h=0),run=run0,now=1)
		}
		if(o<z)return wait(fn);
		done(O,a,o);return O
	}
	return fn()
}
```
また変な仕様のprogramになってやがります。解読してね。callback関数は与えなくてもよい
```js:使用例
let A=Array.from("That that is is that that is not is not is that it it is",a=>a.charCodeAt());
let e=PPMDenc(A,7,16), d=PPMDdec(e);
console.log(A.length,"->",e.length, String.fromCharCode(...d))
```
操作系CodePen
<p class="codepen" data-height="300" data-default-tab="result" data-slug-hash="QWXbdeY" data-pen-title="PPMD compressor" data-user="xezz" style="height:300px;box-sizing:border-box;display:flex;align-items:center;justify-content:center;border:2px solid;margin:1em 0;padding:1em"> <span>See the Pen <a href="https://codepen.io/xezz/pen/QWXbdeY">PPMD compressor</a> by xezz (<a href="https://codepen.io/xezz">@xezz</a>) on <a href="https://codepen.io">CodePen</a>.</span></p><script async src="https://cpwebassets.codepen.io/assets/embed/ei.js"></script>

## 参考文献
[compression.ru](http://www.compression.ru/ds/) 設計図が丸出しなのだ
[改造版](http://ctxmodel.net/files/PPMd/)

## 挫折
memory管理の実装が難し過ぎて手抜き力全開です。本家との互換性? そんなもんある訳が無いのです