fix: interval arithmetic for the gamma function

Author: arnog

CHANGELOG.md

@@ -1,5 +1,11 @@
 ## [Unreleased]
 
+### Bug Fixes
+
+- **Interval-JS compilation for Gamma functions**
+
+### Other Changes
+
 - Updated color palettes
 
 ## 0.51.0 _2026-02-14_

src/compute-engine/interval/elementary.ts

@@ -7,7 +7,10 @@
 import type { Interval, IntervalResult } from './types';
 import { ok, containsZero, isNegative, unwrapOrPropagate } from './util';
 import { sub, mul, div } from './arithmetic';
-import { gamma as scalarGamma } from '../numerics/special-functions';
+import {
+  gamma as scalarGamma,
+  gammaln as scalarGammaln,
+} from '../numerics/special-functions';
 
 /**
  * Square root of an interval (or IntervalResult).
@@ -578,3 +581,44 @@ function _gamma(x: Interval): IntervalResult {
     hi: Math.max(scalarGamma(x.lo), scalarGamma(x.hi)),
   });
 }
+
+/**
+ * Natural logarithm of the absolute value of the gamma function.
+ *
+ * gammaln(x) = ln(|gamma(x)|)
+ *
+ * Has the same poles as gamma (at non-positive integers), but approaches
+ * -Infinity near poles instead of ±Infinity. For positive x, gammaln is
+ * monotonically increasing.
+ */
+export function gammaln(x: Interval | IntervalResult): IntervalResult {
+  const unwrapped = unwrapOrPropagate(x);
+  if (!Array.isArray(unwrapped)) return unwrapped;
+  const [xVal] = unwrapped;
+  return _gammaln(xVal);
+}
+
+function _gammaln(x: Interval): IntervalResult {
+  // Check for poles: gammaln has poles at every non-positive integer.
+  if (x.hi >= 0 && x.lo <= 0) {
+    // Interval crosses or touches zero — pole at 0
+    return { kind: 'singular', at: 0 };
+  }
+  if (x.lo < 0) {
+    // Entirely negative: check if interval spans a negative integer
+    const ceilLo = Math.ceil(x.lo);
+    const floorHi = Math.floor(x.hi);
+    // If any integer in [ceil(lo), floor(hi)] is <= 0, there's a pole
+    if (ceilLo <= floorHi) {
+      return { kind: 'singular', at: ceilLo };
+    }
+    // No pole in interval — gammaln is continuous between poles
+    const gLo = scalarGammaln(x.lo);
+    const gHi = scalarGammaln(x.hi);
+    return ok({ lo: Math.min(gLo, gHi), hi: Math.max(gLo, gHi) });
+  }
+
+  // x.lo > 0: entirely positive
+  // gammaln is monotonically increasing for x > 0
+  return ok({ lo: scalarGammaln(x.lo), hi: scalarGammaln(x.hi) });
+}

src/compute-engine/interval/index.ts

@@ -55,6 +55,8 @@ import {
   mod as _mod,
   remainder as _remainder,
   sign as _sign,
+  gamma as _gamma,
+  gammaln as _gammaln,
 } from './elementary';
 import {
   sin as _sin,
@@ -144,6 +146,8 @@ export {
   mod,
   remainder,
   sign,
+  gamma,
+  gammaln,
 } from './elementary';
 
 // Trigonometric functions
@@ -242,6 +246,8 @@ export const IntervalArithmetic = {
   mod: _mod,
   remainder: _remainder,
   sign: _sign,
+  gamma: _gamma,
+  gammaln: _gammaln,
 
   // Trigonometric
   sin: _sin,

test/compute-engine/compile-interval-js.test.ts

@@ -123,6 +123,20 @@ describe('INTERVAL JS COMPILATION - FUNCTIONS', () => {
     const fn = compile(expr, { to: 'interval-js' });
     expect(fn.code).toContain('_IA.piecewise');
   });
+
+  test('compiles Gamma', () => {
+    const expr = ce.parse('\\Gamma(x)');
+    const fn = compile(expr, { to: 'interval-js' });
+    expect(fn.success).toBe(true);
+    expect(fn.code).toContain('_IA.gamma');
+  });
+
+  test('compiles GammaLn', () => {
+    const expr = ce.box(['GammaLn', 'x']);
+    const fn = compile(expr, { to: 'interval-js' });
+    expect(fn.success).toBe(true);
+    expect(fn.code).toContain('_IA.gammaln');
+  });
 });
 
 describe('INTERVAL JS EXECUTION', () => {
@@ -421,4 +435,35 @@ describe('INTERVAL JS COMPLEX EXPRESSIONS', () => {
     expect(result.value.lo).toBeCloseTo(Math.exp(-1), 5);
     expect(result.value.hi).toBeCloseTo(1, 5);
   });
+
+  test('Gamma function positive values', () => {
+    const expr = ce.parse('\\Gamma(x)');
+    const fn = compile(expr, { to: 'interval-js' });
+
+    // Gamma(2.5) ≈ 1.329
+    const result = fn.run!({ x: { lo: 2.5, hi: 2.5 } });
+    expect(result.kind).toBe('interval');
+    expect(result.value.lo).toBeCloseTo(1.329, 2);
+    expect(result.value.hi).toBeCloseTo(1.329, 2);
+  });
+
+  test('Gamma function detects singularity at zero', () => {
+    const expr = ce.parse('\\Gamma(x)');
+    const fn = compile(expr, { to: 'interval-js' });
+
+    // Interval crossing zero should detect the pole
+    const result = fn.run!({ x: { lo: -0.5, hi: 0.5 } });
+    expect(result.kind).toBe('singular');
+  });
+
+  test('GammaLn function', () => {
+    const expr = ce.box(['GammaLn', 'x']);
+    const fn = compile(expr, { to: 'interval-js' });
+
+    // GammaLn(2.5) ≈ ln(1.329) ≈ 0.284
+    const result = fn.run!({ x: { lo: 2.5, hi: 2.5 } });
+    expect(result.kind).toBe('interval');
+    expect(result.value.lo).toBeCloseTo(0.284, 2);
+    expect(result.value.hi).toBeCloseTo(0.284, 2);
+  });
 });