Press n or j to go to the next uncovered block, b, p or k for the previous block.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 | 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 6010x 6010x 6010x 6010x 6010x 6010x 6010x 6010x 6010x 6010x 6010x 1x 1x 6009x 6010x 1x 1x 6008x 6010x 2x 2x 6006x 6010x 1x 1x 6010x 5x 5x 6000x 6010x 3500x 3500x 6010x 2500x 2500x 2500x 6000x 6000x 6000x 6010x 3002x 3002x 3002x 3002x 2998x 6010x 1000x 1000x 1000x 1000x 1000x 1998x 1998x 1998x 6010x 1190x 1190x 1190x 1190x 808x 808x 808x 808x 808x 808x 6010x 1x 1x 1x 1x 1x | /**
* @license Apache-2.0
*
* Copyright (c) 2026 The Stdlib Authors.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*
* ## Notice
*
* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_81_0/boost/math/special_functions/detail/erf_inv.hpp}. This implementation follows the original, but has been modified for JavaScript.
*
* ```text
* (C) Copyright John Maddock 2006.
*
* Use, modification and distribution are subject to the
* Boost Software License, Version 1.0. (See accompanying file
* LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
* ```
*/
'use strict';
// MODULES //
var isnanf = require( '@stdlib/math/base/assert/is-nanf' );
var sqrtf = require( '@stdlib/math/base/special/sqrtf' );
var lnf = require( '@stdlib/math/base/special/lnf' );
var PINF = require( '@stdlib/constants/float32/pinf' );
var NINF = require( '@stdlib/constants/float32/ninf' );
var f32 = require( '@stdlib/number/float64/base/to-float32' );
var rationalFcnR1 = require( './rational_p1q1.js' );
var rationalFcnR2 = require( './rational_p2q2.js' );
var rationalFcnR3 = require( './rational_p3q3.js' );
var rationalFcnR4 = require( './rational_p4q4.js' );
var rationalFcnR5 = require( './rational_p5q5.js' );
// VARIABLES //
var Y1 = 8.91314744949340820313e-2;
var Y2 = 2.249481201171875;
var Y3 = 8.07220458984375e-1;
var Y4 = 9.3995571136474609375e-1;
var Y5 = 9.8362827301025390625e-1;
var ZERO = f32( 0.0 );
var ONE = f32( 1.0 );
var TWO = f32( 2.0 );
var NEG_ONE = f32( -1.0 );
var HALF = f32( 0.5 );
var QUARTER = f32( 0.25 );
var TEN = f32( 10.0 );
var ONE_TWO_FIVE = f32( 1.125 );
var THREE = f32( 3.0 );
var SIX = f32( 6.0 );
// MAIN //
/**
* Evaluates the inverse complementary error function (single-precision).
*
* Note that
*
* ```tex
* \operatorname{erfc^{-1}}(1-z) = \operatorname{erf^{-1}}(z)
* ```
*
* ## Method
*
* 1. For \\(|x| \leq 0.5\\), we evaluate the inverse error function using the rational approximation
*
* ```tex
* \operatorname{erf^{-1}}(x) = x(x+10)(\mathrm{Y} + \operatorname{R}(x))
* ```
*
* where \\(Y\\) is a constant and \\(\operatorname{R}(x)\\) is optimized for a low absolute error compared to \\(|Y|\\).
*
* <!-- <note> -->
*
* Max error \\(2.001849\mbox{e-}18\\). Maximum deviation found (error term at infinite precision) \\(8.030\mbox{e-}21\\).
*
* <!-- </note> -->
*
* 2. For \\(0.5 > 1-|x| \geq 0\\), we evaluate the inverse error function using the rational approximation
*
* ```tex
* \operatorname{erf^{-1}} = \frac{\sqrt{-2 \cdot \ln(1-x)}}{\mathrm{Y} + \operatorname{R}(1-x)}
* ```
*
* where \\(Y\\) is a constant, and \\(\operatorname{R}(q)\\) is optimized for a low absolute error compared to \\(Y\\).
*
* <!-- <note> -->
*
* Max error \\(7.403372\mbox{e-}17\\). Maximum deviation found (error term at infinite precision) \\(4.811\mbox{e-}20\\).
*
* <!-- </note> -->
*
* 3. For \\(1-|x| < 0.25\\), we have a series of rational approximations all of the general form
*
* ```tex
* p = \sqrt{-\ln(1-x)}
* ```
*
* Accordingly, result is given by
*
* ```tex
* \operatorname{erf^{-1}}(x) = p(\mathrm{Y} + \operatorname{R}(p-B))
* ```
*
* where \\(Y\\) is a constant, \\(B\\) is the lowest value of \\(p\\) for which the approximation is valid, and \\(\operatorname{R}(x-B)\\) is optimized for a low absolute error compared to \\(Y\\).
*
* <!-- <note> -->
*
* Almost all code will only go through the first or maybe second approximation. After that we are dealing with very small input values.
*
* - If \\(p < 3\\), max error \\(1.089051\mbox{e-}20\\).
* - If \\(p < 6\\), max error \\(8.389174\mbox{e-}21\\).
* - If \\(p < 18\\), max error \\(1.481312\mbox{e-}19\\).
* - If \\(p < 44\\), max error \\(5.697761\mbox{e-}20\\).
*
* <!-- </note> -->
*
* <!-- <note> -->
*
* The Boost library can accommodate \\(80\\) and \\(128\\) bit long doubles. JavaScript only supports a \\(64\\) bit double (IEEE 754). Accordingly, the smallest \\(p\\) (in JavaScript at the time of this writing) is \\(\sqrt{-\ln(\sim5\mbox{e-}324)} = 27.284429111150214\\).
*
* <!-- </note> -->
*
* @param {number} x - input value
* @returns {number} function value
*
* @example
* var y = erfcinvf( 0.5 );
* // returns ~0.4769
*
* @example
* var y = erfcinvf( 0.8 );
* // returns ~0.1791
*
* @example
* var y = erfcinvf( 0.0 );
* // returns Infinity
*
* @example
* var y = erfcinvf( 2.0 );
* // returns -Infinity
*
* @example
* var y = erfcinvf( NaN );
* // returns NaN
*/
function erfcinvf( x ) {
var sign;
var qs;
var q;
var g;
var r;
x = f32( x );
// Special case: NaN
if ( isnanf( x ) ) {
return NaN;
}
// Special case: 0
if ( x === ZERO ) {
return PINF;
}
// Special case: 2
if ( x === TWO ) {
return NINF;
}
// Special case: 1
if ( x === ONE ) {
return 0.0;
}
if ( x > TWO || x < ZERO ) {
return NaN;
}
// Argument reduction (reduce to interval [0,1]). If `x` is outside [0,1], we can take advantage of the complementary error function reflection formula: `erfc(-z) = 2 - erfc(z)`, by negating the result once finished.
if ( x > ONE ) {
sign = NEG_ONE;
q = f32( TWO - x );
} else {
sign = ONE;
q = x;
}
x = f32( ONE - q );
// x = 1-q <= 0.5
if ( x <= HALF ) {
g = f32( x * ( x + TEN ) );
r = rationalFcnR1( x );
return f32( sign * ( f32(g*Y1) + f32(g*r) ) );
}
// q >= 0.25
if ( q >= QUARTER ) {
g = sqrtf( f32( -2.0 * lnf( q ) ) );
q = f32( q - QUARTER );
r = rationalFcnR2( q );
return f32( sign * ( g / f32(Y2+r) ) );
}
q = sqrtf( f32( -lnf( q ) ) );
// q < 3
if ( q < THREE ) {
qs = f32( q - ONE_TWO_FIVE );
r = rationalFcnR3( qs );
return f32( sign * ( f32(Y3*q) + f32(r*q) ) );
}
// q < 6
if ( q < SIX ) {
qs = f32( q - THREE );
r = rationalFcnR4( qs );
return f32( sign * ( f32(Y4*q) + f32(r*q) ) );
}
// q < 18
qs = f32( q - SIX );
r = rationalFcnR5( qs );
return f32( sign * ( f32(Y5*q) + f32(r*q) ) );
}
// EXPORTS //
module.exports = erfcinvf;
|