Code Coverage |
||||||||||
Lines |
Functions and Methods |
Classes and Traits |
||||||||
| Total | |
95.52% |
64 / 67 |
|
81.82% |
9 / 11 |
CRAP | |
0.00% |
0 / 1 |
| Gamma | |
95.52% |
64 / 67 |
|
81.82% |
9 / 11 |
27 | |
0.00% |
0 / 1 |
| __construct | n/a |
0 / 0 |
n/a |
0 / 0 |
1 | |||||
| gamma | |
100.00% |
1 / 1 |
|
100.00% |
1 / 1 |
1 | |||
| lanczosApproximationReal | |
100.00% |
8 / 8 |
|
100.00% |
1 / 1 |
3 | |||
| stirlingApproximation | |
100.00% |
1 / 1 |
|
100.00% |
1 / 1 |
1 | |||
| spougeApproximation | |
100.00% |
10 / 10 |
|
100.00% |
1 / 1 |
3 | |||
| logGamma | |
100.00% |
11 / 11 |
|
100.00% |
1 / 1 |
2 | |||
| getGammaInteger | |
100.00% |
1 / 1 |
|
100.00% |
1 / 1 |
1 | |||
| incompleteGammaFirst | |
100.00% |
1 / 1 |
|
100.00% |
1 / 1 |
1 | |||
| incompleteGammaSecond | |
100.00% |
1 / 1 |
|
100.00% |
1 / 1 |
1 | |||
| regularizedGamma | |
100.00% |
5 / 5 |
|
100.00% |
1 / 1 |
5 | |||
| gammaSeriesExpansion | |
90.00% |
9 / 10 |
|
0.00% |
0 / 1 |
3.01 | |||
| gammaFraction | |
88.89% |
16 / 18 |
|
0.00% |
0 / 1 |
5.03 | |||
| 1 | <?php |
| 2 | /** |
| 3 | * Jingga |
| 4 | * |
| 5 | * PHP Version 8.1 |
| 6 | * |
| 7 | * @package phpOMS\Math\Functions |
| 8 | * @copyright Dennis Eichhorn |
| 9 | * @license OMS License 2.0 |
| 10 | * @version 1.0.0 |
| 11 | * @link https://jingga.app |
| 12 | */ |
| 13 | declare(strict_types=1); |
| 14 | |
| 15 | namespace phpOMS\Math\Functions; |
| 16 | |
| 17 | /** |
| 18 | * Gamma function |
| 19 | * |
| 20 | * @package phpOMS\Math\Functions |
| 21 | * @license OMS License 2.0 |
| 22 | * @link https://jingga.app |
| 23 | * @since 1.0.0 |
| 24 | */ |
| 25 | final class Gamma |
| 26 | { |
| 27 | /** |
| 28 | * Constructor. |
| 29 | * |
| 30 | * @since 1.0.0 |
| 31 | * @codeCoverageIgnore |
| 32 | */ |
| 33 | private function __construct() |
| 34 | { |
| 35 | } |
| 36 | |
| 37 | /** |
| 38 | * Gamma function |
| 39 | * |
| 40 | * @param int|float $z Value |
| 41 | * |
| 42 | * @return float |
| 43 | * |
| 44 | * @since 1.0.0 |
| 45 | */ |
| 46 | public static function gamma(int | float $z) : float |
| 47 | { |
| 48 | return \exp(self::logGamma($z)); |
| 49 | } |
| 50 | |
| 51 | /** |
| 52 | * approximation values. |
| 53 | * |
| 54 | * @var float[] |
| 55 | * @since 1.0.0 |
| 56 | */ |
| 57 | private const LANCZOSAPPROXIMATION = [ |
| 58 | 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, |
| 59 | 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7, |
| 60 | ]; |
| 61 | |
| 62 | /** |
| 63 | * Calculate gamma with Lanczos approximation |
| 64 | * |
| 65 | * @param int|float $z Value |
| 66 | * |
| 67 | * @return float |
| 68 | * |
| 69 | * @since 1.0.0 |
| 70 | */ |
| 71 | public static function lanczosApproximationReal(int | float $z) : float |
| 72 | { |
| 73 | if ($z < 0.5) { |
| 74 | return \M_PI / (\sin(\M_PI * $z) * self::lanczosApproximationReal(1 - $z)); |
| 75 | } |
| 76 | |
| 77 | --$z; |
| 78 | $a = self::LANCZOSAPPROXIMATION[0]; |
| 79 | $t = $z + 7.5; |
| 80 | |
| 81 | for ($i = 1; $i < 9; ++$i) { |
| 82 | $a += self::LANCZOSAPPROXIMATION[$i] / ($z + $i); |
| 83 | } |
| 84 | |
| 85 | return \sqrt(2 * \M_PI) * \pow($t, $z + 0.5) * \exp(-$t) * $a; |
| 86 | } |
| 87 | |
| 88 | /** |
| 89 | * Calculate gamma with Stirling approximation |
| 90 | * |
| 91 | * @param int|float $x Value |
| 92 | * |
| 93 | * @return float |
| 94 | * |
| 95 | * @since 1.0.0 |
| 96 | */ |
| 97 | public static function stirlingApproximation(int | float $x) : float |
| 98 | { |
| 99 | return \sqrt(2.0 * \M_PI / $x) * \pow($x / \M_E, $x); |
| 100 | } |
| 101 | |
| 102 | /** |
| 103 | * Calculate gamma with Spouge approximation |
| 104 | * |
| 105 | * @param int|float $z Value |
| 106 | * |
| 107 | * @return float |
| 108 | * |
| 109 | * @since 1.0.0 |
| 110 | */ |
| 111 | public static function spougeApproximation(int | float $z) : float |
| 112 | { |
| 113 | $k1_fact = 1.0; |
| 114 | $c = [\sqrt(2.0 * \M_PI)]; |
| 115 | |
| 116 | for ($k = 1; $k < 12; ++$k) { |
| 117 | $c[$k] = \exp(12 - $k) * \pow(12 - $k, $k - 0.5) / $k1_fact; |
| 118 | $k1_fact *= -$k; |
| 119 | } |
| 120 | |
| 121 | $accm = $c[0]; |
| 122 | for ($k = 1; $k < 12; ++$k) { |
| 123 | $accm += $c[$k] / ($z + $k); |
| 124 | } |
| 125 | |
| 126 | $accm *= \exp(-$z - 12) * \pow($z + 12, $z + 0.5); |
| 127 | |
| 128 | return $accm / $z; |
| 129 | } |
| 130 | |
| 131 | /** |
| 132 | * Log of the gamma function |
| 133 | * |
| 134 | * @param int|float $z Value |
| 135 | * |
| 136 | * @return float |
| 137 | * |
| 138 | * @see Book: Numerical Recipes - 9780521406895 |
| 139 | * |
| 140 | * @since 1.0.0 |
| 141 | */ |
| 142 | public static function logGamma(int | float $z) : float |
| 143 | { |
| 144 | static $approx = [ |
| 145 | 76.18009172947146,-86.50532032941677, |
| 146 | 24.01409824083091,-1.231739572450155, |
| 147 | 0.1208650973866179e-2,-0.5395239384953e-5, |
| 148 | ]; |
| 149 | |
| 150 | $y = $z; |
| 151 | |
| 152 | $temp = $z + 5.5 - ($z + 0.5) * \log($z + 5.5); |
| 153 | $sum = 1.000000000190015; |
| 154 | |
| 155 | for ($i = 0; $i < 6; ++$i) { |
| 156 | $sum += $approx[$i] / ++$y; |
| 157 | } |
| 158 | |
| 159 | return -$temp + \log(\sqrt(2 * \M_PI) * $sum / $z); |
| 160 | } |
| 161 | |
| 162 | /** |
| 163 | * Calculate gamma function value. |
| 164 | * |
| 165 | * Example: (7) |
| 166 | * |
| 167 | * @param int $k Variable |
| 168 | * |
| 169 | * @return int |
| 170 | * |
| 171 | * @since 1.0.0 |
| 172 | */ |
| 173 | public static function getGammaInteger(int $k) : int |
| 174 | { |
| 175 | return Functions::fact($k - 1); |
| 176 | } |
| 177 | |
| 178 | /** |
| 179 | * First or lower incomplete gamma function |
| 180 | * |
| 181 | * @param float $a a |
| 182 | * @param float $x Value |
| 183 | * |
| 184 | * @return float |
| 185 | * |
| 186 | * @since 1.0.0 |
| 187 | */ |
| 188 | public static function incompleteGammaFirst(float $a, float $x) : float |
| 189 | { |
| 190 | return self::regularizedGamma($a, $x) * \exp(self::logGamma($a)); |
| 191 | } |
| 192 | |
| 193 | /** |
| 194 | * Second or upper incomplete gamma function |
| 195 | * |
| 196 | * @param float $a a |
| 197 | * @param float $x Value |
| 198 | * |
| 199 | * @return float |
| 200 | * |
| 201 | * @since 1.0.0 |
| 202 | */ |
| 203 | public static function incompleteGammaSecond(float $a, float $x) : float |
| 204 | { |
| 205 | return \exp(self::logGamma($a)) - self::regularizedGamma($a, $x) * \exp(self::logGamma($a)); |
| 206 | } |
| 207 | |
| 208 | /** |
| 209 | * Incomplete gamma function |
| 210 | * |
| 211 | * @param float $a a |
| 212 | * @param float $x Value |
| 213 | * |
| 214 | * @return float |
| 215 | * |
| 216 | * @since 1.0.0 |
| 217 | */ |
| 218 | public static function regularizedGamma(float $a, float $x) : float |
| 219 | { |
| 220 | if ($x <= 0.0 || $a <= 0.0 || $a > 10000000000.0) { |
| 221 | return 0.0; |
| 222 | } elseif ($x < $a + 1.0) { |
| 223 | return self::gammaSeriesExpansion($a, $x); |
| 224 | } |
| 225 | |
| 226 | return 1.0 - self::gammaFraction($a, $x); |
| 227 | } |
| 228 | |
| 229 | /** |
| 230 | * Gamma series expansion |
| 231 | * |
| 232 | * @param float $a a |
| 233 | * @param float $x Value |
| 234 | * |
| 235 | * @return float |
| 236 | * |
| 237 | * @see JSci |
| 238 | * @author Jaco van Kooten |
| 239 | * @license LGPL 2.1 |
| 240 | * @since 1.0.0 |
| 241 | */ |
| 242 | private static function gammaSeriesExpansion(float $a, float $x) : float |
| 243 | { |
| 244 | $ap = $a; |
| 245 | $del = 1.0 / $a; |
| 246 | $sum = $del; |
| 247 | |
| 248 | for ($i = 1; $i < 150; ++$i) { |
| 249 | ++$ap; |
| 250 | |
| 251 | $del *= $x / $ap; |
| 252 | $sum += $del; |
| 253 | |
| 254 | if ($del < $sum * 2.22e-16) { |
| 255 | return $sum * \exp(-$x + $a * \log($x) - self::logGamma($a)); |
| 256 | } |
| 257 | } |
| 258 | |
| 259 | return 0.0; |
| 260 | } |
| 261 | |
| 262 | /** |
| 263 | * Gamma fraction |
| 264 | * |
| 265 | * @param float $a a |
| 266 | * @param float $x Value |
| 267 | * |
| 268 | * @return float |
| 269 | * |
| 270 | * @see JSci |
| 271 | * @author Jaco van Kooten |
| 272 | * @license LGPL 2.1 |
| 273 | * @since 1.0.0 |
| 274 | */ |
| 275 | private static function gammaFraction(float $a, float $x) : float |
| 276 | { |
| 277 | $b = $x + 1.0 - $a; |
| 278 | $c = 1.0 / 1.18e-37; |
| 279 | $d = 1.0 / $b; |
| 280 | $h = $d; |
| 281 | $del = 0.0; |
| 282 | |
| 283 | for ($i = 1; $i < 150 && \abs($del - 1.0) > 2.22e-16; ++$i) { |
| 284 | $an = - $i * ($i - $a); |
| 285 | $b += 2.0; |
| 286 | $d = $an * $d + $b; |
| 287 | $c = $b + $an / $c; |
| 288 | |
| 289 | if (\abs($c) < 1.18e-37) { |
| 290 | $c = 1.18e-37; |
| 291 | } |
| 292 | |
| 293 | if (\abs($d) < 1.18e-37) { |
| 294 | $d = 1.18e-37; |
| 295 | } |
| 296 | |
| 297 | $d = 1.0 / $d; |
| 298 | $del = $d * $c; |
| 299 | $h *= $del; |
| 300 | } |
| 301 | |
| 302 | return \exp(-$x + $a * \log($x) - self::logGamma($a)) * $h; |
| 303 | } |
| 304 | } |