simple-squiggle

A restricted subset of Squiggle
Log | Files | Refs | README

gamma.js (3797B)

      1 import { deepMap } from '../../utils/collection.js';
      2 import { factory } from '../../utils/factory.js';
      3 import { gammaG, gammaNumber, gammaP } from '../../plain/number/index.js';
      4 var name = 'gamma';
      5 var dependencies = ['typed', 'config', 'multiplyScalar', 'pow', 'BigNumber', 'Complex'];
      6 export var createGamma = /* #__PURE__ */factory(name, dependencies, _ref => {
      7   var {
      8     typed,
      9     config,
     10     multiplyScalar,
     11     pow,
     12     BigNumber: _BigNumber,
     13     Complex: _Complex
     14   } = _ref;
     15 
     16   /**
     17    * Compute the gamma function of a value using Lanczos approximation for
     18    * small values, and an extended Stirling approximation for large values.
     19    *
     20    * For matrices, the function is evaluated element wise.
     21    *
     22    * Syntax:
     23    *
     24    *    math.gamma(n)
     25    *
     26    * Examples:
     27    *
     28    *    math.gamma(5)       // returns 24
     29    *    math.gamma(-0.5)    // returns -3.5449077018110335
     30    *    math.gamma(math.i)  // returns -0.15494982830180973 - 0.49801566811835596i
     31    *
     32    * See also:
     33    *
     34    *    combinations, factorial, permutations
     35    *
     36    * @param {number | Array | Matrix} n   A real or complex number
     37    * @return {number | Array | Matrix}    The gamma of `n`
     38    */
     39   return typed(name, {
     40     number: gammaNumber,
     41     Complex: function Complex(n) {
     42       if (n.im === 0) {
     43         return this(n.re);
     44       } // Lanczos approximation doesn't work well with real part lower than 0.5
     45       // So reflection formula is required
     46 
     47 
     48       if (n.re < 0.5) {
     49         // Euler's reflection formula
     50         // gamma(1-z) * gamma(z) = PI / sin(PI * z)
     51         // real part of Z should not be integer [sin(PI) == 0 -> 1/0 - undefined]
     52         // thanks to imperfect sin implementation sin(PI * n) != 0
     53         // we can safely use it anyway
     54         var _t = new _Complex(1 - n.re, -n.im);
     55 
     56         var r = new _Complex(Math.PI * n.re, Math.PI * n.im);
     57         return new _Complex(Math.PI).div(r.sin()).div(this(_t));
     58       } // Lanczos approximation
     59       // z -= 1
     60 
     61 
     62       n = new _Complex(n.re - 1, n.im); // x = gammaPval[0]
     63 
     64       var x = new _Complex(gammaP[0], 0); // for (i, gammaPval) in enumerate(gammaP):
     65 
     66       for (var i = 1; i < gammaP.length; ++i) {
     67         // x += gammaPval / (z + i)
     68         var gammaPval = new _Complex(gammaP[i], 0);
     69         x = x.add(gammaPval.div(n.add(i)));
     70       } // t = z + gammaG + 0.5
     71 
     72 
     73       var t = new _Complex(n.re + gammaG + 0.5, n.im); // y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x
     74 
     75       var twoPiSqrt = Math.sqrt(2 * Math.PI);
     76       var tpow = t.pow(n.add(0.5));
     77       var expt = t.neg().exp(); // y = [x] * [sqrt(2 * pi)] * [t ** (z + 0.5)] * [exp(-t)]
     78 
     79       return x.mul(twoPiSqrt).mul(tpow).mul(expt);
     80     },
     81     BigNumber: function BigNumber(n) {
     82       if (n.isInteger()) {
     83         return n.isNegative() || n.isZero() ? new _BigNumber(Infinity) : bigFactorial(n.minus(1));
     84       }
     85 
     86       if (!n.isFinite()) {
     87         return new _BigNumber(n.isNegative() ? NaN : Infinity);
     88       }
     89 
     90       throw new Error('Integer BigNumber expected');
     91     },
     92     'Array | Matrix': function ArrayMatrix(n) {
     93       return deepMap(n, this);
     94     }
     95   });
     96   /**
     97    * Calculate factorial for a BigNumber
     98    * @param {BigNumber} n
     99    * @returns {BigNumber} Returns the factorial of n
    100    */
    101 
    102   function bigFactorial(n) {
    103     if (n < 8) {
    104       return new _BigNumber([1, 1, 2, 6, 24, 120, 720, 5040][n]);
    105     }
    106 
    107     var precision = config.precision + (Math.log(n.toNumber()) | 0);
    108 
    109     var Big = _BigNumber.clone({
    110       precision: precision
    111     });
    112 
    113     if (n % 2 === 1) {
    114       return n.times(bigFactorial(new _BigNumber(n - 1)));
    115     }
    116 
    117     var p = n;
    118     var prod = new Big(n);
    119     var sum = n.toNumber();
    120 
    121     while (p > 2) {
    122       p -= 2;
    123       sum += p;
    124       prod = prod.times(sum);
    125     }
    126 
    127     return new _BigNumber(prod.toPrecision(_BigNumber.precision));
    128   }
    129 });