mandelbrot/experiments/root-cause/ref_snap_compare.mjs
2026-09-29 23:44:26 +09:00

56 lines
No EOL
9.8 KiB
JavaScript

import fs from 'node:fs';
import vm from 'node:vm';
const prodHtml=fs.readFileSync(new URL('../../index.html',import.meta.url),'utf8');
const scripts=[...prodHtml.matchAll(/<script[^>]*>([\s\S]*?)<\/script>/g)].map(m=>m[1]);
const kc={};vm.runInNewContext(scripts[0],kc);const kernels=kc.MANDEL_WEBGPU_KERNELS;
function appContext(){
const element={width:800,height:600,clientWidth:800,clientHeight:600,style:{},hidden:false,value:'',textContent:'',disabled:false,addEventListener(){},classList:{toggle(){}},setAttribute(){},removeAttribute(){},showModal(){},close(){},toBlob(){},getBoundingClientRect(){return{left:0,top:0,width:800,height:600}}};
const context={MANDEL_WEBGPU_KERNELS:kernels,document:{querySelector:()=>element,getElementById:()=>element,documentElement:{classList:{toggle(){}}},body:{classList:{toggle(){}},appendChild(){}}},navigator:{hardwareConcurrency:4,deviceMemory:8},location:{search:'',hash:''},history:{replaceState(){}},performance,URLSearchParams,URL:{createObjectURL:()=>`blob:test`,revokeObjectURL(){}},TextEncoder,Blob:class{},console,addEventListener(){},requestAnimationFrame:()=>1,cancelAnimationFrame(){},setTimeout,clearTimeout,setInterval,clearInterval,innerWidth:800,innerHeight:600,devicePixelRatio:1,localStorage:{getItem(){return null},setItem(){},removeItem(){}}};
const cut=scripts[1].indexOf('// ── boot / teardown');
const code=scripts[1].slice(0,cut)+`\nglobalThis.internal={state,fastReferenceWorkerSource,precisionFallbackWorkerSource,fastReferenceBits,pixelPrecisionBits,referencePixelForSource,spanMantExp,adaptiveProbeChoice};})();`;
vm.runInNewContext(code,context);return context.internal;
}
const app=appContext();
function runWorker(source,data){let result;const c={self:{},postMessage:d=>{result=d},performance,console,BigInt,Math,ArrayBuffer,DataView,Float32Array,Uint32Array,Uint8Array};vm.runInNewContext(source,c);c.self.onmessage({data});if(!result)throw new Error('worker returned no result');if(result.type==='error')throw new Error(result.error);return result;}
function buildFast(snap,iter,w,h){const targetBits=app.fastReferenceBits(snap,w),d=runWorker(app.fastReferenceWorkerSource(),{type:'build',id:1,key:'analysis',sourceBits:snap.bits,targetBits,re:snap.re.toString(),im:snap.im.toString(),span:snap.span.toString(),width:w,height:h,iter});d.source={bits:snap.bits,re:BigInt(d.referenceRe),im:BigInt(d.referenceIm),viewRe:snap.re,viewIm:snap.im,span:snap.span,iter};return d;}
function oracle(snap,iter,w,h,indices){const targetBits=app.pixelPrecisionBits(snap,w,40);return runWorker(app.precisionFallbackWorkerSource(),{type:'solve',id:2,bits:snap.bits,targetBits,re:snap.re.toString(),im:snap.im.toString(),span:snap.span.toString(),width:w,height:h,iter,terminalClass:-1,indices:Uint32Array.from(indices).buffer});}
const f=Math.fround,U=5.960464477539063e-8;
const add=(a,b)=>f(f(a)+f(b)), mul=(a,b)=>f(f(a)*f(b));
function cmul(a,b){return [add(mul(a[0],b[0]),-mul(a[1],b[1])),add(mul(a[0],b[1]),mul(a[1],b[0]))]}
const maxabs=v=>Math.max(Math.abs(v[0]),Math.abs(v[1]));
function ld(v,e){return f(v*Math.pow(2,e))}
function refsFrom(buf,n){const dv=new DataView(buf),a=[];for(let i=0;i<=n;i++){const o=i*16;a.push([[dv.getFloat32(o,true),dv.getFloat32(o+4,true)],[dv.getFloat32(o+8,true),dv.getFloat32(o+12,true)]])}return a}
function dsQuick(a,b){const q=f(a+b),e=f(b-f(q-a));return [q,e]}
function dsSum(a,b){const q=f(a+b),bb=f(q-a),e=f(f(a-f(q-bb))+f(b-bb));return [q,e]}
function dsProd(a,b){const q=f(a*b),ca=f(4097*a),ah=f(ca-f(ca-a)),al=f(a-ah),cb=f(4097*b),bh=f(cb-f(cb-b)),bl=f(b-bh);let e=f(f(ah*bh)-q);e=f(e+f(ah*bl));e=f(e+f(al*bh));e=f(e+f(al*bl));return [q,e]}
function dsAdd(a,b){const q=dsSum(a[0],b[0]);return dsQuick(q[0],f(q[1]+f(a[1]+b[1])))}
const dsNeg=a=>[-a[0],-a[1]],dsSub=(a,b)=>dsAdd(a,dsNeg(b));
function dsMul(a,b){const q=dsProd(a[0],b[0]);let e=f(q[1]+f(a[0]*b[1]));e=f(e+f(a[1]*b[0]));e=f(e+f(a[1]*b[1]));return dsQuick(q[0],e)}
function dsScale(a,b){const q=dsProd(a[0],b);return dsQuick(q[0],f(q[1]+f(a[1]*b)))}
const dsPow2=(a,e)=>[ld(a[0],e),ld(a[1],e)],dsVal=a=>f(a[0]+a[1]);
const cdsAdd=(a,b)=>[dsAdd(a[0],b[0]),dsAdd(a[1],b[1])];
const cdsMul=(a,b)=>[dsSub(dsMul(a[0],b[0]),dsMul(a[1],b[1])),dsAdd(dsMul(a[0],b[1]),dsMul(a[1],b[0]))];
const cdsScale=(a,b)=>[dsScale(a[0],b),dsScale(a[1],b)],cdsPow2=(a,e)=>[dsPow2(a[0],e),dsPow2(a[1],e)];
const cdsMag2=a=>dsAdd(dsMul(a[0],a[0]),dsMul(a[1],a[1])),cdsMax=a=>Math.max(Math.abs(dsVal(a[0])),Math.abs(dsVal(a[1])));
function simulateDS(snap,w,h,iter,ctx,index){const px=index%w,py=Math.floor(index/w),rpix=app.referencePixelForSource(ctx.source,snap,w,h),se=app.spanMantExp(snap),mh=f(se.mant),ml=f(se.mant-mh),iw=1/w,iwh=f(iw),iwl=f(iw-iwh),gx=f(px+.5),gy=f(py+.5),ox=f(gx-f(rpix.x)),oy=f(f(rpix.y)-gy),inv=[iwh,iwl],sm=[mh,ml],dx=dsScale(inv,ox),dy=dsScale(inv,oy),d0=[dsMul(sm,dx),dsMul(sm,dy)];let d=d0,ww=[[0,0],[0,0]],scaleExp=se.exp,n=0,m=0;const refs=refsFrom(ctx.refs,ctx.refLen);while(true){const rp=refs[m],r=[[rp[0][0],rp[1][0]],[rp[0][1],rp[1][1]]],delta=cdsPow2(ww,scaleExp),z=cdsAdd(r,delta),mag=cdsMag2(z);if(mag[0]>4||(mag[0]===4&&mag[1]>0))return{kind:'escape',n};if(n>=iter)return{kind:'op',n};if(m>=ctx.refLen)return{kind:'refend',n};const linear=cdsScale(cdsMul(r,ww),2),sq=cdsPow2(cdsMul(ww,ww),scaleExp);ww=cdsAdd(cdsAdd(linear,sq),d);m++;n++;const mm=Math.max(cdsMax(ww),cdsMax(d));if(!Number.isFinite(mm)||mm>=1e30)return{kind:'range',n};if(mm>65536){ww=cdsScale(ww,1/65536);d=cdsScale(d,1/65536);scaleExp+=16}else if(mm>0&&mm<1/65536&&scaleExp>se.exp){ww=cdsScale(ww,65536);d=cdsScale(d,65536);scaleExp-=16}if(scaleExp>126)return{kind:'range',n};}}
function simulate(snap,w,h,iter,ctx,index,guarded){const px=index%w,py=Math.floor(index/w),rpix=app.referencePixelForSource(ctx.source,snap,w,h),se=app.spanMantExp(snap),mh=f(se.mant),ml=f(se.mant-mh),iw=1/w,iwh=f(iw),iwl=f(iw-iwh),gx=f(px+.5),gy=f(py+.5),ox=f(gx-f(rpix.x)),oy=f(f(rpix.y)-gy);let d;
if(guarded){const dxh=mul(ox,iwh),dxl=mul(ox,iwl),dyh=mul(oy,iwh),dyl=mul(oy,iwl);d=[add(mul(mh,dxh),add(mul(mh,dxl),mul(ml,dxh))),add(mul(mh,dyh),add(mul(mh,dyl),mul(ml,dyh)))];}
else {const dx=f(ox/f(w)),dy=f(oy/f(w));d=[mul(mh,dx),mul(mh,dy)];}
let ww=[0,0],scaleExp=se.exp,n=0,m=0,err=4*U*(maxabs(d)+1e-30);const refs=refsFrom(ctx.refs,ctx.refLen);
while(true){const rp=refs[Math.min(m,ctx.refLen)],delta=[ld(ww[0],scaleExp),ld(ww[1],scaleExp)],z=[add(rp[0][0],add(rp[1][0],delta[0])),add(rp[0][1],add(rp[1][1],delta[1]))],mag=add(mul(z[0],z[0]),mul(z[1],z[1]));
if(mag>4){if(!guarded)return{kind:'escape',n};const errAbs=Math.abs(ld(err,scaleExp))+U*(maxabs(z)+maxabs(delta)+1e-30);if(Math.hypot(z[0],z[1])-errAbs>2)return{kind:'escape',n};return{kind:'uncertain',n};}
if(n>=iter)return{kind:'op',n};if(m>=ctx.refLen)return{kind:'refend',n};
const refAbs=maxabs(rp[0])+maxabs(rp[1]),wAbs=maxabs(ww),dAbs=maxabs(d),deltaAbs=maxabs(delta),ww2=cmul(ww,ww),sq=[ld(ww2[0],scaleExp),ld(ww2[1],scaleExp)],sqAbs=maxabs(sq),gain=2*refAbs+2*deltaAbs,roundErr=U*(2*refAbs*wAbs+sqAbs+dAbs+1e-30);if(guarded)err=f(gain*err+roundErr);
const l1=cmul(rp[0],ww),l2=cmul(rp[1],ww),linear=[f(2*f(l1[0]+l2[0])),f(2*f(l1[1]+l2[1]))];ww=[add(add(linear[0],sq[0]),d[0]),add(add(linear[1],sq[1]),d[1])];m++;n++;const mm=Math.max(maxabs(ww),maxabs(d));if(!Number.isFinite(mm)||mm>=1e30||guarded&&(!Number.isFinite(err)||err>1e35))return{kind:'range',n};if(mm>65536){ww=[mul(ww[0],1/65536),mul(ww[1],1/65536)];d=[mul(d[0],1/65536),mul(d[1],1/65536)];if(guarded)err=f(err/65536);scaleExp+=16}else if(mm>0&&mm<1/65536&&scaleExp>se.exp){ww=[mul(ww[0],65536),mul(ww[1],65536)];d=[mul(d[0],65536),mul(d[1],65536)];if(guarded)err=f(err*65536);scaleExp-=16}if(scaleExp>126)return{kind:'range',n};
}
}
function grid(w,h,nx=8,ny=5){const out=[];for(let j=0;j<ny;j++)for(let i=0;i<nx;i++){const x=Math.min(w-1,Math.floor((i+.5)*w/nx)),y=Math.min(h-1,Math.floor((j+.5)*h/ny));out.push(y*w+x)}return out}
const cases=[
{name:'zoom-black',bits:273,re:'-20475306797689005889663451277157868621960315439319427815381950343245479294092514544',im:'966333796787038295532838729476816377447677486901484166221710883227030438393551731',sp:'14873556985132116323960339070653445992212160392759733161044631107644716'},
{name:'specific-black-magnification',bits:273,re:'-1186119846127651495122135785614819020200240938545203313122293422684276395840442090',im:'13357094062527643855275324732358464420434330930839287896478877080241737274649035786',sp:'137791365414599997721476797879067370932659138165656231845353593015922659734848'},
{name:'undrawn-area',bits:307,re:'46025019095185088783974324228976729149454131232585283954369753562823949096514860778965742615',im:'-151398041000872486580541963077907528105710741688397340075259729545288600585999778752201285140',sp:'1361049793375351627487861104794890130303800780964576514408724611446175'}
];
function gridN(w,h,nx=6,ny=4){const out=[];for(let j=0;j<ny;j++)for(let i=0;i<nx;i++){const x=Math.min(w-1,Math.floor((i+.5)*w/nx)),y=Math.min(h-1,Math.floor((j+.5)*h/ny));out.push(y*w+x)}return out}
const its={'zoom-black':1024,'specific-black-magnification':1536,'undrawn-area':4702};const rows=[];function snap(refLen,iter){const choices=[350,512,768,1024,1536,2048].filter((x,i,a)=>x<=refLen&&a.indexOf(x)===i);return choices.length?Math.max(...choices):Math.max(1,Math.min(iter,Math.max(64,Math.floor(refLen/64)*64)))}for(const c of cases){const snapv={bits:c.bits,re:BigInt(c.re),im:BigInt(c.im),span:BigInt(c.sp)},w=1365,h=768,requested=its[c.name],ctx=buildFast(snapv,requested,w,h),cur=snap(ctx.refLen,requested),indices=gridN(w,h,8,6);for(const iter of [cur,ctx.refLen]){const o=oracle(snapv,iter,w,h,indices),meta=new Uint32Array(o.meta);let exactEsc=0,fast={};for(let k=0;k<indices.length;k++){if(((meta[k]>>>28)&3)===1)exactEsc++;const q=simulate(snapv,w,h,iter,ctx,indices[k],false);fast[q.kind]=(fast[q.kind]||0)+1;}rows.push({case:c.name,requested,refLen:ctx.refLen,iter,exactEsc,samples:indices.length,fast})}}console.log(JSON.stringify(rows,null,2));