time-to-botec

erf.js (10649B)

---

 /**
 * @license Apache-2.0
 *
 * Copyright (c) 2018 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 following copyright, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/9.3.0/lib/msun/src/s_erf.c}. The implementation follows the original, but has been modified for JavaScript.
 *
 * ```text
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ```
 */
 
 'use strict';
 
 // MODULES //
 
 var isnan = require( './../../../../base/assert/is-nan' );
 var exp = require( './../../../../base/special/exp' );
 var setLowWord = require( '@stdlib/number/float64/base/set-low-word' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 var polyvalPP = require( './polyval_pp.js' );
 var polyvalQQ = require( './polyval_qq.js' );
 var polyvalPA = require( './polyval_pa.js' );
 var polyvalQA = require( './polyval_qa.js' );
 var polyvalRA = require( './polyval_ra.js' );
 var polyvalSA = require( './polyval_sa.js' );
 var polyvalRB = require( './polyval_rb.js' );
 var polyvalSB = require( './polyval_sb.js' );
 
 
 // VARIABLES //
 
 var TINY = 1.0e-300;
 var VERY_TINY = 2.848094538889218e-306; // 0x00800000, 0x00000000
 
 // 2**-28 = 1/(1<<28) = 1/268435456
 var SMALL = 3.725290298461914e-9;
 
 var ERX = 8.45062911510467529297e-1;  // 0x3FEB0AC1, 0x60000000
 
 var EFX = 1.28379167095512586316e-1;  // 0x3FC06EBA, 0x8214DB69
 var EFX8 = 1.02703333676410069053;    // 0x3FF06EBA, 0x8214DB69
 
 var PPC = 1.28379167095512558561e-1;  // 0x3FC06EBA, 0x8214DB68
 var QQC = 1.0;
 
 var PAC = -2.36211856075265944077e-3; // 0xBF6359B8, 0xBEF77538
 var QAC = 1.0;
 
 var RAC = -9.86494403484714822705e-3; // 0xBF843412, 0x600D6435
 var SAC = 1.0;
 
 var RBC = -9.86494292470009928597e-3; // 0xBF843412, 0x39E86F4A
 var SBC = 1.0;
 
 
 // MAIN //
 
 /**
 * Evaluates the error function.
 *
 * ```tex
 * \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int^{x}_{0} e^{-t^2}\ \mathrm{dt}
 * ```
 *
 * Note that
 *
 * ```tex
 * \begin{align*}
 * \operatorname{erfc}(x) &= 1 - \operatorname{erf}(x) \\
 * \operatorname{erf}(-x) &= -\operatorname{erf}(x) \\
 * \operatorname{erfc}(-x) &= 2 - \operatorname{erfc}(x)
 * \end{align*}
 * ```
 *
 * ## Method
 *
 * 1.  For \\(|x| \in [0, 0.84375)\\),
 *
 *     ```tex
 *     \operatorname{erf}(x) = x + x \cdot \operatorname{R}(x^2)
 *     ```
 *
 *     and
 *
 *     ```tex
 *     \operatorname{erfc}(x) = \begin{cases}
 *     1 - \operatorname{erf}(x) & \textrm{if}\ x \in (-.84375,0.25) \\
 *     0.5 + ((0.5-x)-x \mathrm{R}) & \textrm{if}\ x \in [0.25,0.84375)
 *     \end{cases}
 *     ```
 *
 *     where \\(R = P/Q\\) and where \\(P\\) is an odd polynomial of degree \\(8\\) and \\(Q\\) is an odd polynomial of degree \\(10\\).
 *
 *     ```tex
 *     \biggl| \mathrm{R} - \frac{\operatorname{erf}(x)-x}{x} \biggr| \leq 2^{-57.90}
 *     ```
 *
 *     <!-- <note> -->
 *
 *     The formula is derived by noting
 *
 *     ```tex
 *     \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\biggl(x - \frac{x^3}{3} + \frac{x^5}{10} - \frac{x^7}{42} + \ldots \biggr)
 *     ```
 *
 *     and that
 *
 *     ```tex
 *     \frac{2}{\sqrt{\pi}} = 1.128379167095512573896158903121545171688
 *     ```
 *
 *     is close to unity. The interval is chosen because the fix point of \\(\operatorname{erf}(x)\\) is near \\(0.6174\\) (i.e., \\(\operatorname{erf(x)} = x\\) when \\(x\\) is near \\(0.6174\\)), and, by some experiment, \\(0.84375\\) is chosen to guarantee the error is less than one ulp for \\(\operatorname{erf}(x)\\).
 *
 *     <!-- </note> -->
 *
 * 2.  For \\(|x| \in [0.84375,1.25)\\), let \\(s = |x|-1\\), and \\(c = 0.84506291151\\) rounded to single (\\(24\\) bits)
 *
 *     ```tex
 *     \operatorname{erf}(x) = \operatorname{sign}(x) \cdot \biggl(c + \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)}\biggr)
 *     ```
 *
 *     and
 *
 *     ```tex
 *     \operatorname{erfc}(x) = \begin{cases}
 *     (1-c) - \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)} & \textrm{if}\ x > 0 \\
 *     1 + \biggl(c + \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)}\biggr) & \textrm{if}\ x < 0
 *     \end{cases}
 *     ```
 *
 *     where
 *
 *     ```tex
 *     \biggl|\frac{\mathrm{P1}}{\mathrm{Q1}} - (\operatorname{erf}(|x|)-c)\biggr| \leq 2^{-59.06}
 *     ```
 *
 *     <!-- <note> -->
 *
 *     Here, we use the Taylor series expansion at \\(x = 1\\)
 *
 *     ```tex
 *     \begin{align*}
 *     \operatorname{erf}(1+s) &= \operatorname{erf}(1) + s\cdot \operatorname{poly}(s) \\
 *     &= 0.845.. + \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)}
 *     \end{align*}
 *     ```
 *
 *     using a rational approximation to approximate
 *
 *     ```tex
 *     \operatorname{erf}(1+s) - (c = (\mathrm{single})0.84506291151)
 *     ```
 *
 *     <!-- </note> -->
 *
 *     Note that, for \\(x \in [0.84375,1.25)\\), \\(|\mathrm{P1}/\mathrm{Q1}| < 0.078\\), where
 *
 *     -   \\(\operatorname{P1}(s)\\) is a degree \\(6\\) polynomial in \\(s\\)
 *     -   \\(\operatorname{Q1}(s)\\) is a degree \\(6\\) polynomial in \\(s\\)
 *
 * 3.  For \\(x \in [1.25,1/0.35)\\),
 *
 *     ```tex
 *     \begin{align*}
 *     \operatorname{erfc}(x) &= \frac{1}{x}e^{-x^2-0.5625+(\mathrm{R1}/\mathrm{S1})} \\
 *     \operatorname{erf}(x) &= 1 - \operatorname{erfc}(x)
 *     \end{align*}
 *     ```
 *
 *     where
 *
 *     -   \\(\operatorname{R1}(z)\\) is a degree \\(7\\) polynomial in \\(z\\), where \\(z = 1/x^2\\)
 *     -   \\(\operatorname{S1}(z)\\) is a degree \\(8\\) polynomial in \\(z\\)
 *
 * 4.  For \\(x \in [1/0.35,28)\\),
 *
 *     ```tex
 *     \operatorname{erfc}(x) = \begin{cases}
 *     \frac{1}{x} e^{-x^2-0.5625+(\mathrm{R2}/\mathrm{S2})} & \textrm{if}\ x > 0 \\
 *     2.0 - \frac{1}{x} e^{-x^2-0.5625+(\mathrm{R2}/\mathrm{S2})} & \textrm{if}\ -6 < x < 0 \\
 *     2.0 - \mathrm{tiny} & \textrm{if}\ x \leq -6
 *     \end{cases}
 *     ```
 *
 *     and
 *
 *     ```tex
 *     \operatorname{erf}(x) = \begin{cases}
 *     \operatorname{sign}(x) \cdot (1.0 - \operatorname{erfc}(x)) & \textrm{if}\ x < 6 \\
 *     \operatorname{sign}(x) \cdot (1.0 - \mathrm{tiny}) & \textrm{otherwise}
 *     \end{cases}
 *     ```
 *
 *     where
 *
 *     -   \\(\operatorname{R2}(z)\\) is a degree \\(6\\) polynomial in \\(z\\), where \\(z = 1/x^2\\)
 *     -   \\(\operatorname{S2}(z)\\) is a degree \\(7\\) polynomial in \\(z\\)
 *
 * 5.  For \\(x \in [28, \infty)\\),
 *
 *     ```tex
 *     \begin{align*}
 *     \operatorname{erf}(x) &= \operatorname{sign}(x) \cdot (1 - \mathrm{tiny}) & \textrm{(raise inexact)}
 *     \end{align*}
 *     ```
 *
 *     and
 *
 *     ```tex
 *     \operatorname{erfc}(x) = \begin{cases}
 *     \mathrm{tiny} \cdot \mathrm{tiny} & \textrm{if}\ x > 0\ \textrm{(raise underflow)} \\
 *     2 - \mathrm{tiny} & \textrm{if}\ x < 0
 *     \end{cases}
 *     ```
 *
 * ## Special Cases
 *
 * ```tex
 * \begin{align*}
 * \operatorname{erf}(0) &= 0 \\
 * \operatorname{erf}(-0) &= -0 \\
 * \operatorname{erf}(\infty) &= 1 \\
 * \operatorname{erf}(-\infty) &= -1 \\
 * \operatorname{erfc}(0) &= 1 \\
 * \operatorname{erfc}(\infty) &= 0 \\
 * \operatorname{erfc}(-\infty) &= 2 \\
 * \operatorname{erf}(\mathrm{NaN}) &= \mathrm{NaN} \\
 * \operatorname{erfc}(\mathrm{NaN}) &= \mathrm{NaN}
 * \end{align*}
 * ```
 *
 *
 * ## Notes
 *
 * -   To compute \\(\exp(-x^2-0.5625+(\mathrm{R}/\mathrm{S}))\\), let \\(s\\) be a single precision number and \\(s := x\\); then
 *
 *     ```tex
 *     -x^2 = -s^2 + (s-x)(s+x)
 *     ```
 *
 *     and
 *
 *     ```tex
 *     e^{-x^2-0.5626+(\mathrm{R}/\mathrm{S})} = e^{-s^2-0.5625} e^{(s-x)(s+x)+(\mathrm{R}/\mathrm{S})}
 *     ```
 *
 * -   `#4` and `#5` make use of the asymptotic series
 *
 *     ```tex
 *     \operatorname{erfc}(x) \approx \frac{e^{-x^2}}{x\sqrt{\pi}} (1 + \operatorname{poly}(1/x^2))
 *     ```
 *
 *     We use a rational approximation to approximate
 *
 *     ```tex
 *     g(s) = f(1/x^2) = \ln(\operatorname{erfc}(x) \cdot x) - x^2 + 0.5625
 *     ```
 *
 * -   The error bound for \\(\mathrm{R1}/\mathrm{S1}\\) is
 *
 *     ```tex
 *     |\mathrm{R1}/\mathrm{S1} - f(x)| < 2^{-62.57}
 *     ```
 *
 *     and for \\(\mathrm{R2}/\mathrm{S2}\\) is
 *
 *     ```tex
 *     |\mathrm{R2}/\mathrm{S2} - f(x)| < 2^{-61.52}
 *     ```
 *
 * @param {number} x - input value
 * @returns {number} function value
 *
 * @example
 * var y = erf( 2.0 );
 * // returns ~0.9953
 *
 * @example
 * var y = erf( -1.0 );
 * // returns ~-0.8427
 *
 * @example
 * var y = erf( -0.0 );
 * // returns -0.0
 *
 * @example
 * var y = erf( NaN );
 * // returns NaN
 */
 function erf( x ) {
 	var sign;
 	var ax;
 	var z;
 	var r;
 	var s;
 	var y;
 	var p;
 	var q;
 
 	// Special case: NaN
 	if ( isnan( x ) ) {
 		return NaN;
 	}
 	// Special case: +infinity
 	if ( x === PINF ) {
 		return 1.0;
 	}
 	// Special case: -infinity
 	if ( x === NINF ) {
 		return -1.0;
 	}
 	// Special case: +-0
 	if ( x === 0.0 ) {
 		return x;
 	}
 	if ( x < 0.0 ) {
 		sign = true;
 		ax = -x;
 	} else {
 		sign = false;
 		ax = x;
 	}
 	// |x| < 0.84375
 	if ( ax < 0.84375 ) {
 		if ( ax < SMALL ) {
 			if ( ax < VERY_TINY ) {
 				// Avoid underflow:
 				return 0.125 * ( (8.0*x) + (EFX8*x) );
 			}
 			return x + (EFX*x);
 		}
 		z = x * x;
 		r = PPC + ( z*polyvalPP( z ) );
 		s = QQC + ( z*polyvalQQ( z ) );
 		y = r / s;
 		return x + (x*y);
 	}
 	// 0.84375 <= |x| < 1.25
 	if ( ax < 1.25 ) {
 		s = ax - 1.0;
 		p = PAC + ( s*polyvalPA( s ) );
 		q = QAC + ( s*polyvalQA( s ) );
 		if ( sign ) {
 			return -ERX - (p/q);
 		}
 		return ERX + (p/q);
 	}
 	// +inf > |x| >= 6
 	if ( ax >= 6.0 ) {
 		if ( sign ) {
 			return TINY - 1.0; // raise inexact
 		}
 		return 1.0 - TINY; // raise inexact
 	}
 	s = 1.0 / (ax*ax);
 
 	// |x| < 1/0.35 ~ 2.857143
 	if ( ax < 2.857142857142857 ) {
 		r = RAC + ( s*polyvalRA( s ) );
 		s = SAC + ( s*polyvalSA( s ) );
 	}
 	// |x| >= 1/0.35 ~ 2.857143
 	else {
 		r = RBC + ( s*polyvalRB( s ) );
 		s = SBC + ( s*polyvalSB( s ) );
 	}
 	z = setLowWord( ax, 0 ); // pseudo-single (20-bit) precision x
 	r = exp( -(z*z) - 0.5625 ) * exp( ( (z-ax) * (z+ax) ) + (r/s) );
 	if ( sign ) {
 		return (r/ax) - 1.0;
 	}
 	return 1.0 - (r/ax);
 }
 
 
 // EXPORTS //
 
 module.exports = erf;