Fix erf and erfc function + impl erfcinv

This commit is contained in:
Dennis Eichhorn 2023-08-14 16:37:19 +00:00
parent 1a7fe13e66
commit bc1cd635fe

View File

@ -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.