diff --git a/Math/Functions/Functions.php b/Math/Functions/Functions.php index 0af7c4755..88c4855a3 100755 --- a/Math/Functions/Functions.php +++ b/Math/Functions/Functions.php @@ -243,6 +243,73 @@ final class Functions return \abs(self::mod($value - $start, $length)); } + private const ERF_COF = [ + -1.3026537197817094, 6.4196979235649026e-1, + 1.9476473204185836e-2,-9.561514786808631e-3,-9.46595344482036e-4, + 3.66839497852761e-4,4.2523324806907e-5,-2.0278578112534e-5, + -1.624290004647e-6,1.303655835580e-6,1.5626441722e-8,-8.5238095915e-8, + 6.529054439e-9,5.059343495e-9,-9.91364156e-10,-2.27365122e-10, + 9.6467911e-11, 2.394038e-12,-6.886027e-12,8.94487e-13, 3.13092e-13, + -1.12708e-13,3.81e-16,7.106e-15,-1.523e-15,-9.4e-17,1.21e-16,-2.8e-17 + ]; + + public static function getErf(float $x) : float + { + return $x > 0.0 + ? 1.0 - self::erfccheb($x) + : self::erfccheb(-$x) - 1.0; + } + + public static function getErfc(float $x) : float + { + return $x > 0.0 + ? self::erfccheb($x) + : 2.0 - self::erfccheb(-$x); + } + + private static function erfccheb(float $z) : float + { + $d = 0.; + $dd = 0.; + + $ncof = \count(self::ERF_COF); + + if ($z < 0.) { + throw new \InvalidArgumentException("erfccheb requires nonnegative argument"); + } + + $t = 2. / (2. + $z); + $ty = 4. * $t - 2.; + + for ($j = $ncof - 1; $j > 0; --$j) { + $tmp = $d; + $d = $ty * $d - $dd + self::ERF_COF[$j]; + $dd = $tmp; + } + + return $t * \exp(-$z * $z + 0.5*(self::ERF_COF[0] + $ty * $d) - $dd); + } + + public static function getInvErfc(float $p) : float + { + if ($p >= 2.0) { + return -100.; + } elseif ($p <= 0.0) { + return 100.; + } + + $pp = ($p < 1.0) ? $p : 2. - $p; + $t = sqrt(-2. * \log($pp / 2.)); + $x = -0.70711 * ((2.30753 + $t * 0.27061)/(1. + $t * (0.99229 + $t * 0.04481)) - $t); + + for ($j = 0; $j < 2; ++$j) { + $err = self::getErfc($x) - $pp; + $x += $err / (1.12837916709551257 * \exp(-($x * $x)) - $x * $err); + } + + return ($p < 1.0? $x : -$x); + } + /** * Calculate the value of the error function (gauss error function) * @@ -255,6 +322,7 @@ final class Functions * * @since 1.0.0 */ + /* public static function getErf(float $value) : float { if (\abs($value) > 2.2) { @@ -280,6 +348,7 @@ final class Functions return 2 / \sqrt(\M_PI) * $sum; } + */ /** * Calculate the value of the complementary error fanction @@ -293,6 +362,7 @@ final class Functions * * @since 1.0.0 */ + /* public static function getErfc(float $value) : float { if (\abs($value) <= 2.2) { @@ -323,6 +393,33 @@ final class Functions return 1 / \sqrt(\M_PI) * \exp(-$value * $value) * $q2; } + */ + + public static function getErfcInv(float $value) : float + { + if ($value >= 2) { + return \PHP_FLOAT_MIN; + } elseif ($value <= 0) { + return \PHP_FLOAT_MAX; + } elseif ($value === 1.0) { + return 0.0; + } + + if ($ge = ($value >= 1)) { + $value = 2 - $value; + } + + $t = \sqrt(-2 * \log($value / 2.0)); + $x = -0.70711 * ((2.30753 + $t * 0.27061) / (1. + $t * (0.99229 + $t * 0.04481)) - $t); + + $err = self::getErfc($x) - $value; + $x = $err / (1.12837916709551257 * \exp(-($x ** 2)) - $x * $err); + + $err = self::getErfc($x) - $value; + $x = $err / (1.12837916709551257 * \exp(-($x ** 2)) - $x * $err); + + return $ge ? -$x : $x; + } /** * Generalized hypergeometric function.