gamma.js 3.8 KB

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