phpOMS/Math/Functions/Gamma.php

<?php
/**
 * Orange Management
 *
 * PHP Version 7.4
 *
 * @package   phpOMS\Math\Functions
 * @copyright Dennis Eichhorn
 * @license   OMS License 1.0
 * @version   1.0.0
 * @link      https://orange-management.org
 */
declare(strict_types=1);

namespace phpOMS\Math\Functions;

/**
 * Gamma function
 *
 * @package phpOMS\Math\Functions
 * @license OMS License 1.0
 * @link    https://orange-management.org
 * @since   1.0.0
 */
final class Gamma
{
    /**
     * Constructor.
     *
     * @since 1.0.0
     * @codeCoverageIgnore
     */
    private function __construct()
    {
    }

    /**
     * approximation values.
     *
     * @var float[]
     * @since 1.0.0
     */
    private const LANCZOSAPPROXIMATION = [
        0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059,
        12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
    ];

    /**
     * Calculate gamma with Lanczos approximation
     *
     * @param int|float $z Value
     *
     * @return float
     *
     * @since 1.0.0
     */
    public static function lanczosApproximationReal($z) : float
    {
        if ($z < 0.5) {
            return \M_PI / (\sin(\M_PI * $z) * self::lanczosApproximationReal(1 - $z));
        }

        --$z;
        $a = self::LANCZOSAPPROXIMATION[0];
        $t = $z + 7.5;

        for ($i = 1; $i < 9; ++$i) {
            $a += self::LANCZOSAPPROXIMATION[$i] / ($z + $i);
        }

        return \sqrt(2 * \M_PI) * \pow($t, $z + 0.5) * \exp(-$t) * $a;
    }

    /**
     * Calculate gamma with Stirling approximation
     *
     * @param int|float $x Value
     *
     * @return float
     *
     * @since 1.0.0
     */
    public static function stirlingApproximation($x) : float
    {
        return \sqrt(2.0 * \M_PI / $x) * \pow($x / \M_E, $x);
    }

    /**
     * Calculate gamma with Spouge approximation
     *
     * @param int|float $z Value
     *
     * @return float
     *
     * @since 1.0.0
     */
    public static function spougeApproximation($z) : float
    {
        $k1_fact = 1.0;
        $c       = [\sqrt(2.0 * \M_PI)];

        for ($k = 1; $k < 12; ++$k) {
            $c[$k]    = \exp(12 - $k) * \pow(12 - $k, $k - 0.5) / $k1_fact;
            $k1_fact *= -$k;
        }

        $accm = $c[0];
        for ($k = 1; $k < 12; ++$k) {
            $accm += $c[$k] / ($z + $k);
        }

        $accm *= \exp(-$z - 12) * \pow($z + 12, $z + 0.5);

        return $accm / $z;
    }

    /**
     * Log of the gamma function
     *
     * @param int|float $z Value
     *
     * @return float
     *
     * @see Book: Numerical Recipes - 9780521406895
     *
     * @since 1.0.0
     */
    public static function logGamma($z) : float
    {
        static $approx = [
            76.18009172947146,-86.50532032941677,
            24.01409824083091,-1.231739572450155,
            0.1208650973866179e-2,-0.5395239384953e-5,
        ];

        $y = $z;
        $x = $z;

        $temp = $x + 5.5 - ($x + 0.5) * \log($x + 5.5);
        $sum  = 1.000000000190015;

        for ($i = 0; $i < 6; ++$i) {
            $sum += $approx[$i] / ++$y;
        }

        return -$temp + \log(\sqrt(2 * \M_PI * $sum / $x));
    }

    /**
     * Calculate gamma function value.
     *
     * Example: (7)
     *
     * @param int $k Variable
     *
     * @return int
     *
     * @since 1.0.0
     */
    public static function getGammaInteger(int $k) : int
    {
        return Functions::fact($k - 1);
    }
}