From 218cfac4c43072d1335f4ca15cba21bf971522c9 Mon Sep 17 00:00:00 2001 From: AminMohamed-3 <124639423+AminMohamed-3@users.noreply.github.com> Date: Wed, 9 Sep 2026 01:04:26 +0300 Subject: [PATCH 1/4] fix(stats): propagate NaN through normal CDF and quantile kernels --- src/flopscope/stats/_erf.py | 381 +++++++++++++++++----------------- src/flopscope/stats/_ndtri.py | 249 +++++++++++----------- 2 files changed, 316 insertions(+), 314 deletions(-) diff --git a/src/flopscope/stats/_erf.py b/src/flopscope/stats/_erf.py index fc222a807a..9b73760279 100644 --- a/src/flopscope/stats/_erf.py +++ b/src/flopscope/stats/_erf.py @@ -1,190 +1,191 @@ -"""Vectorized error function using the fdlibm/Sun rational approximation. - -Uses the same algorithm and coefficients as the C99 math.erf / glibc / fdlibm -``s_erf.c`` implementation, which is the gold-standard for double-precision erf. - -Four regions: - |x| < 0.84375 : erf(x) = x + x*P(x^2)/Q(x^2) - 0.84375<=|x|<1.25 : erf(x) = erx + P(|x|-1)/Q(|x|-1) - 1.25 <= |x| < ~6 : erfc(x) = exp(-x^2-0.5625+R(1/x^2)/S(1/x^2)) / |x| - |x| >= 6 : erf(x) = sign(x) * 1.0 - -Accuracy: matches C99 math.erf / scipy.special.erf to ~1 ULP. - -References ----------- -.. [1] Sun Microsystems, "Freely Distributable LIBM (fdlibm) s_erf.c", - https://www.netlib.org/fdlibm/s_erf.c - Algorithm and coefficients by Sun Microsystems (see copyright below). - -.. [2] glibc implementation (derived from fdlibm): - https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_erf.c - -Copyright (original fdlibm source) ------------------------------------ -Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. -Developed at SunSoft, a Sun Microsystems, Inc. business. -Permission to use, copy, modify, and distribute this software is freely -granted, provided that this notice is preserved. -""" - -from __future__ import annotations - -import numpy as np - -# --------------------------------------------------------------------------- -# Coefficients from fdlibm s_erf.c (Sun / glibc / OpenBSD) -# --------------------------------------------------------------------------- - -# erx = erf(1) to extended precision -_erx = 8.45062911510467529297e-01 - -# Region 1: |x| < 0.84375 -- erf(x) = x + x * pp(x^2)/qq(x^2) -_pp0 = 1.28379167095512558561e-01 -_pp1 = -3.25042107247001499370e-01 -_pp2 = -2.84817495755985104766e-02 -_pp3 = -5.77027029648944159157e-03 -_pp4 = -2.37630166566501626084e-05 - -_qq1 = 3.97917223959155352819e-01 -_qq2 = 6.50222499887672944485e-02 -_qq3 = 5.08130628187576562776e-03 -_qq4 = 1.32494738004321644526e-04 -_qq5 = -3.96022827877536812320e-06 - -# Region 2: 0.84375 <= |x| < 1.25 -- erf(x) = erx + P(s)/Q(s), s = |x| - 1 -_pa0 = -2.36211856075265944077e-03 -_pa1 = 4.14856118683748331666e-01 -_pa2 = -3.72207876035701323847e-01 -_pa3 = 3.18346619901161753674e-01 -_pa4 = -1.10894694282396677476e-01 -_pa5 = 3.54783043195201877747e-02 -_pa6 = -2.16637559983254089680e-03 - -_qa1 = 1.06420880400844228286e-01 -_qa2 = 5.40397917702171048937e-01 -_qa3 = 7.18286544141962539399e-02 -_qa4 = 1.26171219808761642112e-01 -_qa5 = 1.36370839120290507362e-02 -_qa6 = 1.19844998467991074170e-02 - -# Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) -_ra0 = -9.86494403484714822705e-03 -_ra1 = -6.93858572707181764372e-01 -_ra2 = -1.05586262253232909814e01 -_ra3 = -6.23753324503260060396e01 -_ra4 = -1.62396669462573071767e02 -_ra5 = -1.84605092906711035994e02 -_ra6 = -8.12874355063065934246e01 -_ra7 = -9.81432934416914548592e00 - -_sa1 = 1.96512716674392571292e01 -_sa2 = 1.37657754143519702237e02 -_sa3 = 4.34565877475229228608e02 -_sa4 = 6.45387271733267880594e02 -_sa5 = 4.29008140027567833386e02 -_sa6 = 1.08635005541779435134e02 -_sa7 = 6.57024977031928170135e00 -_sa8 = -6.04244152148580987438e-02 - -# Region 3b: 1/0.35 (~2.857) <= |x| < 6 -_rb0 = -9.86494292470009928597e-03 -_rb1 = -7.99283237680523006574e-01 -_rb2 = -1.77579549177547519889e01 -_rb3 = -1.60636384855557935030e02 -_rb4 = -6.37566443368389085394e02 -_rb5 = -1.02509513161107724954e03 -_rb6 = -4.83519191608651397019e02 - -_sb1 = 3.03380607875625778203e01 -_sb2 = 3.25792512996573918826e02 -_sb3 = 1.53672958608443695994e03 -_sb4 = 3.19985821950859553908e03 -_sb5 = 2.55305040643316442583e03 -_sb6 = 4.74528541206955367215e02 -_sb7 = -2.24409524465858183362e01 - - -def _erf(x): - """Vectorized error function matching scipy.special.erf to ~1 ULP. - - Algorithm and coefficients from fdlibm ``s_erf.c`` [1]_. - """ - x = np.asarray(x, dtype=np.float64) - scalar = x.ndim == 0 - x = np.atleast_1d(x) - - out = np.empty_like(x) - ax = np.abs(x) - sign = np.sign(x) - - # --- Region 1: |x| < 0.84375 --- - m1 = ax < 0.84375 - if np.any(m1): - xm = x[m1] - s = xm * xm - r = _pp0 + s * (_pp1 + s * (_pp2 + s * (_pp3 + s * _pp4))) - S = 1.0 + s * (_qq1 + s * (_qq2 + s * (_qq3 + s * (_qq4 + s * _qq5)))) - out[m1] = xm + xm * (r / S) - - # --- Region 2: 0.84375 <= |x| < 1.25 --- - m2 = (ax >= 0.84375) & (ax < 1.25) - if np.any(m2): - s = ax[m2] - 1.0 - P = _pa0 + s * ( - _pa1 + s * (_pa2 + s * (_pa3 + s * (_pa4 + s * (_pa5 + s * _pa6)))) - ) - Q = 1.0 + s * ( - _qa1 + s * (_qa2 + s * (_qa3 + s * (_qa4 + s * (_qa5 + s * _qa6)))) - ) - out[m2] = sign[m2] * (_erx + P / Q) - - # --- Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) --- - m3a = (ax >= 1.25) & (ax < (1.0 / 0.35)) - if np.any(m3a): - axm = ax[m3a] - s = 1.0 / (axm * axm) - R = _ra0 + s * ( - _ra1 - + s * (_ra2 + s * (_ra3 + s * (_ra4 + s * (_ra5 + s * (_ra6 + s * _ra7))))) - ) - S = 1.0 + s * ( - _sa1 - + s - * ( - _sa2 - + s - * (_sa3 + s * (_sa4 + s * (_sa5 + s * (_sa6 + s * (_sa7 + s * _sa8))))) - ) - ) - erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm - out[m3a] = sign[m3a] * (1.0 - erfc_val) - - # --- Region 3b: 1/0.35 <= |x| < 6 --- - m3b = (ax >= (1.0 / 0.35)) & (ax < 6.0) - if np.any(m3b): - axm = ax[m3b] - s = 1.0 / (axm * axm) - R = _rb0 + s * ( - _rb1 + s * (_rb2 + s * (_rb3 + s * (_rb4 + s * (_rb5 + s * _rb6)))) - ) - S = 1.0 + s * ( - _sb1 - + s * (_sb2 + s * (_sb3 + s * (_sb4 + s * (_sb5 + s * (_sb6 + s * _sb7))))) - ) - erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm - out[m3b] = sign[m3b] * (1.0 - erfc_val) - - # --- Region 4: |x| >= 6 --- - m4 = ax >= 6.0 - if np.any(m4): - out[m4] = sign[m4] * 1.0 - - if scalar: - return float(out[0]) - return out - - -def _erfc(x): - """Vectorized complementary error function: erfc(x) = 1 - erf(x).""" - return 1.0 - _erf(x) +"""Vectorized error function using the fdlibm/Sun rational approximation. + +Uses the same algorithm and coefficients as the C99 math.erf / glibc / fdlibm +``s_erf.c`` implementation, which is the gold-standard for double-precision erf. + +Four regions: + |x| < 0.84375 : erf(x) = x + x*P(x^2)/Q(x^2) + 0.84375<=|x|<1.25 : erf(x) = erx + P(|x|-1)/Q(|x|-1) + 1.25 <= |x| < ~6 : erfc(x) = exp(-x^2-0.5625+R(1/x^2)/S(1/x^2)) / |x| + |x| >= 6 : erf(x) = sign(x) * 1.0 + +Accuracy: matches C99 math.erf / scipy.special.erf to ~1 ULP. + +References +---------- +.. [1] Sun Microsystems, "Freely Distributable LIBM (fdlibm) s_erf.c", + https://www.netlib.org/fdlibm/s_erf.c + Algorithm and coefficients by Sun Microsystems (see copyright below). + +.. [2] glibc implementation (derived from fdlibm): + https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_erf.c + +Copyright (original fdlibm source) +----------------------------------- +Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +Developed at SunSoft, a Sun Microsystems, Inc. business. +Permission to use, copy, modify, and distribute this software is freely +granted, provided that this notice is preserved. +""" + +from __future__ import annotations + +import numpy as np + +# --------------------------------------------------------------------------- +# Coefficients from fdlibm s_erf.c (Sun / glibc / OpenBSD) +# --------------------------------------------------------------------------- + +# erx = erf(1) to extended precision +_erx = 8.45062911510467529297e-01 + +# Region 1: |x| < 0.84375 -- erf(x) = x + x * pp(x^2)/qq(x^2) +_pp0 = 1.28379167095512558561e-01 +_pp1 = -3.25042107247001499370e-01 +_pp2 = -2.84817495755985104766e-02 +_pp3 = -5.77027029648944159157e-03 +_pp4 = -2.37630166566501626084e-05 + +_qq1 = 3.97917223959155352819e-01 +_qq2 = 6.50222499887672944485e-02 +_qq3 = 5.08130628187576562776e-03 +_qq4 = 1.32494738004321644526e-04 +_qq5 = -3.96022827877536812320e-06 + +# Region 2: 0.84375 <= |x| < 1.25 -- erf(x) = erx + P(s)/Q(s), s = |x| - 1 +_pa0 = -2.36211856075265944077e-03 +_pa1 = 4.14856118683748331666e-01 +_pa2 = -3.72207876035701323847e-01 +_pa3 = 3.18346619901161753674e-01 +_pa4 = -1.10894694282396677476e-01 +_pa5 = 3.54783043195201877747e-02 +_pa6 = -2.16637559983254089680e-03 + +_qa1 = 1.06420880400844228286e-01 +_qa2 = 5.40397917702171048937e-01 +_qa3 = 7.18286544141962539399e-02 +_qa4 = 1.26171219808761642112e-01 +_qa5 = 1.36370839120290507362e-02 +_qa6 = 1.19844998467991074170e-02 + +# Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) +_ra0 = -9.86494403484714822705e-03 +_ra1 = -6.93858572707181764372e-01 +_ra2 = -1.05586262253232909814e01 +_ra3 = -6.23753324503260060396e01 +_ra4 = -1.62396669462573071767e02 +_ra5 = -1.84605092906711035994e02 +_ra6 = -8.12874355063065934246e01 +_ra7 = -9.81432934416914548592e00 + +_sa1 = 1.96512716674392571292e01 +_sa2 = 1.37657754143519702237e02 +_sa3 = 4.34565877475229228608e02 +_sa4 = 6.45387271733267880594e02 +_sa5 = 4.29008140027567833386e02 +_sa6 = 1.08635005541779435134e02 +_sa7 = 6.57024977031928170135e00 +_sa8 = -6.04244152148580987438e-02 + +# Region 3b: 1/0.35 (~2.857) <= |x| < 6 +_rb0 = -9.86494292470009928597e-03 +_rb1 = -7.99283237680523006574e-01 +_rb2 = -1.77579549177547519889e01 +_rb3 = -1.60636384855557935030e02 +_rb4 = -6.37566443368389085394e02 +_rb5 = -1.02509513161107724954e03 +_rb6 = -4.83519191608651397019e02 + +_sb1 = 3.03380607875625778203e01 +_sb2 = 3.25792512996573918826e02 +_sb3 = 1.53672958608443695994e03 +_sb4 = 3.19985821950859553908e03 +_sb5 = 2.55305040643316442583e03 +_sb6 = 4.74528541206955367215e02 +_sb7 = -2.24409524465858183362e01 + + +def _erf(x): + """Vectorized error function matching scipy.special.erf to ~1 ULP. + + Algorithm and coefficients from fdlibm ``s_erf.c`` [1]_. + """ + x = np.asarray(x, dtype=np.float64) + scalar = x.ndim == 0 + x = np.atleast_1d(x) + + # NaN matches none of the region masks below, so initialize its result. + out = np.full_like(x, np.nan) + ax = np.abs(x) + sign = np.sign(x) + + # --- Region 1: |x| < 0.84375 --- + m1 = ax < 0.84375 + if np.any(m1): + xm = x[m1] + s = xm * xm + r = _pp0 + s * (_pp1 + s * (_pp2 + s * (_pp3 + s * _pp4))) + S = 1.0 + s * (_qq1 + s * (_qq2 + s * (_qq3 + s * (_qq4 + s * _qq5)))) + out[m1] = xm + xm * (r / S) + + # --- Region 2: 0.84375 <= |x| < 1.25 --- + m2 = (ax >= 0.84375) & (ax < 1.25) + if np.any(m2): + s = ax[m2] - 1.0 + P = _pa0 + s * ( + _pa1 + s * (_pa2 + s * (_pa3 + s * (_pa4 + s * (_pa5 + s * _pa6)))) + ) + Q = 1.0 + s * ( + _qa1 + s * (_qa2 + s * (_qa3 + s * (_qa4 + s * (_qa5 + s * _qa6)))) + ) + out[m2] = sign[m2] * (_erx + P / Q) + + # --- Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) --- + m3a = (ax >= 1.25) & (ax < (1.0 / 0.35)) + if np.any(m3a): + axm = ax[m3a] + s = 1.0 / (axm * axm) + R = _ra0 + s * ( + _ra1 + + s * (_ra2 + s * (_ra3 + s * (_ra4 + s * (_ra5 + s * (_ra6 + s * _ra7))))) + ) + S = 1.0 + s * ( + _sa1 + + s + * ( + _sa2 + + s + * (_sa3 + s * (_sa4 + s * (_sa5 + s * (_sa6 + s * (_sa7 + s * _sa8))))) + ) + ) + erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm + out[m3a] = sign[m3a] * (1.0 - erfc_val) + + # --- Region 3b: 1/0.35 <= |x| < 6 --- + m3b = (ax >= (1.0 / 0.35)) & (ax < 6.0) + if np.any(m3b): + axm = ax[m3b] + s = 1.0 / (axm * axm) + R = _rb0 + s * ( + _rb1 + s * (_rb2 + s * (_rb3 + s * (_rb4 + s * (_rb5 + s * _rb6)))) + ) + S = 1.0 + s * ( + _sb1 + + s * (_sb2 + s * (_sb3 + s * (_sb4 + s * (_sb5 + s * (_sb6 + s * _sb7))))) + ) + erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm + out[m3b] = sign[m3b] * (1.0 - erfc_val) + + # --- Region 4: |x| >= 6 --- + m4 = ax >= 6.0 + if np.any(m4): + out[m4] = sign[m4] * 1.0 + + if scalar: + return float(out[0]) + return out + + +def _erfc(x): + """Vectorized complementary error function: erfc(x) = 1 - erf(x).""" + return 1.0 - _erf(x) diff --git a/src/flopscope/stats/_ndtri.py b/src/flopscope/stats/_ndtri.py index 44321d87e7..9740369e6b 100644 --- a/src/flopscope/stats/_ndtri.py +++ b/src/flopscope/stats/_ndtri.py @@ -1,124 +1,125 @@ -"""Inverse standard normal CDF (ndtri) via a rational approximation + Newton -refinement. - -Accuracy: ~1e-12 against scipy.special.ndtri. -""" - -from __future__ import annotations - -import numpy as np - -from flopscope.stats._erf import _erf - -# --------------------------------------------------------------------------- -# Rational approximation coefficients -# --------------------------------------------------------------------------- - -_A = ( - -3.969683028665376e01, - 2.209460984245205e02, - -2.759285104469687e02, - 1.383577518672690e02, - -3.066479806614716e01, - 2.506628277459239e00, -) -_B = ( - -5.447609879822406e01, - 1.615858368580409e02, - -1.556989798598866e02, - 6.680131188771972e01, - -1.328068155288572e01, -) -_C = ( - -7.784894002430293e-03, - -3.223964580411365e-01, - -2.400758277161838e00, - -2.549732539343734e00, - 4.374664141464968e00, - 2.938163982698783e00, -) -_D = ( - 7.784695709041462e-03, - 3.224671290700398e-01, - 2.445134137142996e00, - 3.754408661907416e00, -) - -_P_LOW = 0.02425 -_P_HIGH = 1.0 - _P_LOW - -_SQRT2 = np.sqrt(2.0) -_INV_SQRT_2PI = 1.0 / np.sqrt(2.0 * np.pi) - - -def _norm_pdf_internal(x): - """Standard normal PDF (internal helper, no budget deduction).""" - return _INV_SQRT_2PI * np.exp(-0.5 * x * x) - - -def _norm_cdf_internal(x): - """Standard normal CDF (internal helper, no budget deduction).""" - return 0.5 * (1.0 + _erf(x / _SQRT2)) - - -def _ndtri(p): - """Inverse standard normal CDF. - - Maps probability *p* in (0, 1) to the quantile *x* such that - Phi(x) = p, where Phi is the standard normal CDF. - - Uses a rational approximation with one Newton-Raphson refinement step - for ~1e-12 accuracy. - - Edge cases: p=0 -> -inf, p=1 -> +inf, p<0 or p>1 -> nan. - """ - p = np.asarray(p, dtype=np.float64) - scalar = p.ndim == 0 - p = np.atleast_1d(p) - - out = np.empty_like(p) - - # Edge cases - out[p == 0.0] = -np.inf - out[p == 1.0] = np.inf - out[(p < 0.0) | (p > 1.0)] = np.nan - - # Lower region: 0 < p < P_LOW - m_low = (p > 0.0) & (p < _P_LOW) - if np.any(m_low): - q = np.sqrt(-2.0 * np.log(p[m_low])) - num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] - den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 - out[m_low] = num / den - - # Central region: P_LOW <= p <= P_HIGH - m_mid = (p >= _P_LOW) & (p <= _P_HIGH) - if np.any(m_mid): - q = p[m_mid] - 0.5 - r = q * q - num = ((((_A[0] * r + _A[1]) * r + _A[2]) * r + _A[3]) * r + _A[4]) * r + _A[5] - den = ((((_B[0] * r + _B[1]) * r + _B[2]) * r + _B[3]) * r + _B[4]) * r + 1.0 - out[m_mid] = q * num / den - - # Upper region: P_HIGH < p < 1 - m_high = (p > _P_HIGH) & (p < 1.0) - if np.any(m_high): - q = np.sqrt(-2.0 * np.log(1.0 - p[m_high])) - num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] - den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 - out[m_high] = -(num / den) - - # Newton refinement (one step) for the interior - m_interior = (p > 0.0) & (p < 1.0) - if np.any(m_interior): - x0 = out[m_interior] - phi = _norm_cdf_internal(x0) - pdf = _norm_pdf_internal(x0) - # Avoid division by zero in pdf - safe = pdf > 1e-300 - correction = np.where(safe, (phi - p[m_interior]) / pdf, 0.0) - out[m_interior] = x0 - correction - - if scalar: - return float(out[0]) - return out +"""Inverse standard normal CDF (ndtri) via a rational approximation + Newton +refinement. + +Accuracy: ~1e-12 against scipy.special.ndtri. +""" + +from __future__ import annotations + +import numpy as np + +from flopscope.stats._erf import _erf + +# --------------------------------------------------------------------------- +# Rational approximation coefficients +# --------------------------------------------------------------------------- + +_A = ( + -3.969683028665376e01, + 2.209460984245205e02, + -2.759285104469687e02, + 1.383577518672690e02, + -3.066479806614716e01, + 2.506628277459239e00, +) +_B = ( + -5.447609879822406e01, + 1.615858368580409e02, + -1.556989798598866e02, + 6.680131188771972e01, + -1.328068155288572e01, +) +_C = ( + -7.784894002430293e-03, + -3.223964580411365e-01, + -2.400758277161838e00, + -2.549732539343734e00, + 4.374664141464968e00, + 2.938163982698783e00, +) +_D = ( + 7.784695709041462e-03, + 3.224671290700398e-01, + 2.445134137142996e00, + 3.754408661907416e00, +) + +_P_LOW = 0.02425 +_P_HIGH = 1.0 - _P_LOW + +_SQRT2 = np.sqrt(2.0) +_INV_SQRT_2PI = 1.0 / np.sqrt(2.0 * np.pi) + + +def _norm_pdf_internal(x): + """Standard normal PDF (internal helper, no budget deduction).""" + return _INV_SQRT_2PI * np.exp(-0.5 * x * x) + + +def _norm_cdf_internal(x): + """Standard normal CDF (internal helper, no budget deduction).""" + return 0.5 * (1.0 + _erf(x / _SQRT2)) + + +def _ndtri(p): + """Inverse standard normal CDF. + + Maps probability *p* in (0, 1) to the quantile *x* such that + Phi(x) = p, where Phi is the standard normal CDF. + + Uses a rational approximation with one Newton-Raphson refinement step + for ~1e-12 accuracy. + + Edge cases: p=0 -> -inf, p=1 -> +inf, p<0 or p>1 -> nan. + """ + p = np.asarray(p, dtype=np.float64) + scalar = p.ndim == 0 + p = np.atleast_1d(p) + + # NaN matches neither the edge cases nor an approximation region. + out = np.full_like(p, np.nan) + + # Edge cases + out[p == 0.0] = -np.inf + out[p == 1.0] = np.inf + out[(p < 0.0) | (p > 1.0)] = np.nan + + # Lower region: 0 < p < P_LOW + m_low = (p > 0.0) & (p < _P_LOW) + if np.any(m_low): + q = np.sqrt(-2.0 * np.log(p[m_low])) + num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] + den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 + out[m_low] = num / den + + # Central region: P_LOW <= p <= P_HIGH + m_mid = (p >= _P_LOW) & (p <= _P_HIGH) + if np.any(m_mid): + q = p[m_mid] - 0.5 + r = q * q + num = ((((_A[0] * r + _A[1]) * r + _A[2]) * r + _A[3]) * r + _A[4]) * r + _A[5] + den = ((((_B[0] * r + _B[1]) * r + _B[2]) * r + _B[3]) * r + _B[4]) * r + 1.0 + out[m_mid] = q * num / den + + # Upper region: P_HIGH < p < 1 + m_high = (p > _P_HIGH) & (p < 1.0) + if np.any(m_high): + q = np.sqrt(-2.0 * np.log(1.0 - p[m_high])) + num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] + den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 + out[m_high] = -(num / den) + + # Newton refinement (one step) for the interior + m_interior = (p > 0.0) & (p < 1.0) + if np.any(m_interior): + x0 = out[m_interior] + phi = _norm_cdf_internal(x0) + pdf = _norm_pdf_internal(x0) + # Avoid division by zero in pdf + safe = pdf > 1e-300 + correction = np.where(safe, (phi - p[m_interior]) / pdf, 0.0) + out[m_interior] = x0 - correction + + if scalar: + return float(out[0]) + return out From 6bb33aced47122f00595c0863debcdf5fe31f023 Mon Sep 17 00:00:00 2001 From: AminMohamed-3 <124639423+AminMohamed-3@users.noreply.github.com> Date: Wed, 9 Sep 2026 01:06:04 +0300 Subject: [PATCH 2/4] test(stats): cover NaN propagation and unchanged FLOP charges --- tests/test_stats_nan_propagation.py | 79 +++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) create mode 100644 tests/test_stats_nan_propagation.py diff --git a/tests/test_stats_nan_propagation.py b/tests/test_stats_nan_propagation.py new file mode 100644 index 0000000000..3907ade023 --- /dev/null +++ b/tests/test_stats_nan_propagation.py @@ -0,0 +1,79 @@ +"""NaN propagation must not depend on recycled allocation contents.""" + +import numpy as np +import pytest +from scipy import special, stats + +from flopscope import BudgetContext +from flopscope.stats import lognorm, norm, truncnorm +from flopscope.stats._erf import _erf +from flopscope.stats._ndtri import _ndtri + + +@pytest.fixture +def finite_empty_buffers(monkeypatch): + """Make uninitialized outputs deterministic rather than allocator-dependent. + + Filling empty buffers with a finite sentinel is valid: empty_like makes + no promise about its contents. This avoids false passes when recycled + buffers happen to contain NaN already. + """ + original = np.empty_like + + def empty_with_sentinel(*args, **kwargs): + result = original(*args, **kwargs) + result.fill(0.125) + return result + + monkeypatch.setattr(np, "empty_like", empty_with_sentinel) + + +@pytest.mark.parametrize( + "actual,expected,values", + [ + (_erf, special.erf, [-np.inf, -2, 0, np.nan, 1, np.inf]), + (_ndtri, special.ndtri, [-0.1, 0, 0.25, np.nan, 0.75, 1, 1.1]), + (norm.cdf, stats.norm.cdf, [-np.inf, -2, 0, np.nan, 1, np.inf]), + (norm.ppf, stats.norm.ppf, [-0.1, 0, 0.25, np.nan, 0.75, 1, 1.1]), + ( + lambda x: lognorm.ppf(x, 0.5), + lambda x: stats.lognorm.ppf(x, 0.5), + [0, 0.25, np.nan, 0.75, 1], + ), + ( + lambda x: truncnorm.cdf(x, -2, 2), + lambda x: stats.truncnorm.cdf(x, -2, 2), + [-3, -1, np.nan, 1, 3], + ), + ( + lambda x: truncnorm.ppf(x, -2, 2), + lambda x: stats.truncnorm.ppf(x, -2, 2), + [0, 0.25, np.nan, 0.75, 1], + ), + ], +) +def test_nan_propagates_with_mixed_valid_values( + finite_empty_buffers, actual, expected, values +): + values = np.asarray(values, dtype=np.float64) + np.testing.assert_allclose( + np.asarray(actual(values)), + expected(values), + rtol=1e-11, + atol=1e-11, + equal_nan=True, + ) + + +@pytest.mark.parametrize("fn", [_erf, _ndtri, norm.cdf, norm.ppf]) +def test_scalar_nan_propagates(finite_empty_buffers, fn): + assert np.isnan(np.asarray(fn(np.nan))) + + +@pytest.mark.parametrize("fn", [norm.cdf, norm.ppf]) +def test_nan_input_keeps_same_shape_flop_cost(fn): + with BudgetContext(flop_budget=10**6) as finite_budget: + fn(np.array([0.25, 0.5, 0.75])) + with BudgetContext(flop_budget=10**6) as nan_budget: + fn(np.array([0.25, np.nan, 0.75])) + assert nan_budget.flops_used == finite_budget.flops_used > 0 From 9c9a1fdfe59fe3786672cf6d1b7e71f265764ea8 Mon Sep 17 00:00:00 2001 From: AminMohamed-3 <124639423+AminMohamed-3@users.noreply.github.com> Date: Wed, 9 Sep 2026 01:07:12 +0300 Subject: [PATCH 3/4] style(stats): preserve upstream LF line endings --- src/flopscope/stats/_erf.py | 382 +++++++++++++++++----------------- src/flopscope/stats/_ndtri.py | 250 +++++++++++----------- 2 files changed, 316 insertions(+), 316 deletions(-) diff --git a/src/flopscope/stats/_erf.py b/src/flopscope/stats/_erf.py index 9b73760279..f1cb5ec751 100644 --- a/src/flopscope/stats/_erf.py +++ b/src/flopscope/stats/_erf.py @@ -1,191 +1,191 @@ -"""Vectorized error function using the fdlibm/Sun rational approximation. - -Uses the same algorithm and coefficients as the C99 math.erf / glibc / fdlibm -``s_erf.c`` implementation, which is the gold-standard for double-precision erf. - -Four regions: - |x| < 0.84375 : erf(x) = x + x*P(x^2)/Q(x^2) - 0.84375<=|x|<1.25 : erf(x) = erx + P(|x|-1)/Q(|x|-1) - 1.25 <= |x| < ~6 : erfc(x) = exp(-x^2-0.5625+R(1/x^2)/S(1/x^2)) / |x| - |x| >= 6 : erf(x) = sign(x) * 1.0 - -Accuracy: matches C99 math.erf / scipy.special.erf to ~1 ULP. - -References ----------- -.. [1] Sun Microsystems, "Freely Distributable LIBM (fdlibm) s_erf.c", - https://www.netlib.org/fdlibm/s_erf.c - Algorithm and coefficients by Sun Microsystems (see copyright below). - -.. [2] glibc implementation (derived from fdlibm): - https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_erf.c - -Copyright (original fdlibm source) ------------------------------------ -Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. -Developed at SunSoft, a Sun Microsystems, Inc. business. -Permission to use, copy, modify, and distribute this software is freely -granted, provided that this notice is preserved. -""" - -from __future__ import annotations - -import numpy as np - -# --------------------------------------------------------------------------- -# Coefficients from fdlibm s_erf.c (Sun / glibc / OpenBSD) -# --------------------------------------------------------------------------- - -# erx = erf(1) to extended precision -_erx = 8.45062911510467529297e-01 - -# Region 1: |x| < 0.84375 -- erf(x) = x + x * pp(x^2)/qq(x^2) -_pp0 = 1.28379167095512558561e-01 -_pp1 = -3.25042107247001499370e-01 -_pp2 = -2.84817495755985104766e-02 -_pp3 = -5.77027029648944159157e-03 -_pp4 = -2.37630166566501626084e-05 - -_qq1 = 3.97917223959155352819e-01 -_qq2 = 6.50222499887672944485e-02 -_qq3 = 5.08130628187576562776e-03 -_qq4 = 1.32494738004321644526e-04 -_qq5 = -3.96022827877536812320e-06 - -# Region 2: 0.84375 <= |x| < 1.25 -- erf(x) = erx + P(s)/Q(s), s = |x| - 1 -_pa0 = -2.36211856075265944077e-03 -_pa1 = 4.14856118683748331666e-01 -_pa2 = -3.72207876035701323847e-01 -_pa3 = 3.18346619901161753674e-01 -_pa4 = -1.10894694282396677476e-01 -_pa5 = 3.54783043195201877747e-02 -_pa6 = -2.16637559983254089680e-03 - -_qa1 = 1.06420880400844228286e-01 -_qa2 = 5.40397917702171048937e-01 -_qa3 = 7.18286544141962539399e-02 -_qa4 = 1.26171219808761642112e-01 -_qa5 = 1.36370839120290507362e-02 -_qa6 = 1.19844998467991074170e-02 - -# Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) -_ra0 = -9.86494403484714822705e-03 -_ra1 = -6.93858572707181764372e-01 -_ra2 = -1.05586262253232909814e01 -_ra3 = -6.23753324503260060396e01 -_ra4 = -1.62396669462573071767e02 -_ra5 = -1.84605092906711035994e02 -_ra6 = -8.12874355063065934246e01 -_ra7 = -9.81432934416914548592e00 - -_sa1 = 1.96512716674392571292e01 -_sa2 = 1.37657754143519702237e02 -_sa3 = 4.34565877475229228608e02 -_sa4 = 6.45387271733267880594e02 -_sa5 = 4.29008140027567833386e02 -_sa6 = 1.08635005541779435134e02 -_sa7 = 6.57024977031928170135e00 -_sa8 = -6.04244152148580987438e-02 - -# Region 3b: 1/0.35 (~2.857) <= |x| < 6 -_rb0 = -9.86494292470009928597e-03 -_rb1 = -7.99283237680523006574e-01 -_rb2 = -1.77579549177547519889e01 -_rb3 = -1.60636384855557935030e02 -_rb4 = -6.37566443368389085394e02 -_rb5 = -1.02509513161107724954e03 -_rb6 = -4.83519191608651397019e02 - -_sb1 = 3.03380607875625778203e01 -_sb2 = 3.25792512996573918826e02 -_sb3 = 1.53672958608443695994e03 -_sb4 = 3.19985821950859553908e03 -_sb5 = 2.55305040643316442583e03 -_sb6 = 4.74528541206955367215e02 -_sb7 = -2.24409524465858183362e01 - - -def _erf(x): - """Vectorized error function matching scipy.special.erf to ~1 ULP. - - Algorithm and coefficients from fdlibm ``s_erf.c`` [1]_. - """ - x = np.asarray(x, dtype=np.float64) - scalar = x.ndim == 0 - x = np.atleast_1d(x) - - # NaN matches none of the region masks below, so initialize its result. - out = np.full_like(x, np.nan) - ax = np.abs(x) - sign = np.sign(x) - - # --- Region 1: |x| < 0.84375 --- - m1 = ax < 0.84375 - if np.any(m1): - xm = x[m1] - s = xm * xm - r = _pp0 + s * (_pp1 + s * (_pp2 + s * (_pp3 + s * _pp4))) - S = 1.0 + s * (_qq1 + s * (_qq2 + s * (_qq3 + s * (_qq4 + s * _qq5)))) - out[m1] = xm + xm * (r / S) - - # --- Region 2: 0.84375 <= |x| < 1.25 --- - m2 = (ax >= 0.84375) & (ax < 1.25) - if np.any(m2): - s = ax[m2] - 1.0 - P = _pa0 + s * ( - _pa1 + s * (_pa2 + s * (_pa3 + s * (_pa4 + s * (_pa5 + s * _pa6)))) - ) - Q = 1.0 + s * ( - _qa1 + s * (_qa2 + s * (_qa3 + s * (_qa4 + s * (_qa5 + s * _qa6)))) - ) - out[m2] = sign[m2] * (_erx + P / Q) - - # --- Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) --- - m3a = (ax >= 1.25) & (ax < (1.0 / 0.35)) - if np.any(m3a): - axm = ax[m3a] - s = 1.0 / (axm * axm) - R = _ra0 + s * ( - _ra1 - + s * (_ra2 + s * (_ra3 + s * (_ra4 + s * (_ra5 + s * (_ra6 + s * _ra7))))) - ) - S = 1.0 + s * ( - _sa1 - + s - * ( - _sa2 - + s - * (_sa3 + s * (_sa4 + s * (_sa5 + s * (_sa6 + s * (_sa7 + s * _sa8))))) - ) - ) - erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm - out[m3a] = sign[m3a] * (1.0 - erfc_val) - - # --- Region 3b: 1/0.35 <= |x| < 6 --- - m3b = (ax >= (1.0 / 0.35)) & (ax < 6.0) - if np.any(m3b): - axm = ax[m3b] - s = 1.0 / (axm * axm) - R = _rb0 + s * ( - _rb1 + s * (_rb2 + s * (_rb3 + s * (_rb4 + s * (_rb5 + s * _rb6)))) - ) - S = 1.0 + s * ( - _sb1 - + s * (_sb2 + s * (_sb3 + s * (_sb4 + s * (_sb5 + s * (_sb6 + s * _sb7))))) - ) - erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm - out[m3b] = sign[m3b] * (1.0 - erfc_val) - - # --- Region 4: |x| >= 6 --- - m4 = ax >= 6.0 - if np.any(m4): - out[m4] = sign[m4] * 1.0 - - if scalar: - return float(out[0]) - return out - - -def _erfc(x): - """Vectorized complementary error function: erfc(x) = 1 - erf(x).""" - return 1.0 - _erf(x) +"""Vectorized error function using the fdlibm/Sun rational approximation. + +Uses the same algorithm and coefficients as the C99 math.erf / glibc / fdlibm +``s_erf.c`` implementation, which is the gold-standard for double-precision erf. + +Four regions: + |x| < 0.84375 : erf(x) = x + x*P(x^2)/Q(x^2) + 0.84375<=|x|<1.25 : erf(x) = erx + P(|x|-1)/Q(|x|-1) + 1.25 <= |x| < ~6 : erfc(x) = exp(-x^2-0.5625+R(1/x^2)/S(1/x^2)) / |x| + |x| >= 6 : erf(x) = sign(x) * 1.0 + +Accuracy: matches C99 math.erf / scipy.special.erf to ~1 ULP. + +References +---------- +.. [1] Sun Microsystems, "Freely Distributable LIBM (fdlibm) s_erf.c", + https://www.netlib.org/fdlibm/s_erf.c + Algorithm and coefficients by Sun Microsystems (see copyright below). + +.. [2] glibc implementation (derived from fdlibm): + https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_erf.c + +Copyright (original fdlibm source) +----------------------------------- +Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +Developed at SunSoft, a Sun Microsystems, Inc. business. +Permission to use, copy, modify, and distribute this software is freely +granted, provided that this notice is preserved. +""" + +from __future__ import annotations + +import numpy as np + +# --------------------------------------------------------------------------- +# Coefficients from fdlibm s_erf.c (Sun / glibc / OpenBSD) +# --------------------------------------------------------------------------- + +# erx = erf(1) to extended precision +_erx = 8.45062911510467529297e-01 + +# Region 1: |x| < 0.84375 -- erf(x) = x + x * pp(x^2)/qq(x^2) +_pp0 = 1.28379167095512558561e-01 +_pp1 = -3.25042107247001499370e-01 +_pp2 = -2.84817495755985104766e-02 +_pp3 = -5.77027029648944159157e-03 +_pp4 = -2.37630166566501626084e-05 + +_qq1 = 3.97917223959155352819e-01 +_qq2 = 6.50222499887672944485e-02 +_qq3 = 5.08130628187576562776e-03 +_qq4 = 1.32494738004321644526e-04 +_qq5 = -3.96022827877536812320e-06 + +# Region 2: 0.84375 <= |x| < 1.25 -- erf(x) = erx + P(s)/Q(s), s = |x| - 1 +_pa0 = -2.36211856075265944077e-03 +_pa1 = 4.14856118683748331666e-01 +_pa2 = -3.72207876035701323847e-01 +_pa3 = 3.18346619901161753674e-01 +_pa4 = -1.10894694282396677476e-01 +_pa5 = 3.54783043195201877747e-02 +_pa6 = -2.16637559983254089680e-03 + +_qa1 = 1.06420880400844228286e-01 +_qa2 = 5.40397917702171048937e-01 +_qa3 = 7.18286544141962539399e-02 +_qa4 = 1.26171219808761642112e-01 +_qa5 = 1.36370839120290507362e-02 +_qa6 = 1.19844998467991074170e-02 + +# Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) +_ra0 = -9.86494403484714822705e-03 +_ra1 = -6.93858572707181764372e-01 +_ra2 = -1.05586262253232909814e01 +_ra3 = -6.23753324503260060396e01 +_ra4 = -1.62396669462573071767e02 +_ra5 = -1.84605092906711035994e02 +_ra6 = -8.12874355063065934246e01 +_ra7 = -9.81432934416914548592e00 + +_sa1 = 1.96512716674392571292e01 +_sa2 = 1.37657754143519702237e02 +_sa3 = 4.34565877475229228608e02 +_sa4 = 6.45387271733267880594e02 +_sa5 = 4.29008140027567833386e02 +_sa6 = 1.08635005541779435134e02 +_sa7 = 6.57024977031928170135e00 +_sa8 = -6.04244152148580987438e-02 + +# Region 3b: 1/0.35 (~2.857) <= |x| < 6 +_rb0 = -9.86494292470009928597e-03 +_rb1 = -7.99283237680523006574e-01 +_rb2 = -1.77579549177547519889e01 +_rb3 = -1.60636384855557935030e02 +_rb4 = -6.37566443368389085394e02 +_rb5 = -1.02509513161107724954e03 +_rb6 = -4.83519191608651397019e02 + +_sb1 = 3.03380607875625778203e01 +_sb2 = 3.25792512996573918826e02 +_sb3 = 1.53672958608443695994e03 +_sb4 = 3.19985821950859553908e03 +_sb5 = 2.55305040643316442583e03 +_sb6 = 4.74528541206955367215e02 +_sb7 = -2.24409524465858183362e01 + + +def _erf(x): + """Vectorized error function matching scipy.special.erf to ~1 ULP. + + Algorithm and coefficients from fdlibm ``s_erf.c`` [1]_. + """ + x = np.asarray(x, dtype=np.float64) + scalar = x.ndim == 0 + x = np.atleast_1d(x) + + # NaN matches none of the region masks below, so initialize its result. + out = np.full_like(x, np.nan) + ax = np.abs(x) + sign = np.sign(x) + + # --- Region 1: |x| < 0.84375 --- + m1 = ax < 0.84375 + if np.any(m1): + xm = x[m1] + s = xm * xm + r = _pp0 + s * (_pp1 + s * (_pp2 + s * (_pp3 + s * _pp4))) + S = 1.0 + s * (_qq1 + s * (_qq2 + s * (_qq3 + s * (_qq4 + s * _qq5)))) + out[m1] = xm + xm * (r / S) + + # --- Region 2: 0.84375 <= |x| < 1.25 --- + m2 = (ax >= 0.84375) & (ax < 1.25) + if np.any(m2): + s = ax[m2] - 1.0 + P = _pa0 + s * ( + _pa1 + s * (_pa2 + s * (_pa3 + s * (_pa4 + s * (_pa5 + s * _pa6)))) + ) + Q = 1.0 + s * ( + _qa1 + s * (_qa2 + s * (_qa3 + s * (_qa4 + s * (_qa5 + s * _qa6)))) + ) + out[m2] = sign[m2] * (_erx + P / Q) + + # --- Region 3a: 1.25 <= |x| < 1/0.35 (~2.857) --- + m3a = (ax >= 1.25) & (ax < (1.0 / 0.35)) + if np.any(m3a): + axm = ax[m3a] + s = 1.0 / (axm * axm) + R = _ra0 + s * ( + _ra1 + + s * (_ra2 + s * (_ra3 + s * (_ra4 + s * (_ra5 + s * (_ra6 + s * _ra7))))) + ) + S = 1.0 + s * ( + _sa1 + + s + * ( + _sa2 + + s + * (_sa3 + s * (_sa4 + s * (_sa5 + s * (_sa6 + s * (_sa7 + s * _sa8))))) + ) + ) + erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm + out[m3a] = sign[m3a] * (1.0 - erfc_val) + + # --- Region 3b: 1/0.35 <= |x| < 6 --- + m3b = (ax >= (1.0 / 0.35)) & (ax < 6.0) + if np.any(m3b): + axm = ax[m3b] + s = 1.0 / (axm * axm) + R = _rb0 + s * ( + _rb1 + s * (_rb2 + s * (_rb3 + s * (_rb4 + s * (_rb5 + s * _rb6)))) + ) + S = 1.0 + s * ( + _sb1 + + s * (_sb2 + s * (_sb3 + s * (_sb4 + s * (_sb5 + s * (_sb6 + s * _sb7))))) + ) + erfc_val = np.exp(-axm * axm - 0.5625 + R / S) / axm + out[m3b] = sign[m3b] * (1.0 - erfc_val) + + # --- Region 4: |x| >= 6 --- + m4 = ax >= 6.0 + if np.any(m4): + out[m4] = sign[m4] * 1.0 + + if scalar: + return float(out[0]) + return out + + +def _erfc(x): + """Vectorized complementary error function: erfc(x) = 1 - erf(x).""" + return 1.0 - _erf(x) diff --git a/src/flopscope/stats/_ndtri.py b/src/flopscope/stats/_ndtri.py index 9740369e6b..e7eea1a7bd 100644 --- a/src/flopscope/stats/_ndtri.py +++ b/src/flopscope/stats/_ndtri.py @@ -1,125 +1,125 @@ -"""Inverse standard normal CDF (ndtri) via a rational approximation + Newton -refinement. - -Accuracy: ~1e-12 against scipy.special.ndtri. -""" - -from __future__ import annotations - -import numpy as np - -from flopscope.stats._erf import _erf - -# --------------------------------------------------------------------------- -# Rational approximation coefficients -# --------------------------------------------------------------------------- - -_A = ( - -3.969683028665376e01, - 2.209460984245205e02, - -2.759285104469687e02, - 1.383577518672690e02, - -3.066479806614716e01, - 2.506628277459239e00, -) -_B = ( - -5.447609879822406e01, - 1.615858368580409e02, - -1.556989798598866e02, - 6.680131188771972e01, - -1.328068155288572e01, -) -_C = ( - -7.784894002430293e-03, - -3.223964580411365e-01, - -2.400758277161838e00, - -2.549732539343734e00, - 4.374664141464968e00, - 2.938163982698783e00, -) -_D = ( - 7.784695709041462e-03, - 3.224671290700398e-01, - 2.445134137142996e00, - 3.754408661907416e00, -) - -_P_LOW = 0.02425 -_P_HIGH = 1.0 - _P_LOW - -_SQRT2 = np.sqrt(2.0) -_INV_SQRT_2PI = 1.0 / np.sqrt(2.0 * np.pi) - - -def _norm_pdf_internal(x): - """Standard normal PDF (internal helper, no budget deduction).""" - return _INV_SQRT_2PI * np.exp(-0.5 * x * x) - - -def _norm_cdf_internal(x): - """Standard normal CDF (internal helper, no budget deduction).""" - return 0.5 * (1.0 + _erf(x / _SQRT2)) - - -def _ndtri(p): - """Inverse standard normal CDF. - - Maps probability *p* in (0, 1) to the quantile *x* such that - Phi(x) = p, where Phi is the standard normal CDF. - - Uses a rational approximation with one Newton-Raphson refinement step - for ~1e-12 accuracy. - - Edge cases: p=0 -> -inf, p=1 -> +inf, p<0 or p>1 -> nan. - """ - p = np.asarray(p, dtype=np.float64) - scalar = p.ndim == 0 - p = np.atleast_1d(p) - - # NaN matches neither the edge cases nor an approximation region. - out = np.full_like(p, np.nan) - - # Edge cases - out[p == 0.0] = -np.inf - out[p == 1.0] = np.inf - out[(p < 0.0) | (p > 1.0)] = np.nan - - # Lower region: 0 < p < P_LOW - m_low = (p > 0.0) & (p < _P_LOW) - if np.any(m_low): - q = np.sqrt(-2.0 * np.log(p[m_low])) - num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] - den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 - out[m_low] = num / den - - # Central region: P_LOW <= p <= P_HIGH - m_mid = (p >= _P_LOW) & (p <= _P_HIGH) - if np.any(m_mid): - q = p[m_mid] - 0.5 - r = q * q - num = ((((_A[0] * r + _A[1]) * r + _A[2]) * r + _A[3]) * r + _A[4]) * r + _A[5] - den = ((((_B[0] * r + _B[1]) * r + _B[2]) * r + _B[3]) * r + _B[4]) * r + 1.0 - out[m_mid] = q * num / den - - # Upper region: P_HIGH < p < 1 - m_high = (p > _P_HIGH) & (p < 1.0) - if np.any(m_high): - q = np.sqrt(-2.0 * np.log(1.0 - p[m_high])) - num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] - den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 - out[m_high] = -(num / den) - - # Newton refinement (one step) for the interior - m_interior = (p > 0.0) & (p < 1.0) - if np.any(m_interior): - x0 = out[m_interior] - phi = _norm_cdf_internal(x0) - pdf = _norm_pdf_internal(x0) - # Avoid division by zero in pdf - safe = pdf > 1e-300 - correction = np.where(safe, (phi - p[m_interior]) / pdf, 0.0) - out[m_interior] = x0 - correction - - if scalar: - return float(out[0]) - return out +"""Inverse standard normal CDF (ndtri) via a rational approximation + Newton +refinement. + +Accuracy: ~1e-12 against scipy.special.ndtri. +""" + +from __future__ import annotations + +import numpy as np + +from flopscope.stats._erf import _erf + +# --------------------------------------------------------------------------- +# Rational approximation coefficients +# --------------------------------------------------------------------------- + +_A = ( + -3.969683028665376e01, + 2.209460984245205e02, + -2.759285104469687e02, + 1.383577518672690e02, + -3.066479806614716e01, + 2.506628277459239e00, +) +_B = ( + -5.447609879822406e01, + 1.615858368580409e02, + -1.556989798598866e02, + 6.680131188771972e01, + -1.328068155288572e01, +) +_C = ( + -7.784894002430293e-03, + -3.223964580411365e-01, + -2.400758277161838e00, + -2.549732539343734e00, + 4.374664141464968e00, + 2.938163982698783e00, +) +_D = ( + 7.784695709041462e-03, + 3.224671290700398e-01, + 2.445134137142996e00, + 3.754408661907416e00, +) + +_P_LOW = 0.02425 +_P_HIGH = 1.0 - _P_LOW + +_SQRT2 = np.sqrt(2.0) +_INV_SQRT_2PI = 1.0 / np.sqrt(2.0 * np.pi) + + +def _norm_pdf_internal(x): + """Standard normal PDF (internal helper, no budget deduction).""" + return _INV_SQRT_2PI * np.exp(-0.5 * x * x) + + +def _norm_cdf_internal(x): + """Standard normal CDF (internal helper, no budget deduction).""" + return 0.5 * (1.0 + _erf(x / _SQRT2)) + + +def _ndtri(p): + """Inverse standard normal CDF. + + Maps probability *p* in (0, 1) to the quantile *x* such that + Phi(x) = p, where Phi is the standard normal CDF. + + Uses a rational approximation with one Newton-Raphson refinement step + for ~1e-12 accuracy. + + Edge cases: p=0 -> -inf, p=1 -> +inf, p<0 or p>1 -> nan. + """ + p = np.asarray(p, dtype=np.float64) + scalar = p.ndim == 0 + p = np.atleast_1d(p) + + # NaN matches neither the edge cases nor an approximation region. + out = np.full_like(p, np.nan) + + # Edge cases + out[p == 0.0] = -np.inf + out[p == 1.0] = np.inf + out[(p < 0.0) | (p > 1.0)] = np.nan + + # Lower region: 0 < p < P_LOW + m_low = (p > 0.0) & (p < _P_LOW) + if np.any(m_low): + q = np.sqrt(-2.0 * np.log(p[m_low])) + num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] + den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 + out[m_low] = num / den + + # Central region: P_LOW <= p <= P_HIGH + m_mid = (p >= _P_LOW) & (p <= _P_HIGH) + if np.any(m_mid): + q = p[m_mid] - 0.5 + r = q * q + num = ((((_A[0] * r + _A[1]) * r + _A[2]) * r + _A[3]) * r + _A[4]) * r + _A[5] + den = ((((_B[0] * r + _B[1]) * r + _B[2]) * r + _B[3]) * r + _B[4]) * r + 1.0 + out[m_mid] = q * num / den + + # Upper region: P_HIGH < p < 1 + m_high = (p > _P_HIGH) & (p < 1.0) + if np.any(m_high): + q = np.sqrt(-2.0 * np.log(1.0 - p[m_high])) + num = ((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5] + den = (((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1.0 + out[m_high] = -(num / den) + + # Newton refinement (one step) for the interior + m_interior = (p > 0.0) & (p < 1.0) + if np.any(m_interior): + x0 = out[m_interior] + phi = _norm_cdf_internal(x0) + pdf = _norm_pdf_internal(x0) + # Avoid division by zero in pdf + safe = pdf > 1e-300 + correction = np.where(safe, (phi - p[m_interior]) / pdf, 0.0) + out[m_interior] = x0 - correction + + if scalar: + return float(out[0]) + return out From ef72437bcfb1b3ba21ddcb5b7348092f245f3196 Mon Sep 17 00:00:00 2001 From: AminMohamed-3 <124639423+AminMohamed-3@users.noreply.github.com> Date: Wed, 9 Sep 2026 01:07:37 +0300 Subject: [PATCH 4/4] style(tests): preserve LF line endings