[Summary view]

Differences Between: [Versions 400 and 401] [Versions 401 and 402] [Versions 401 and 403]

   1  <?php
   2  
   3  namespace PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;
   4  
   5  use PhpOffice\PhpSpreadsheet\Calculation\Functions;
   6  use PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;
   7  
   8  abstract class GammaBase
   9  {
  10      private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;
  11  
  12      private const EPS = 2.22e-16;
  13  
  14      private const MAX_VALUE = 1.2e308;
  15  
  16      private const SQRT2PI = 2.5066282746310005024157652848110452530069867406099;
  17  
  18      private const MAX_ITERATIONS = 256;
  19  
  20      protected static function calculateDistribution(float $value, float $a, float $b, bool $cumulative)
  21      {
  22          if ($cumulative) {
  23              return self::incompleteGamma($a, $value / $b) / self::gammaValue($a);
  24          }
  25  
  26          return (1 / ($b ** $a * self::gammaValue($a))) * $value ** ($a - 1) * exp(0 - ($value / $b));
  27      }
  28  
  29      protected static function calculateInverse(float $probability, float $alpha, float $beta)
  30      {
  31          $xLo = 0;
  32          $xHi = $alpha * $beta * 5;
  33  
  34          $dx = 1024;
  35          $x = $xNew = 1;
  36          $i = 0;
  37  
  38          while ((abs($dx) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
  39              // Apply Newton-Raphson step
  40              $result = self::calculateDistribution($x, $alpha, $beta, true);
  41              $error = $result - $probability;
  42  
  43              if ($error == 0.0) {
  44                  $dx = 0;
  45              } elseif ($error < 0.0) {
  46                  $xLo = $x;
  47              } else {
  48                  $xHi = $x;
  49              }
  50  
  51              $pdf = self::calculateDistribution($x, $alpha, $beta, false);
  52              // Avoid division by zero
  53              if ($pdf !== 0.0) {
  54                  $dx = $error / $pdf;
  55                  $xNew = $x - $dx;
  56              }
  57  
  58              // If the NR fails to converge (which for example may be the
  59              // case if the initial guess is too rough) we apply a bisection
  60              // step to determine a more narrow interval around the root.
  61              if (($xNew < $xLo) || ($xNew > $xHi) || ($pdf == 0.0)) {
  62                  $xNew = ($xLo + $xHi) / 2;
  63                  $dx = $xNew - $x;
  64              }
  65              $x = $xNew;
  66          }
  67  
  68          if ($i === self::MAX_ITERATIONS) {
  69              return ExcelError::NA();
  70          }
  71  
  72          return $x;
  73      }
  74  
  75      //
  76      //    Implementation of the incomplete Gamma function
  77      //
  78      public static function incompleteGamma(float $a, float $x): float
  79      {
  80          static $max = 32;
  81          $summer = 0;
  82          for ($n = 0; $n <= $max; ++$n) {
  83              $divisor = $a;
  84              for ($i = 1; $i <= $n; ++$i) {
  85                  $divisor *= ($a + $i);
  86              }
  87              $summer += ($x ** $n / $divisor);
  88          }
  89  
  90          return $x ** $a * exp(0 - $x) * $summer;
  91      }
  92  
  93      //
  94      //    Implementation of the Gamma function
  95      //
  96      public static function gammaValue(float $value): float
  97      {
  98          if ($value == 0.0) {
  99              return 0;
 100          }
 101  
 102          static $p0 = 1.000000000190015;
 103          static $p = [
 104              1 => 76.18009172947146,
 105              2 => -86.50532032941677,
 106              3 => 24.01409824083091,
 107              4 => -1.231739572450155,
 108              5 => 1.208650973866179e-3,
 109              6 => -5.395239384953e-6,
 110          ];
 111  
 112          $y = $x = $value;
 113          $tmp = $x + 5.5;
 114          $tmp -= ($x + 0.5) * log($tmp);
 115  
 116          $summer = $p0;
 117          for ($j = 1; $j <= 6; ++$j) {
 118              $summer += ($p[$j] / ++$y);
 119          }
 120  
 121          return exp(0 - $tmp + log(self::SQRT2PI * $summer / $x));
 122      }
 123  
 124      /**
 125       * logGamma function.
 126       *
 127       * @version 1.1
 128       *
 129       * @author Jaco van Kooten
 130       *
 131       * Original author was Jaco van Kooten. Ported to PHP by Paul Meagher.
 132       *
 133       * The natural logarithm of the gamma function. <br />
 134       * Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz <br />
 135       * Applied Mathematics Division <br />
 136       * Argonne National Laboratory <br />
 137       * Argonne, IL 60439 <br />
 138       * <p>
 139       * References:
 140       * <ol>
 141       * <li>W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural
 142       *     Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.</li>
 143       * <li>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.</li>
 144       * <li>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.</li>
 145       * </ol>
 146       * </p>
 147       * <p>
 148       * From the original documentation:
 149       * </p>
 150       * <p>
 151       * This routine calculates the LOG(GAMMA) function for a positive real argument X.
 152       * Computation is based on an algorithm outlined in references 1 and 2.
 153       * The program uses rational functions that theoretically approximate LOG(GAMMA)
 154       * to at least 18 significant decimal digits. The approximation for X > 12 is from
 155       * reference 3, while approximations for X < 12.0 are similar to those in reference
 156       * 1, but are unpublished. The accuracy achieved depends on the arithmetic system,
 157       * the compiler, the intrinsic functions, and proper selection of the
 158       * machine-dependent constants.
 159       * </p>
 160       * <p>
 161       * Error returns: <br />
 162       * The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
 163       * The computation is believed to be free of underflow and overflow.
 164       * </p>
 165       *
 166       * @return float MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
 167       */
 168  
 169      // Log Gamma related constants
 170      private const  LG_D1 = -0.5772156649015328605195174;
 171  
 172      private const LG_D2 = 0.4227843350984671393993777;
 173  
 174      private const LG_D4 = 1.791759469228055000094023;
 175  
 176      private const LG_P1 = [
 177          4.945235359296727046734888,
 178          201.8112620856775083915565,
 179          2290.838373831346393026739,
 180          11319.67205903380828685045,
 181          28557.24635671635335736389,
 182          38484.96228443793359990269,
 183          26377.48787624195437963534,
 184          7225.813979700288197698961,
 185      ];
 186  
 187      private const LG_P2 = [
 188          4.974607845568932035012064,
 189          542.4138599891070494101986,
 190          15506.93864978364947665077,
 191          184793.2904445632425417223,
 192          1088204.76946882876749847,
 193          3338152.967987029735917223,
 194          5106661.678927352456275255,
 195          3074109.054850539556250927,
 196      ];
 197  
 198      private const LG_P4 = [
 199          14745.02166059939948905062,
 200          2426813.369486704502836312,
 201          121475557.4045093227939592,
 202          2663432449.630976949898078,
 203          29403789566.34553899906876,
 204          170266573776.5398868392998,
 205          492612579337.743088758812,
 206          560625185622.3951465078242,
 207      ];
 208  
 209      private const LG_Q1 = [
 210          67.48212550303777196073036,
 211          1113.332393857199323513008,
 212          7738.757056935398733233834,
 213          27639.87074403340708898585,
 214          54993.10206226157329794414,
 215          61611.22180066002127833352,
 216          36351.27591501940507276287,
 217          8785.536302431013170870835,
 218      ];
 219  
 220      private const LG_Q2 = [
 221          183.0328399370592604055942,
 222          7765.049321445005871323047,
 223          133190.3827966074194402448,
 224          1136705.821321969608938755,
 225          5267964.117437946917577538,
 226          13467014.54311101692290052,
 227          17827365.30353274213975932,
 228          9533095.591844353613395747,
 229      ];
 230  
 231      private const LG_Q4 = [
 232          2690.530175870899333379843,
 233          639388.5654300092398984238,
 234          41355999.30241388052042842,
 235          1120872109.61614794137657,
 236          14886137286.78813811542398,
 237          101680358627.2438228077304,
 238          341747634550.7377132798597,
 239          446315818741.9713286462081,
 240      ];
 241  
 242      private const LG_C = [
 243          -0.001910444077728,
 244          8.4171387781295e-4,
 245          -5.952379913043012e-4,
 246          7.93650793500350248e-4,
 247          -0.002777777777777681622553,
 248          0.08333333333333333331554247,
 249          0.0057083835261,
 250      ];
 251  
 252      // Rough estimate of the fourth root of logGamma_xBig
 253      private const LG_FRTBIG = 2.25e76;
 254  
 255      private const PNT68 = 0.6796875;
 256  
 257      // Function cache for logGamma
 258      private static $logGammaCacheResult = 0.0;
 259  
 260      private static $logGammaCacheX = 0.0;
 261  
 262      public static function logGamma(float $x): float
 263      {
 264          if ($x == self::$logGammaCacheX) {
 265              return self::$logGammaCacheResult;
 266          }
 267  
 268          $y = $x;
 269          if ($y > 0.0 && $y <= self::LOG_GAMMA_X_MAX_VALUE) {
 270              if ($y <= self::EPS) {
 271                  $res = -log($y);
 272              } elseif ($y <= 1.5) {
 273                  $res = self::logGamma1($y);
 274              } elseif ($y <= 4.0) {
 275                  $res = self::logGamma2($y);
 276              } elseif ($y <= 12.0) {
 277                  $res = self::logGamma3($y);
 278              } else {
 279                  $res = self::logGamma4($y);
 280              }
 281          } else {
 282              // --------------------------
 283              //    Return for bad arguments
 284              // --------------------------
 285              $res = self::MAX_VALUE;
 286          }
 287  
 288          // ------------------------------
 289          //    Final adjustments and return
 290          // ------------------------------
 291          self::$logGammaCacheX = $x;
 292          self::$logGammaCacheResult = $res;
 293  
 294          return $res;
 295      }
 296  
 297      private static function logGamma1(float $y)
 298      {
 299          // ---------------------
 300          //    EPS .LT. X .LE. 1.5
 301          // ---------------------
 302          if ($y < self::PNT68) {
 303              $corr = -log($y);
 304              $xm1 = $y;
 305          } else {
 306              $corr = 0.0;
 307              $xm1 = $y - 1.0;
 308          }
 309  
 310          $xden = 1.0;
 311          $xnum = 0.0;
 312          if ($y <= 0.5 || $y >= self::PNT68) {
 313              for ($i = 0; $i < 8; ++$i) {
 314                  $xnum = $xnum * $xm1 + self::LG_P1[$i];
 315                  $xden = $xden * $xm1 + self::LG_Q1[$i];
 316              }
 317  
 318              return $corr + $xm1 * (self::LG_D1 + $xm1 * ($xnum / $xden));
 319          }
 320  
 321          $xm2 = $y - 1.0;
 322          for ($i = 0; $i < 8; ++$i) {
 323              $xnum = $xnum * $xm2 + self::LG_P2[$i];
 324              $xden = $xden * $xm2 + self::LG_Q2[$i];
 325          }
 326  
 327          return $corr + $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
 328      }
 329  
 330      private static function logGamma2(float $y)
 331      {
 332          // ---------------------
 333          //    1.5 .LT. X .LE. 4.0
 334          // ---------------------
 335          $xm2 = $y - 2.0;
 336          $xden = 1.0;
 337          $xnum = 0.0;
 338          for ($i = 0; $i < 8; ++$i) {
 339              $xnum = $xnum * $xm2 + self::LG_P2[$i];
 340              $xden = $xden * $xm2 + self::LG_Q2[$i];
 341          }
 342  
 343          return $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
 344      }
 345  
 346      protected static function logGamma3(float $y)
 347      {
 348          // ----------------------
 349          //    4.0 .LT. X .LE. 12.0
 350          // ----------------------
 351          $xm4 = $y - 4.0;
 352          $xden = -1.0;
 353          $xnum = 0.0;
 354          for ($i = 0; $i < 8; ++$i) {
 355              $xnum = $xnum * $xm4 + self::LG_P4[$i];
 356              $xden = $xden * $xm4 + self::LG_Q4[$i];
 357          }
 358  
 359          return self::LG_D4 + $xm4 * ($xnum / $xden);
 360      }
 361  
 362      protected static function logGamma4(float $y)
 363      {
 364          // ---------------------------------
 365          //    Evaluate for argument .GE. 12.0
 366          // ---------------------------------
 367          $res = 0.0;
 368          if ($y <= self::LG_FRTBIG) {
 369              $res = self::LG_C[6];
 370              $ysq = $y * $y;
 371              for ($i = 0; $i < 6; ++$i) {
 372                  $res = $res / $ysq + self::LG_C[$i];
 373              }
 374              $res /= $y;
 375              $corr = log($y);
 376              $res = $res + log(self::SQRT2PI) - 0.5 * $corr;
 377              $res += $y * ($corr - 1.0);
 378          }
 379  
 380          return $res;
 381      }
 382  }