/* CausalScope #20: dependency-free statistical computation module. */
(function(root, factory) {
  const engine = factory();
  if (typeof module !== 'undefined' && module.exports) module.exports = engine;
  root.CausalCore = engine;
})(typeof globalThis !== 'undefined' ? globalThis : this, function() {
  'use strict';
  const clamp = (v,lo,hi) => Math.min(hi, Math.max(lo, v));
  const finite = x => typeof x === 'number' && Number.isFinite(x);
  // Abramowitz & Stegun 7.1.26, max absolute error ~1.5e-7 for erf.
  function erf(x) {
    const sign=x<0?-1:1, a=Math.abs(x), t=1/(1+0.3275911*a);
    const val=1-(((((1.061405429*t-1.453152027)*t)+1.421413741)*t-0.284496736)*t+0.254829592)*t*Math.exp(-a*a);
    return sign*val;
  }
  function normCDF(z){return 0.5*(1+erf(z/Math.SQRT2));}
  function normInv(p) {
    if (!(p>0&&p<1)) throw new Error('Normal quantile requires a probability between 0 and 1.');
    const a=[-39.69683028665376,220.9460984245205,-275.9285104469687,138.3577518672690,-30.66479806614716,2.506628277459239];
    const b=[-54.47609879822406,161.5858368580409,-155.6989798598866,66.80131188771972,-13.28068155288572];
    const c=[-0.007784894002430293,-0.3223964580411365,-2.400758277161838,-2.549732539343734,4.374664141464968,2.938163982698783];
    const d=[0.007784695709041462,0.3224671290700398,2.445134137142996,3.754408661907416];
    let q,r;
    if(p<0.02425){q=Math.sqrt(-2*Math.log(p));return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5])/((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);}
    if(p>0.97575){q=Math.sqrt(-2*Math.log(1-p));return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5])/((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);}
    q=p-0.5;r=q*q;
    return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q/(((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1);
  }
  function wilson(k,n,confidence=0.95){
    if(!Number.isInteger(k)||!Number.isInteger(n)||n<=0||k<0||k>n) throw new Error('Successes must be 0..n and sample size must be positive.');
    const z=normInv((1+confidence)/2), p=k/n, z2=z*z,denom=1+z2/n;
    const center=(p+z2/(2*n))/denom, margin=z*Math.sqrt((p*(1-p)+z2/(4*n))/n)/denom;
    return [clamp(center-margin,0,1),clamp(center+margin,0,1)];
  }
  function evaluate(experiment, alpha=0.05) {
    const nA=Number(experiment.nA), nB=Number(experiment.nB),kA=Number(experiment.kA),kB=Number(experiment.kB);
    if (![nA,nB,kA,kB].every(Number.isInteger)||nA<=0||nB<=0||kA<0||kB<0||kA>nA||kB>nB) throw new Error('Enter whole, nonnegative conversions no greater than their sample sizes.');
    if(!finite(alpha)||alpha<=0||alpha>=1)throw new Error('Alpha must be between 0 and 1.');
    const pA=kA/nA,pB=kB/nB,delta=pB-pA,pooled=(kA+kB)/(nA+nB),se=Math.sqrt(pooled*(1-pooled)*(1/nA+1/nB));
    const z=se>0?delta/se:0, pValue=se>0?clamp(2*(1-normCDF(Math.abs(z))),0,1):1;
    const conf=1-alpha,wa=wilson(kA,nA,conf),wb=wilson(kB,nB,conf);
    // Newcombe-Wilson interval for the difference of independent proportions.
    const low=delta-Math.sqrt((pB-wb[0])**2+(wa[1]-pA)**2);
    const high=delta+Math.sqrt((wb[1]-pB)**2+(pA-wa[0])**2);
    const expectedA=(nA*(kA+kB))/(nA+nB),expectedB=(nB*(kA+kB))/(nA+nB);
    const sparse=Math.min(expectedA,nA-expectedA,expectedB,nB-expectedB)<5;
    const srmZ=(nA-nB)/Math.sqrt(nA+nB),srmP=2*(1-normCDF(Math.abs(srmZ)));
    return {nA,nB,kA,kB,pA,pB,delta,relative:pA>0?delta/pA:null,confidence:conf,ci:[clamp(low,-1,1),clamp(high,-1,1)],z,pValue,detected:pValue<alpha,positive:delta>0,negative:delta<0,sparse,srmP,allocationWarning:srmP<0.001,absoluteConversionsAt10000:delta*10000};
  }
  function sampleSize(pA,relativeMDE,alpha=0.05,power=0.8) {
    if (![pA,relativeMDE,alpha,power].every(finite)||pA<=0||pA>=1||relativeMDE<=0||alpha<=0||alpha>=1||power<=0||power>=1)throw new Error('Check planning parameters.');
    const pB=pA*(1+relativeMDE);
    if(pB>=1)throw new Error('Expected treatment rate must be below 100%. Reduce the MDE.');
    const delta=pB-pA,avg=(pA+pB)/2,zA=normInv(1-alpha/2),zB=normInv(power);
    const n=Math.pow(zA*Math.sqrt(2*avg*(1-avg))+zB*Math.sqrt(pA*(1-pA)+pB*(1-pB)),2)/(delta*delta);
    return {perArm:Math.ceil(n),total:2*Math.ceil(n),pA,pB,absoluteMDE:delta,relativeMDE,alpha,power};
  }
  function seeded(seed){let state=(Number(seed)>>>0)||1;return ()=>{state=(1664525*state+1013904223)>>>0;return state/4294967296;};}
  function normal(random){return Math.sqrt(-2*Math.log(Math.max(1e-12,random())))*Math.cos(2*Math.PI*random());}
  // Normal approximations are intentionally labeled as simulations of planning assumptions, not real outcomes.
  function simulate(experiment, runs=1200,seed=202620){
    const {nA,nB,kA,kB}=experiment, rA=kA/nA,rB=kB/nB;
    if(!(nA>0&&nB>0))throw new Error('Set valid sample sizes.');
    const random=seeded(seed), deltas=[];let positive=0;
    for(let i=0;i<runs;i++){
      const x=clamp(rA+Math.sqrt(rA*(1-rA)/nA)*normal(random),0,1);
      const y=clamp(rB+Math.sqrt(rB*(1-rB)/nB)*normal(random),0,1);
      const d=y-x;deltas.push(d);if(d>0)positive++;
    }
    deltas.sort((a,b)=>a-b);
    return {runs,seed,sharePositive:positive/runs,p05:deltas[Math.floor(runs*0.05)],p50:deltas[Math.floor(runs*0.5)],p95:deltas[Math.floor(runs*0.95)],deltas};
  }
  return {erf,normCDF,normInv,wilson,evaluate,sampleSize,seeded,simulate};
});