<?php

declare(strict_types=1);

namespace Dochia;

/**
 * Normal/lognormal numeric helpers (no extensions, no dependencies).
 *
 * erf uses Abramowitz & Stegun 7.1.26 (|error| < 1.5e-7), which is ample for
 * day-level probabilities and quantiles; golden tests compare with the Python
 * reference at a tolerant epsilon.
 */
final class Erf
{
    public static function erf(float $x): float
    {
        $sign = $x < 0 ? -1.0 : 1.0;
        $x = abs($x);
        $a1 = 0.254829592;
        $a2 = -0.284496736;
        $a3 = 1.421413741;
        $a4 = -1.453152027;
        $a5 = 1.061405429;
        $p = 0.3275911;
        $t = 1.0 / (1.0 + $p * $x);
        $y = 1.0 - ((((($a5 * $t + $a4) * $t + $a3) * $t + $a2) * $t + $a1) * $t) * exp(-$x * $x);
        return $sign * $y;
    }

    public static function normCdf(float $z): float
    {
        return 0.5 * (1.0 + self::erf($z / M_SQRT2));
    }

    public static function normPdf(float $z): float
    {
        return exp(-0.5 * $z * $z) / sqrt(2.0 * M_PI);
    }

    /** Inverse normal CDF by bisection (matches the Python reference to ~1e-12). */
    public static function normPpf(float $p): float
    {
        if ($p <= 0.0) {
            return -INF;
        }
        if ($p >= 1.0) {
            return INF;
        }
        $lo = -12.0;
        $hi = 12.0;
        for ($i = 0; $i < 200; $i++) {
            $mid = 0.5 * ($lo + $hi);
            if (self::normCdf($mid) < $p) {
                $lo = $mid;
            } else {
                $hi = $mid;
            }
        }
        return 0.5 * ($lo + $hi);
    }
}