gamma.ts

// Adapted from npm's `gamma` to work with `decimal.js`.

import { Decimal } from "decimal.js"

const g = 7

const p = [
  0.99999999999980993, 676.5203681218851, -1259.1392167224028,
  771.32342877765313, -176.61502916214059, 12.507343278686905,
  -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
].map((value) => new Decimal(value))

const g_ln = new Decimal(607 / 128)

const p_ln = [
  0.99999999999999709182, 57.156235665862923517, -59.597960355475491248,
  14.136097974741747174, -0.49191381609762019978, 0.33994649984811888699e-4,
  0.46523628927048575665e-4, -0.98374475304879564677e-4,
  0.15808870322491248884e-3, -0.21026444172410488319e-3,
  0.2174396181152126432e-3, -0.16431810653676389022e-3,
  0.84418223983852743293e-4, -0.2619083840158140867e-4,
  0.36899182659531622704e-5,
].map((value) => new Decimal(value))

const PI = new Decimal(
  "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632789",
)

// Spouge approximation (suitable for large arguments)
function lngamma(z: Decimal) {
  if (z.lessThan(0)) return new Decimal(NaN)

  let x = p_ln[0]!

  for (let i = p_ln.length - 1; i > 0; --i) {
    x = x.plus(p_ln[i]!.dividedBy(z.plus(i)))
  }

  const t = z.plus(g_ln).plus(0.5)

  return z
    .plus(0.5)
    .times(t.ln())
    .plus(0.5 * Math.log(2 * Math.PI))
    .minus(t)
    .plus(x.ln())
    .minus(z.ln())
}

export function gamma(x: Decimal): Decimal {
  if (x.lt(0.5)) {
    return PI.dividedBy(PI.times(x).times(gamma(x.neg().plus(1))))
  } else if (x.gt(100)) {
    return lngamma(x).exp()
  } else {
    x = x.minus(1)

    let q = p[0]!

    for (let i = 1; i < g + 2; i++) {
      q = q.plus(p[i]!.dividedBy(x.plus(i)))
    }

    const t = x.plus(g + 0.5)

    return PI.times(2)
      .sqrt()
      .times(t.pow(x.plus(0.5)))
      .times(t.neg().exp())
      .times(q)
  }
}