erf.js 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152
  1. /* eslint-disable no-loss-of-precision */
  2. import { deepMap } from '../../utils/collection.js';
  3. import { sign } from '../../utils/number.js';
  4. import { factory } from '../../utils/factory.js';
  5. var name = 'erf';
  6. var dependencies = ['typed'];
  7. export var createErf = /* #__PURE__ */factory(name, dependencies, _ref => {
  8. var {
  9. typed
  10. } = _ref;
  11. /**
  12. * Compute the erf function of a value using a rational Chebyshev
  13. * approximations for different intervals of x.
  14. *
  15. * This is a translation of W. J. Cody's Fortran implementation from 1987
  16. * ( https://www.netlib.org/specfun/erf ). See the AMS publication
  17. * "Rational Chebyshev Approximations for the Error Function" by W. J. Cody
  18. * for an explanation of this process.
  19. *
  20. * For matrices, the function is evaluated element wise.
  21. *
  22. * Syntax:
  23. *
  24. * math.erf(x)
  25. *
  26. * Examples:
  27. *
  28. * math.erf(0.2) // returns 0.22270258921047847
  29. * math.erf(-0.5) // returns -0.5204998778130465
  30. * math.erf(4) // returns 0.9999999845827421
  31. *
  32. * @param {number | Array | Matrix} x A real number
  33. * @return {number | Array | Matrix} The erf of `x`
  34. */
  35. return typed('name', {
  36. number: function number(x) {
  37. var y = Math.abs(x);
  38. if (y >= MAX_NUM) {
  39. return sign(x);
  40. }
  41. if (y <= THRESH) {
  42. return sign(x) * erf1(y);
  43. }
  44. if (y <= 4.0) {
  45. return sign(x) * (1 - erfc2(y));
  46. }
  47. return sign(x) * (1 - erfc3(y));
  48. },
  49. 'Array | Matrix': typed.referToSelf(self => n => deepMap(n, self))
  50. // TODO: For complex numbers, use the approximation for the Faddeeva function
  51. // from "More Efficient Computation of the Complex Error Function" (AMS)
  52. });
  53. /**
  54. * Approximates the error function erf() for x <= 0.46875 using this function:
  55. * n
  56. * erf(x) = x * sum (p_j * x^(2j)) / (q_j * x^(2j))
  57. * j=0
  58. */
  59. function erf1(y) {
  60. var ysq = y * y;
  61. var xnum = P[0][4] * ysq;
  62. var xden = ysq;
  63. var i;
  64. for (i = 0; i < 3; i += 1) {
  65. xnum = (xnum + P[0][i]) * ysq;
  66. xden = (xden + Q[0][i]) * ysq;
  67. }
  68. return y * (xnum + P[0][3]) / (xden + Q[0][3]);
  69. }
  70. /**
  71. * Approximates the complement of the error function erfc() for
  72. * 0.46875 <= x <= 4.0 using this function:
  73. * n
  74. * erfc(x) = e^(-x^2) * sum (p_j * x^j) / (q_j * x^j)
  75. * j=0
  76. */
  77. function erfc2(y) {
  78. var xnum = P[1][8] * y;
  79. var xden = y;
  80. var i;
  81. for (i = 0; i < 7; i += 1) {
  82. xnum = (xnum + P[1][i]) * y;
  83. xden = (xden + Q[1][i]) * y;
  84. }
  85. var result = (xnum + P[1][7]) / (xden + Q[1][7]);
  86. var ysq = parseInt(y * 16) / 16;
  87. var del = (y - ysq) * (y + ysq);
  88. return Math.exp(-ysq * ysq) * Math.exp(-del) * result;
  89. }
  90. /**
  91. * Approximates the complement of the error function erfc() for x > 4.0 using
  92. * this function:
  93. *
  94. * erfc(x) = (e^(-x^2) / x) * [ 1/sqrt(pi) +
  95. * n
  96. * 1/(x^2) * sum (p_j * x^(-2j)) / (q_j * x^(-2j)) ]
  97. * j=0
  98. */
  99. function erfc3(y) {
  100. var ysq = 1 / (y * y);
  101. var xnum = P[2][5] * ysq;
  102. var xden = ysq;
  103. var i;
  104. for (i = 0; i < 4; i += 1) {
  105. xnum = (xnum + P[2][i]) * ysq;
  106. xden = (xden + Q[2][i]) * ysq;
  107. }
  108. var result = ysq * (xnum + P[2][4]) / (xden + Q[2][4]);
  109. result = (SQRPI - result) / y;
  110. ysq = parseInt(y * 16) / 16;
  111. var del = (y - ysq) * (y + ysq);
  112. return Math.exp(-ysq * ysq) * Math.exp(-del) * result;
  113. }
  114. });
  115. /**
  116. * Upper bound for the first approximation interval, 0 <= x <= THRESH
  117. * @constant
  118. */
  119. var THRESH = 0.46875;
  120. /**
  121. * Constant used by W. J. Cody's Fortran77 implementation to denote sqrt(pi)
  122. * @constant
  123. */
  124. var SQRPI = 5.6418958354775628695e-1;
  125. /**
  126. * Coefficients for each term of the numerator sum (p_j) for each approximation
  127. * interval (see W. J. Cody's paper for more details)
  128. * @constant
  129. */
  130. var P = [[3.16112374387056560e00, 1.13864154151050156e02, 3.77485237685302021e02, 3.20937758913846947e03, 1.85777706184603153e-1], [5.64188496988670089e-1, 8.88314979438837594e00, 6.61191906371416295e01, 2.98635138197400131e02, 8.81952221241769090e02, 1.71204761263407058e03, 2.05107837782607147e03, 1.23033935479799725e03, 2.15311535474403846e-8], [3.05326634961232344e-1, 3.60344899949804439e-1, 1.25781726111229246e-1, 1.60837851487422766e-2, 6.58749161529837803e-4, 1.63153871373020978e-2]];
  131. /**
  132. * Coefficients for each term of the denominator sum (q_j) for each approximation
  133. * interval (see W. J. Cody's paper for more details)
  134. * @constant
  135. */
  136. var Q = [[2.36012909523441209e01, 2.44024637934444173e02, 1.28261652607737228e03, 2.84423683343917062e03], [1.57449261107098347e01, 1.17693950891312499e02, 5.37181101862009858e02, 1.62138957456669019e03, 3.29079923573345963e03, 4.36261909014324716e03, 3.43936767414372164e03, 1.23033935480374942e03], [2.56852019228982242e00, 1.87295284992346047e00, 5.27905102951428412e-1, 6.05183413124413191e-2, 2.33520497626869185e-3]];
  137. /**
  138. * Maximum/minimum safe numbers to input to erf() (in ES6+, this number is
  139. * Number.[MAX|MIN]_SAFE_INTEGER). erf() for all numbers beyond this limit will
  140. * return 1
  141. */
  142. var MAX_NUM = Math.pow(2, 53);