simple-squiggle

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

gamma.js (4059B)

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