// Exact pedagogical models. Qubit q0 is the most significant bit.
export const C=(r=0,i=0)=>[r,i];
export const add=(a,b)=>[a[0]+b[0],a[1]+b[1]];
export const sub=(a,b)=>[a[0]-b[0],a[1]-b[1]];
export const mul=(a,b)=>[a[0]*b[0]-a[1]*b[1],a[0]*b[1]+a[1]*b[0]];
export const scale=(a,s)=>[a[0]*s,a[1]*s];
export const conj=a=>[a[0],-a[1]];
export const abs2=a=>a[0]*a[0]+a[1]*a[1];
export const phase=t=>[Math.cos(t),Math.sin(t)];
export const norm=v=>Math.sqrt(v.reduce((s,z)=>s+abs2(z),0));
export function normalize(v){const n=norm(v);if(!Number.isFinite(n)||n<1e-14)throw Error("Enter a nonzero, finite state.");return v.map(z=>scale(z,1/n));}
export const inner=(a,b)=>a.reduce((s,z,j)=>add(s,mul(conj(z),b[j])),C());
export const ket=(n,x=0)=>Array.from({length:2**n},(_,j)=>C(j===x?1:0));
export const tensor=(a,b)=>a.flatMap(x=>b.map(y=>mul(x,y)));
export const probabilities=v=>v.map(abs2);
export const bits=(x,n)=>x.toString(2).padStart(n,"0");
export const popcount=x=>{let s=0;while(x){s+=x&1;x>>>=1;}return s;};
export const H=[[C(1/Math.sqrt(2)),C(1/Math.sqrt(2))],[C(1/Math.sqrt(2)),C(-1/Math.sqrt(2))]];
export const X=[[C(),C(1)],[C(1),C()]];
export const Y=[[C(),C(0,-1)],[C(0,1),C()]];
export const Z=[[C(1),C()],[C(),C(-1)]];
export const I=[[C(1),C()],[C(),C(1)]];
export function gate(name,t=0){
 const c=Math.cos(t/2),s=Math.sin(t/2);
 return ({H,X,Y,Z,I,S:[[C(1),C()],[C(),C(0,1)]],Sdg:[[C(1),C()],[C(),C(0,-1)]],
 T:[[C(1),C()],[C(),phase(Math.PI/4)]],
 Rx:[[C(c),C(0,-s)],[C(0,-s),C(c)]],
 Ry:[[C(c),C(-s)],[C(s),C(c)]],
 Rz:[[phase(-t/2),C()],[C(),phase(t/2)]],
 P:[[C(1),C()],[C(),phase(t)]]})[name]??(()=>{throw Error("Unknown gate "+name);})();
}
export function apply1(v,n,q,g,controls=[]){
 if(!Number.isInteger(q)||q<0||q>=n||controls.includes(q)||new Set(controls).size!==controls.length||controls.some(c=>c<0||c>=n))throw Error("Invalid gate wires.");
 const out=v.map(z=>z.slice()),mask=2**(n-1-q),cm=controls.reduce((m,c)=>m|2**(n-1-c),0);
 for(let a=0;a<v.length;a++)if(!(a&mask)&&(a&cm)===cm){let b=a|mask;
 out[a]=add(mul(g[0][0],v[a]),mul(g[0][1],v[b]));
 out[b]=add(mul(g[1][0],v[a]),mul(g[1][1],v[b]));}
 return out;
}
export function swap(v,n,a,b){
 const ma=2**(n-1-a),mb=2**(n-1-b);
 return v.map((_,j)=>v[Boolean(j&ma)!==Boolean(j&mb)?j^ma^mb:j].slice());
}
export function rng(seed=1){let a=Number(seed)>>>0;return ()=>{a+=0x6D2B79F5;let t=a;t=Math.imul(t^t>>>15,t|1);t^=t+Math.imul(t^t>>>7,t|61);return ((t^t>>>14)>>>0)/4294967296;};}
export function draw(p,random=Math.random){let u=random()*p.reduce((a,b)=>a+b,0);for(let j=0;j<p.length;j++){u-=p[j];if(u<0)return j;}return p.length-1;}
export function sample(p,shots=1024,seed=1){const count=p.map(()=>0),r=typeof seed==="function"?seed:rng(seed);for(let k=0;k<shots;k++)count[draw(p,r)]++;return count;}
export function measure(v,n,q,r=Math.random){
 const mask=2**(n-1-q),p1=v.reduce((s,z,j)=>s+((j&mask)?abs2(z):0),0),value=r()<p1?1:0,p=value?p1:1-p1;
 return {value,probability:p,state:v.map((z,j)=>Boolean(j&mask)===Boolean(value)?scale(z,1/Math.sqrt(p)):C())};
}
export function circuit(n,ops,seed=1,initial=ket(n)){
 let v=initial,record={},history=[{state:v,record:{},label:"Preparation"}],r=rng(seed);
 for(const op of ops){
  let label=op.name;
  if(op.condition&&record[op.condition.q]===undefined)throw Error("Measure q"+op.condition.q+" before using its classical result.");
  if(op.condition&&record[op.condition.q]!==op.condition.value){history.push({state:v,record:{...record},label:"Condition false: "+label});continue;}
  if(op.name==="M"){const m=measure(v,n,op.q,r);v=m.state;record[op.q]=m.value;label="M(q"+op.q+") = "+m.value;}
  else if(op.name==="SWAP")v=swap(v,n,op.q,op.control);
  else v=apply1(v,n,op.q,gate(op.name==="CX"?"X":op.name==="CZ"?"Z":op.name==="CCX"?"X":op.name==="CP"?"P":op.name,op.angle||0),op.name==="CCX"?[op.control,op.control2]:["CX","CZ","CP"].includes(op.name)?[op.control]:[]);
  history.push({state:v,record:{...record},label});
 }
 return {state:v,record,history};
}
export function circuitShots(n,ops,shots,seed=1){
 if(!ops.some(o=>o.name==="M"))return sample(probabilities(circuit(n,ops,seed).state),shots,seed);
 const counts=Array(2**n).fill(0),r=rng(seed);
 for(let k=0;k<shots;k++){const v=circuit(n,ops,Math.floor(r()*2**32)).state;counts[draw(probabilities(v),r)]++;}
 return counts;
}
export const basisVectors={
 Z:[[C(1),C()],[C(),C(1)]],
 X:[[C(Math.SQRT1_2),C(Math.SQRT1_2)],[C(Math.SQRT1_2),C(-Math.SQRT1_2)]],
 Y:[[C(Math.SQRT1_2),C(0,Math.SQRT1_2)],[C(Math.SQRT1_2),C(0,-Math.SQRT1_2)]]
};
export function projection(raw,basis="Z"){const v=normalize(raw),a=basisVectors[basis].map(b=>inner(b,v));return {state:v,norm:norm(raw),amplitudes:a,probabilities:a.map(abs2)};}
export function bloch(v){const z=mul(conj(v[0]),v[1]);return [2*z[0],2*z[1],abs2(v[0])-abs2(v[1])];}
export const blochState=(theta,phi)=>[C(Math.cos(theta/2)),scale(phase(phi),Math.sin(theta/2))];
export const outer=v=>v.map(a=>v.map(b=>mul(a,conj(b))));
export function matmul(a,b){return a.map(row=>b[0].map((_,j)=>row.reduce((s,x,k)=>add(s,mul(x,b[k][j])),C())));}
export const dagger=a=>a[0].map((_,j)=>a.map(row=>conj(row[j])));
export const trace=a=>a.reduce((s,row,j)=>add(s,row[j]),C());
export function reduced(v,n,keep){
 const out=Array.from({length:2**keep.length},()=>Array.from({length:2**keep.length},()=>C()));
 const masks=keep.map(q=>2**(n-1-q)),environment=Array.from({length:n},(_,j)=>j).filter(q=>!keep.includes(q));
 function index(x,qs){return qs.reduce((a,q)=>a*2+((x>>>(n-1-q))&1),0);}
 for(let a=0;a<v.length;a++)for(let b=0;b<v.length;b++)if(index(a,environment)===index(b,environment)){
 const j=index(a,keep),k=index(b,keep);out[j][k]=add(out[j][k],mul(v[a],conj(v[b])));}
 return out;
}
export const purity=rho=>trace(matmul(rho,rho))[0];
export const rhoBloch=rho=>[2*rho[0][1][0],-2*rho[0][1][1],rho[0][0][0]-rho[1][1][0]];
export function entropyQubit(rho){let r=Math.min(1,Math.hypot(...rhoBloch(rho)));return binaryEntropy((1+r)/2);}
export function binaryEntropy(p){if(p<=0||p>=1)return 0;return -p*Math.log2(p)-(1-p)*Math.log2(1-p);}
export function channel(rho,name,p){
 if(!(p>=0&&p<=1))throw Error("Channel strength must lie between zero and one.");
 let ks;
 if(name==="amplitude")ks=[[[C(1),C()],[C(),C(Math.sqrt(1-p))]],[[C(),C(Math.sqrt(p))],[C(),C()]]];
 else if(name==="phase")ks=[I.map(row=>row.map(z=>scale(z,Math.sqrt(1-p)))),Z.map(row=>row.map(z=>scale(z,Math.sqrt(p))))];
 else if(name==="bit")ks=[I.map(row=>row.map(z=>scale(z,Math.sqrt(1-p)))),X.map(row=>row.map(z=>scale(z,Math.sqrt(p))))];
 else if(name==="depolarizing")ks=[I.map(row=>row.map(z=>scale(z,Math.sqrt(1-3*p/4)))),...[X,Y,Z].map(g=>g.map(row=>row.map(z=>scale(z,Math.sqrt(p/4)))))];
 else throw Error("Unknown channel.");
 const out=rho.map(row=>row.map(()=>C()));
 for(const k of ks){const t=matmul(matmul(k,rho),dagger(k));for(let a=0;a<2;a++)for(let b=0;b<2;b++)out[a][b]=add(out[a][b],t[a][b]);}
 return out;
}
export function chsh(angles,visibility=1,shots=1024,seed=1){
 const [a,ap,b,bp]=angles,es=[a-b,a-bp,ap-b,ap-bp].map(t=>visibility*Math.cos(t)),r=rng(seed);
 const counts=es.map(e=>sample([(1+e)/4,(1-e)/4,(1-e)/4,(1+e)/4],shots,r));
 const estimates=counts.map(c=>(c[0]+c[3]-c[1]-c[2])/shots);
 return {correlations:es,counts,estimates,S:es[0]+es[1]+es[2]-es[3],observed:estimates[0]+estimates[1]+estimates[2]-estimates[3],standardError:Math.sqrt(es.reduce((s,e)=>s+(1-e*e)/shots,0))};
}
export function teleport(input){
 let v=tensor(normalize(input),ket(2));
 v=apply1(v,3,1,H);v=apply1(v,3,2,X,[1]);const shared=v;
 v=apply1(v,3,1,X,[0]);v=apply1(v,3,0,H);
 const branches=[];
 for(let a=0;a<2;a++)for(let b=0;b<2;b++){
 const raw=[v[4*a+2*b],v[4*a+2*b+1]],p=raw.reduce((s,z)=>s+abs2(z),0),before=normalize(raw);
 let after=before;if(b)after=apply1(after,1,0,X);if(a)after=apply1(after,1,0,Z);
 branches.push({a,b,p,before,after,fidelity:abs2(inner(normalize(input),after))});
 }
 return {shared,bellMeasured:v,branches};
}
export function dense(a,b){
 let v=apply1(ket(2),2,0,H);v=apply1(v,2,1,X,[0]);const shared=v;
 if(b)v=apply1(v,2,0,X);if(a)v=apply1(v,2,0,Z);const encoded=v;
 v=apply1(v,2,1,X,[0]);v=apply1(v,2,0,H);return {shared,encoded,decoded:v,probabilities:probabilities(v)};
}
export function bb84(shots,eveRate=0,noise=0,seed=1){
 const r=rng(seed),rows=[];let sifted=0,errors=0;
 for(let k=0;k<shots;k++){
 let bit=r()<.5?0:1,ab=r()<.5?0:1,bb=r()<.5?0:1,v=basisVectors[ab?"X":"Z"][bit],eb=null,ev=null;
 if(r()<eveRate){eb=r()<.5?0:1;ev=draw(basisVectors[eb?"X":"Z"].map(b=>abs2(inner(b,v))),r);v=basisVectors[eb?"X":"Z"][ev];}
 let out=draw(basisVectors[bb?"X":"Z"].map(b=>abs2(inner(b,v))),r);if(r()<noise)out^=1;
 const keep=ab===bb;if(keep){sifted++;if(out!==bit)errors++;}
 rows.push({bit,ab,bb,eb,ev,out,keep,error:keep&&out!==bit});
 }
 const expectedQber=noise+(1-2*noise)*eveRate/4;
 return {rows,sifted,errors,qber:sifted?errors/sifted:null,expectedQber,asymptoticFraction:Math.max(0,1-2*binaryEntropy(expectedQber))};
}
export function phaseOracle(v,f){return v.map((z,x)=>scale(z,f(x)?-1:1));}
export function hadamards(v,n){for(let q=0;q<n;q++)v=apply1(v,n,q,H);return v;}
export function oracleAlgorithm(n,type,secret=1){
 let f;
 if(type==="Deutsch") {n=1;const table=secret&3;f=x=>(table>>x)&1;}
 else if(type==="DJ-constant")f=()=>secret&1;
 else f=x=>popcount(x&secret)&1;
 const start=ket(n),uniform=hadamards(start,n),marked=phaseOracle(uniform,f),out=hadamards(marked,n);
 return {n,table:Array.from({length:2**n},(_,x)=>f(x)),history:[start,uniform,marked,out],state:out,probabilities:probabilities(out)};
}
export function simon(n,s,shots=12,seed=1){
 if(s<=0||s>=2**n)throw Error("Simon's two-to-one promise requires a nonzero n-bit mask.");
 const N=2**n,f=x=>Math.min(x,x^s),pairs=Array.from({length:N},(_,x)=>({x,fx:f(x)}));
 // Trace out the output register of sum_x |x>|f(x)>/sqrt(N).
 // H^n on each input pair yields the exact marginal distribution.
 const p=Array(N).fill(0);
 for(let y=0;y<N;y++)for(let x=0;x<N;x++)for(const xp of [x,x^s]){
 const sign=((popcount(x&y)+popcount(xp&y))&1)?-1:1;p[y]+=sign/(N*N);}
 const r=rng(seed),observations=Array.from({length:shots},()=>draw(p,r)),candidates=Array.from({length:N},(_,x)=>x).filter(x=>observations.every(y=>(popcount(x&y)&1)===0));
 let rows=observations.slice(),rank=0;
 for(let bit=n-1;bit>=0;bit--){const at=rows.findIndex((x,j)=>j>=rank&&(x&(1<<bit)));if(at<0)continue;[rows[at],rows[rank]]=[rows[rank],rows[at]];for(let j=0;j<rows.length;j++)if(j!==rank&&(rows[j]&(1<<bit)))rows[j]^=rows[rank];rank++;}
 return {pairs,probabilities:p,observations,candidates,rank,rref:rows.slice(0,rank)};
}
export function grover(n,marked,iterations){
 const N=2**n,targets=new Set(marked);if(!targets.size||[...targets].some(x=>x<0||x>=N))throw Error("Select a valid marked item.");
 let v=Array.from({length:N},()=>C(1/Math.sqrt(N))),history=[v];
 for(let k=0;k<iterations;k++){v=v.map((z,x)=>scale(z,targets.has(x)?-1:1));const mean=scale(v.reduce(add,C()),1/N);v=v.map(z=>sub(scale(mean,2),z));history.push(v);}
 return {state:v,history,probabilities:probabilities(v),success:[...targets].reduce((s,x)=>s+abs2(v[x]),0),theory:Math.sin((2*iterations+1)*Math.asin(Math.sqrt(targets.size/N)))**2};
}
export function qft(v,inverse=false){
 const N=v.length,sgn=inverse?-1:1;
 return Array.from({length:N},(_,y)=>scale(v.reduce((s,z,x)=>add(s,mul(z,phase(sgn*2*Math.PI*x*y/N))),C()),1/Math.sqrt(N)));
}
export function qpe(m,phi){
 const N=2**m,before=Array.from({length:N},(_,x)=>scale(phase(2*Math.PI*phi*x),1/Math.sqrt(N))),state=qft(before,true);
 return {before,state,probabilities:probabilities(state)};
}
export function gcd(a,b){while(b){[a,b]=[b,a%b];}return a;}
export function modpow(a,b,n){let v=1;while(b>0){if(b&1)v=(v*a)%n;a=a*a%n;b=Math.floor(b/2);}return v;}
export function order(a,N){if(gcd(a,N)!==1)return null;let v=1;for(let r=1;r<=N;r++){v=v*a%N;if(v===1)return r;}return null;}
export function continuedFractions(a,b){
 let p0=0,p1=1,q0=1,q1=0,out=[];
 while(b){let c=Math.floor(a/b),p=c*p1+p0,q=c*q1+q0;out.push({coefficient:c,p,q});[p0,p1,q0,q1]=[p1,p,q1,q];[a,b]=[b,a%b];}
 return out;
}
export function shor(N,a,m=8){
 if(gcd(a,N)!==1)return {gcd:gcd(a,N)};
 const Q=2**m,groups=new Map();
 for(let x=0;x<Q;x++){const y=modpow(a,x,N);if(!groups.has(y))groups.set(y,[]);groups.get(y).push(x);}
 const p=Array(Q).fill(0);
 // Exact marginal after the inverse QFT, summing incoherently over the work register.
 for(const xs of groups.values())for(let y=0;y<Q;y++){
 let z=C();for(const x of xs)z=add(z,phase(-2*Math.PI*x*y/Q));p[y]+=abs2(z)/(Q*Q);}
 return {N,a,m,Q,groups:[...groups.entries()],probabilities:p,referenceOrder:order(a,N)};
}
export function factorFromOutcome(N,a,m,y){
 const cf=continuedFractions(y,2**m),tried=[],factors=[];
 for(const c of cf){if(c.q<2||c.q>=N)continue;
 for(let mult=1;mult<=Math.min(8,Math.floor(N/c.q));mult++){
 const r=c.q*mult;if(tried.includes(r))continue;tried.push(r);
 if(modpow(a,r,N)!==1||r%2)continue;const t=modpow(a,r/2,N);
 for(const f of [gcd(t-1,N),gcd(t+1,N)])if(f>1&&f<N&&!factors.includes(f))factors.push(f);
 }}
 return {cf,tried,factors};
}
export function repetition(p,d=3){let sum=0,bin=1;for(let k=0;k<=d;k++){if(k>=(d+1)/2)sum+=bin*p**k*(1-p)**(d-k);bin=bin*(d-k)/(k+1);}return sum;}
export function repetitionSyndrome(error,d=3){
 const e=Array.from({length:d},(_,j)=>(error>>(d-1-j))&1),syndrome=e.slice(1).map((x,j)=>x^e[j]);
 // Reconstruct the two compatible patterns and choose the smaller Hamming weight.
 const c=[0];for(const s of syndrome)c.push(c.at(-1)^s);
 if(c.reduce((a,b)=>a+b,0)>d/2)for(let j=0;j<d;j++)c[j]^=1;
 const residual=e.map((x,j)=>x^c[j]);
 return {error:e,syndrome,correction:c,residual,logicalFailure:residual.every(Boolean)};
}
export function pauliCommutes(a,b){let flips=0;for(let j=0;j<a.length;j++)if(a[j]!=="I"&&b[j]!=="I"&&a[j]!==b[j])flips++;return flips%2===0;}
// Commuting checks for the distance-three rotated planar code.
export function rotatedSurface3(){
 const xs=["XXIIIIIII","IXXIXXIII","IIIXXIXXI","IIIIIIIXX"];
 const zs=["IIIZIIZII","ZZIZZIIII","IIIIZZIZZ","IIZIIZIII"];
 return [...xs.map((p,j)=>({name:"X"+j,p,type:"X"})),...zs.map((p,j)=>({name:"Z"+j,p,type:"Z"}))];
}
export function syndromeFor(error,checks=rotatedSurface3()){return checks.map(c=>pauliCommutes(error,c.p)?0:1);}
export function resource({p=.001,threshold=.01,A=.1,d=9,locations=1e6,logical=10,cycle=1e-6,depth=1e5,factories=0}){
 const logicalError=A*(p/threshold)**((d+1)/2),physicalPerPatch=2*d*d-1,patches=logical+factories;
 return {logicalError,unionBound:Math.min(1,locations*logicalError),physicalPerPatch,physical:physicalPerPatch*patches,time:depth*d*cycle,below:p<threshold};
}
export function vqe(theta,phi=0){const v=blochState(theta,phi),r=bloch(v);return {state:v,energy:r[2]+.5*r[0],ground:-Math.sqrt(1.25),gradient:-Math.sin(theta)+.5*Math.cos(theta)*Math.cos(phi)};}
export function qaoa(gamma,beta){
 const n=3,N=8,edges=[[0,1],[1,2],[0,2]],cost=x=>edges.reduce((s,[a,b])=>s+(((x>>(n-1-a))&1)!==((x>>(n-1-b))&1)?1:0),0);
 let v=Array.from({length:N},(_,x)=>scale(phase(-gamma*cost(x)),1/Math.sqrt(N)));
 for(let q=0;q<n;q++)v=apply1(v,n,q,gate("Rx",2*beta));
 const p=probabilities(v);return {state:v,probabilities:p,costs:Array.from({length:N},(_,x)=>cost(x)),expectation:p.reduce((s,z,x)=>s+z*cost(x),0),optimum:2};
}
export function hamiltonian(t,steps){
 // H=X+Z. Exact exp(-itH) because H^2=2I.
 const a=Math.sqrt(2)*t,c=Math.cos(a),s=Math.sin(a)/Math.sqrt(2),exact=[C(c,-s),C(0,-s)];
 let approx=ket(1),dt=t/steps;
 for(let j=0;j<steps;j++){approx=apply1(approx,1,0,gate("Rz",2*dt));approx=apply1(approx,1,0,gate("Rx",2*dt));}
 return {exact,approx,fidelity:abs2(inner(exact,approx)),error:norm(exact.map((z,j)=>sub(z,approx[j])))};
}
export function mitigation(lambda,shots=1000,seed=1){
 // Z expectation of a known +1 state under an exponential decay model.
 const scales=[1,2,3],means=scales.map(s=>Math.exp(-lambda*s)),r=rng(seed);
 const estimates=means.map(e=>{const c=sample([(1+e)/2,(1-e)/2],shots,r);return (c[0]-c[1])/shots;});
 return {scales,means,estimates,linear:2*estimates[0]-estimates[1],quadratic:3*estimates[0]-3*estimates[1]+estimates[2],linearExact:2*means[0]-means[1],quadraticExact:3*means[0]-3*means[1]+means[2],varianceLinear:(4*(1-means[0]**2)+(1-means[1]**2))/shots};
}
