mossadal/math-parser

<?php
/*
* @author      Frank Wikström <frank@mossadal.se>
* @copyright   2016 Frank Wikström
* @license     http://www.opensource.org/licenses/lgpl-license.php LGPL
*/

namespace MathParser\Extensions;

class Math {

    /**
     * Compute greates common denominator, using the Euclidean algorithm
     *
     * Compute and return gcd($a, $b)
     *
     * @param int $a
     * @param int $b
     * @return int
     */
    public static function gcd($a, $b)
    {
        $sign = 1;
        if ($a < 0) $sign = -$sign;
        if ($b < 0) $sign = -$sign;

        while ($b != 0)
        {
            $m = $a % $b;
            $a = $b;
            $b = $m;
        }
        return $sign*abs($a);
    }


    /**
     * Compute log(Gamma($a)) where $a is a positive real number
     *
     * For large values of $a ($a > 171), use Stirling asympotic expansion,
     * otherwise use the Lanczos approximation
     *
     * @param float $a
     * @throws InvalidArgumentException if $a < 0
     * @return float
     */
    public static function logGamma($a) {
        if($a < 0)
            throw new \InvalidArgumentException("Log gamma calls should be >0.");

        if ($a >= 171)    // Lanczos approximation w/ the given coefficients is accurate to 15 digits for 0 <= real(z) <= 171
            return self::logStirlingApproximation($a);
        else
            return log(self::lanczosApproximation($a));
    }

    /**
     * Compute log(Gamma($x)) using Stirling asympotic expansion
     *
     * @param float $x
     * @return float
     */
    private static function logStirlingApproximation($x) {
        $t = 0.5*log(2*pi()) - 0.5*log($x) + $x*(log($x))-$x;

        $x2 = $x * $x;
        $x3 = $x2 * $x;
        $x4 = $x3 * $x;

        $err_term = log(1 + (1.0/(12*$x)) + (1.0/(288*$x2)) - (139.0/(51840*$x3))
            - (571.0/(2488320*$x4)));

        $res = $t + $err_term;
        return $res;
    }

    /**
     * Compute factorial n! for an integer $n using iteration
     *
     * @param int $num
     * @return int
     */
    public static function Factorial($num) {
        if ($num < 0) throw new \InvalidArgumentException("Fatorial calls should be >0.");

        $rval=1;
        for ($i = 1; $i <= $num; $i++)
            $rval = $rval * $i;
        return $rval;
    }

    /**
     * Compute semi-factorial n!! for an integer $n using iteration
     *
     * @param int $num
     * @return int
     */
    public static function SemiFactorial($num) {
        if ($num < 0) throw new \InvalidArgumentException("Semifactorial calls should be >0.");

        $rval=1;
        while ($num >= 2) {
            $rval =$rval * $num;
            $num = $num-2;
        }
        return $rval;
    }

    /**
     * Compute log(Gamma($x)) using Lanczos approximation
     *
     * @param float $x
     * @return float
     */
    private static function lanczosApproximation($x) {
        $g = 7;
        $p = array(0.99999999999980993, 676.5203681218851, -1259.1392167224028,
            771.32342877765313, -176.61502916214059, 12.507343278686905,
            -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7);

        if (abs($x - floor($x)) < 1e-16)
        {
            // if we're real close to an integer, let's just compute the factorial integerly.

            if ($x >= 1)
                return self::Factorial($x - 1);
            else
                return INF;
        }
        else
        {
            $x -= 1;

            $y = $p[0];

            for ($i=1; $i < $g+2; $i++)
            {
                $y = $y + $p[$i]/($x + $i);
            }
            $t = $x + $g + 0.5;


            $res_fr = sqrt(2*pi()) * exp((($x+0.5)*log($t))-$t)*$y;

            return $res_fr;
        }
    }
}