feat(compile) added IA functions

Author: arnog

CHANGELOG.md

@@ -17,6 +17,11 @@
   preambles using four-corner `exp(exp * ln(base))` evaluation with special cases
   for point-integer exponents and `(-1)^n`.
 
+- **`Factorial`, `Gamma`, `GammaLn` for GLSL/WGSL interval targets**: Added
+  `ia_factorial` (via `ia_gamma(x+1)`) to both GPU targets. Added `ia_gamma`
+  (Lanczos approximation) and `ia_gammaln` (Stirling asymptotic) to the WGSL
+  target, matching existing GLSL implementations.
+
 ### Bug Fixes
 
 - **Sum/Product with symbolic bounds compiled incorrectly**: Expressions like
@@ -56,12 +61,15 @@
   variables are excluded. This also fixes `Block` expressions where locally
   assigned variables (via `Assign` or `Declare`) were reported as unknowns.
 
-- **Interval `piecewise` with constant branches**: Fixed `piecewise()` returning
-  raw `Interval` objects instead of `IntervalResult` when branches evaluated to
-  constants (e.g. `_IA.point(1)`). This caused `If` expressions like
-  `\text{if}\; x \geq 0 \;\text{then}\; 1 \;\text{else}\; 0` to return
-  `undefined` kind on definite conditions and `'entire'` on indeterminate
-  conditions when compiled to `interval-js`.
+- **`Integrate` with symbolic bounds compiled incorrectly**: Same issue as
+  Sum/Product — `compileIntegrate()` used `normalizeIndexingSet()` which
+  converted symbolic bounds to `NaN`. Now uses `extractLimits()` and compiles
+  bounds as expressions.
+
+- **Interval `piecewise` test fix**: Fixed test that incorrectly accessed
+  `result.lo` directly instead of unwrapping the `IntervalResult` envelope
+  (`result.value.lo`). The `piecewise()` function correctly returns
+  `IntervalResult` objects.
 
 ## 0.51.1 _2026-02-15_
 

COMPUTE_ENGINE.md

@@ -1,10 +1,10 @@
-# Compute Engine Integration
+**NOTE** The following document has been authored by the maintainers of a
+scientific plotting library that integrates the Compute Engine for LaTeX
+expression compilation. It provides context on how CE is used within the
+plotting system and the section "Requests for CE Maintainers" outlines specific
+features and fixes needed from CE to support the plotting use case.
 
-**NOTE** This document has been authored by the maintainers of a scientific
-plotting library that integrates the Compute Engine for LaTeX expression
-compilation. It provides context on how CE is used within the plotting system
-and the section "Requests for CE Maintainers" outlines specific features and
-fixes needed from CE to support the plotting use case.
+# Compute Engine Integration
 
 How the Compute Engine (CE) is used for compiling LaTeX expressions into
 executable functions for the plotting system.
@@ -600,6 +600,31 @@ converted after adding LaTeX string support to `VectorFunction2DInput`:
 - `VectorFunction2DInput` in `types.ts` — now accepts `string` alongside
   `(x, y) => [number, number]` and `{ kind: "js", fn }`
 
+### 8. Recursive `_gpu_gamma` in `interval-glsl` preamble
+
+**Problem:** The CE's `interval-glsl` compilation target emits a monolithic
+~29KB preamble containing the full interval arithmetic library. This preamble
+includes a `_gpu_gamma(float z)` function that uses the reflection formula
+`Gamma(z) = pi / (sin(pi*z) * Gamma(1-z))` — a recursive call. GLSL forbids
+recursion, so any shader that includes this preamble fails to compile. The
+preamble is always emitted in full regardless of whether the expression actually
+uses the gamma function, so even simple expressions like `x^2 + y^2 - 1` are
+affected.
+
+**Discovery:** GPU interval arithmetic for implicit curves produced no visual
+output. Manual shader compilation in the browser console revealed the GLSL
+compiler error pointing to the recursive `_gpu_gamma` call.
+
+**Workaround:** `sanitizeIntervalPreamble()` in `shader-templates.ts` detects
+the recursive `_gpu_gamma` function via regex and replaces it with a
+non-recursive Lanczos approximation (`NON_RECURSIVE_GPU_GAMMA`) that handles
+both the `z >= 0.5` and `z < 0.5` branches inline without recursion.
+
+**Upstream fix:** The `interval-glsl` preamble should use non-recursive function
+implementations. Either replace the recursive gamma with a Lanczos/Stirling
+approximation, or emit the preamble selectively (only include functions that the
+compiled expression actually references).
+
 ## Conversion Patterns
 
 50 of 51 functions in `grid_paper.html` were converted from JS to LaTeX/CE. The
@@ -847,23 +872,29 @@ plotting integration. Ordered by impact.
 
 ### Medium Priority
 
-5. **Warn on `success: false` fallback**: When compilation fails but `run` is
+5. **Fix recursive `_gpu_gamma` in `interval-glsl` preamble (gap #8)**: The
+   preamble uses `Gamma(z) = pi / (sin(pi*z) * Gamma(1-z))` — a recursive call
+   that GLSL forbids. Replace with a non-recursive Lanczos/Stirling
+   approximation. Ideally, also emit the preamble selectively (only include
+   functions the expression actually uses) to reduce shader size.
+
+6. **Warn on `success: false` fallback**: When compilation fails but `run` is
    still set (interpreter fallback), emit a console warning. The current silent
    behavior makes it very hard to detect compilation failures.
 
-6. **Add `SphericalHarmonic(l, m, theta, phi)` and
+7. **Add `SphericalHarmonic(l, m, theta, phi)` and
    `AssociatedLegendreP(n, m, x)`**: Would allow general spherical harmonics
    without manual expansion for each (l, m) pair.
 
-7. **Support `\prod` in interval-js**: Currently only `\sum` compiles to
+8. **Support `\prod` in interval-js**: Currently only `\sum` compiles to
    interval-js. Product accumulation with `_IA.mul` would be analogous.
 
 ### Low Priority (Nice to Have)
 
-8. **GLSL compilation for Bessel, Airy, Zeta, LambertW**: These are JS-only
+9. **GLSL compilation for Bessel, Airy, Zeta, LambertW**: These are JS-only
    today. GLSL preamble-based implementations would enable GPU-accelerated
    rendering of plots using these functions.
 
-9. **Subscripted variable names in blocks**: Allow `r_1 \coloneq expr` to define
-   a variable named `r_1` rather than parsing as `Subscript(r, 1)`. This is
-   common in mathematical notation for intermediate values.
+10. **Subscripted variable names in blocks**: Allow `r_1 \coloneq expr` to
+    define a variable named `r_1` rather than parsing as `Subscript(r, 1)`. This
+    is common in mathematical notation for intermediate values.

src/compute-engine/compilation/interval-glsl-target.ts

@@ -1075,6 +1075,16 @@ IntervalResult ia_gammaln(IntervalResult x) {
   return ia_gammaln(x.value);
 }
 
+// Factorial via gamma: n! = gamma(n+1)
+IntervalResult ia_factorial(vec2 x) {
+  return ia_gamma(ia_add(ia_ok(x), ia_point(1.0)));
+}
+
+IntervalResult ia_factorial(IntervalResult x) {
+  if (ia_is_error(x.status)) return x;
+  return ia_factorial(x.value);
+}
+
 // Boolean interval comparisons
 // Returns 1.0 = true, 0.0 = false, 0.5 = maybe
 const float IA_TRUE = 1.0;
@@ -1224,6 +1234,7 @@ const INTERVAL_GLSL_FUNCTIONS: CompiledFunctions<Expression> = {
   // Special functions
   Gamma: (args, compile) => `ia_gamma(${compile(args[0])})`,
   GammaLn: (args, compile) => `ia_gammaln(${compile(args[0])})`,
+  Factorial: (args, compile) => `ia_factorial(${compile(args[0])})`,
 
   // Elementary functions
   Abs: (args, compile) => `ia_abs(${compile(args[0])})`,

src/compute-engine/compilation/interval-wgsl-target.ts

@@ -953,6 +953,101 @@ fn ia_asech(x: IntervalResult) -> IntervalResult {
   return ia_asech_v(x.value);
 }
 
+// Gamma function using Lanczos approximation (g=7, n=9 coefficients)
+// Poles at non-positive integers; minimum at x ≈ 1.4616
+fn _gpu_gamma(z_in: f32) -> f32 {
+  let PI = 3.14159265358979;
+  var z = z_in;
+  if (z < 0.5) {
+    return PI / (sin(PI * z) * _gpu_gamma(1.0 - z));
+  }
+  z -= 1.0;
+  var x = 0.99999999999980993;
+  x += 676.5203681218851 / (z + 1.0);
+  x += -1259.1392167224028 / (z + 2.0);
+  x += 771.32342877765313 / (z + 3.0);
+  x += -176.61502916214059 / (z + 4.0);
+  x += 12.507343278686905 / (z + 5.0);
+  x += -0.13857109526572012 / (z + 6.0);
+  x += 9.9843695780195716e-6 / (z + 7.0);
+  x += 1.5056327351493116e-7 / (z + 8.0);
+  let t = z + 7.5;
+  return sqrt(2.0 * PI) * pow(t, z + 0.5) * exp(-t) * x;
+}
+
+// Interval gamma function
+// Handles poles at non-positive integers and the minimum at x ≈ 1.4616
+fn ia_gamma_v(x: vec2f) -> IntervalResult {
+  let GAMMA_MIN_X = 1.4616321;
+  let GAMMA_MIN_Y = 0.8856032;
+
+  // Check for poles: interval crosses or touches zero
+  if (x.x <= 0.0 && x.y >= 0.0) {
+    return ia_singular(0.0);
+  }
+
+  // Entirely negative: check if interval spans a negative integer
+  if (x.x < 0.0) {
+    let ceilLo = ceil(x.x);
+    let floorHi = floor(x.y);
+    if (ceilLo <= floorHi) {
+      return ia_singular(ceilLo);
+    }
+    // No pole — both endpoints between same consecutive negative integers
+    let gLo = _gpu_gamma(x.x);
+    let gHi = _gpu_gamma(x.y);
+    return ia_ok(vec2f(min(gLo, gHi) - IA_EPS, max(gLo, gHi) + IA_EPS));
+  }
+
+  // Entirely positive
+  if (x.x >= GAMMA_MIN_X) {
+    // Monotonically increasing
+    return ia_ok(vec2f(_gpu_gamma(x.x) - IA_EPS, _gpu_gamma(x.y) + IA_EPS));
+  }
+  if (x.y <= GAMMA_MIN_X) {
+    // Monotonically decreasing
+    return ia_ok(vec2f(_gpu_gamma(x.y) - IA_EPS, _gpu_gamma(x.x) + IA_EPS));
+  }
+  // Crosses the minimum
+  let gMax = max(_gpu_gamma(x.x), _gpu_gamma(x.y));
+  return ia_ok(vec2f(GAMMA_MIN_Y - IA_EPS, gMax + IA_EPS));
+}
+
+fn ia_gamma(x: IntervalResult) -> IntervalResult {
+  if (ia_is_error(x.status)) { return x; }
+  return ia_gamma_v(x.value);
+}
+
+// Log-gamma using Stirling asymptotic expansion, z > 0
+fn _gpu_gammaln(z: f32) -> f32 {
+  let z3 = z * z * z;
+  return z * log(z) - z - 0.5 * log(z)
+    + 0.5 * log(2.0 * 3.14159265358979)
+    + 1.0 / (12.0 * z)
+    - 1.0 / (360.0 * z3)
+    + 1.0 / (1260.0 * z3 * z * z);
+}
+
+// Interval log-gamma — monotonically increasing for x > 0
+fn ia_gammaln_v(x: vec2f) -> IntervalResult {
+  if (x.y <= 0.0) { return ia_empty(); }
+  if (x.x > 0.0) {
+    return ia_ok(vec2f(_gpu_gammaln(x.x) - IA_EPS, _gpu_gammaln(x.y) + IA_EPS));
+  }
+  // Partial: clipped at lo
+  return ia_partial(vec2f(0.0, _gpu_gammaln(x.y) + IA_EPS), IA_PARTIAL_LO);
+}
+
+fn ia_gammaln(x: IntervalResult) -> IntervalResult {
+  if (ia_is_error(x.status)) { return x; }
+  return ia_gammaln_v(x.value);
+}
+
+// Factorial via gamma: n! = gamma(n+1)
+fn ia_factorial(x: IntervalResult) -> IntervalResult {
+  return ia_gamma(ia_add(x, ia_point(1.0)));
+}
+
 // Boolean interval comparisons
 // Returns 1.0 = true, 0.0 = false, 0.5 = maybe
 const IA_TRUE: f32 = 1.0;
@@ -1103,6 +1198,11 @@ const INTERVAL_WGSL_FUNCTIONS: CompiledFunctions<Expression> = {
   },
   Negate: (args, compile) => `ia_negate(${compile(args[0])})`,
 
+  // Special functions
+  Gamma: (args, compile) => `ia_gamma(${compile(args[0])})`,
+  GammaLn: (args, compile) => `ia_gammaln(${compile(args[0])})`,
+  Factorial: (args, compile) => `ia_factorial(${compile(args[0])})`,
+
   // Elementary functions
   Abs: (args, compile) => `ia_abs(${compile(args[0])})`,
   Ceil: (args, compile) => `ia_ceil(${compile(args[0])})`,

src/compute-engine/compilation/javascript-target.ts

@@ -63,7 +63,6 @@ import {
   variance,
 } from '../numerics/statistics';
 import { monteCarloEstimate } from '../numerics/monte-carlo';
-import { normalizeIndexingSet } from '../library/utils';
 
 import { BaseCompiler } from './base-compiler';
 import type {
@@ -1383,13 +1382,23 @@ function compileSumProduct(
  * Compile integration function
  */
 function compileIntegrate(args, _, target: CompileTarget<Expression>): string {
-  const { index, lower, upper } = normalizeIndexingSet(args[1]);
+  const { index, lowerExpr, upperExpr, lowerNum, upperNum } =
+    extractLimits(args[1]);
   const f = BaseCompiler.compile(args[0], {
     ...target,
     var: (id) => (id === index ? id : target.var(id)),
   });
 
-  return `_SYS.integrate((${index}) => (${f}), ${lower}, ${upper})`;
+  const lo =
+    lowerNum !== undefined
+      ? String(lowerNum)
+      : BaseCompiler.compile(lowerExpr, target);
+  const hi =
+    upperNum !== undefined
+      ? String(upperNum)
+      : BaseCompiler.compile(upperExpr, target);
+
+  return `_SYS.integrate((${index}) => (${f}), ${lo}, ${hi})`;
 }
 
 /**

test/compute-engine/compile-sum-product.test.ts

@@ -178,6 +178,37 @@ describe('COMPILE Sum - symbolic bounds', () => {
     expect(val.lo).toBeCloseTo(120, 10);
     expect(val.hi).toBeCloseTo(120, 10);
   });
+
+  // Symbolic LOWER bound tests
+  test('JS: sum_{k=m}^{10} k with m=3 => 3+4+...+10 = 52', () => {
+    const expr = ce.parse('\\sum_{k=m}^{10} k');
+    const result = compile(expr);
+    expect(result.success).toBe(true);
+    expect(result.run!({ m: 3 })).toBe(52);
+  });
+
+  test('JS: sum_{k=m}^{n} k with m=1 n=4 => 10', () => {
+    const expr = ce.parse('\\sum_{k=m}^{n} k');
+    const result = compile(expr);
+    expect(result.success).toBe(true);
+    expect(result.run!({ m: 1, n: 4 })).toBe(10);
+  });
+
+  test('interval-js: sum_{k=m}^{10} k with m=3 => 52', () => {
+    const expr = ce.parse('\\sum_{k=m}^{10} k');
+    const result = compile(expr, { to: 'interval-js' });
+    expect(result.success).toBe(true);
+    const val = unwrapInterval(result.run!({ m: 3 }));
+    expect(val.lo).toBeCloseTo(52, 10);
+    expect(val.hi).toBeCloseTo(52, 10);
+  });
+
+  test('JS: prod_{k=m}^{5} k with m=3 => 3*4*5 = 60', () => {
+    const expr = ce.parse('\\prod_{k=m}^{5} k');
+    const result = compile(expr);
+    expect(result.success).toBe(true);
+    expect(result.run!({ m: 3 })).toBe(60);
+  });
 });
 
 describe('COMPILE Product - interval-js', () => {
@@ -208,3 +239,22 @@ describe('COMPILE Product - interval-js', () => {
     expect(val.hi).toBe(1);
   });
 });
+
+describe('COMPILE Integrate - symbolic bounds', () => {
+  test('JS: int_0^a x dx with a=2 => 2', () => {
+    const expr = ce.parse('\\int_0^a x \\, dx');
+    const result = compile(expr);
+    expect(result.success).toBe(true);
+    // int_0^2 x dx = x^2/2 |_0^2 = 2
+    const val = result.run!({ a: 2 });
+    expect(val).toBeCloseTo(2, 1);
+  });
+
+  test('JS: int_a^b x dx with a=0 b=1 => 0.5', () => {
+    const expr = ce.parse('\\int_a^b x \\, dx');
+    const result = compile(expr);
+    expect(result.success).toBe(true);
+    const val = result.run!({ a: 0, b: 1 });
+    expect(val).toBeCloseTo(0.5, 1);
+  });
+});

test/compute-engine/compile-which.test.ts

@@ -109,16 +109,16 @@ describe('COMPILE Which', () => {
       expect(result.code).toContain('_IA.piecewise');
 
       // Test execution with point intervals
-      // piecewise returns either {lo, hi} (plain interval) or {kind, value}
+      // piecewise returns IntervalResult: {kind: 'interval', value: {lo, hi}}
       const positiveResult = result.run!({ x: 5 }) as any;
-      // When condition is definitely true, piecewise returns the true branch
-      // which is _IA.point(1) = {lo: 1, hi: 1}
-      expect(positiveResult.lo).toBe(1);
-      expect(positiveResult.hi).toBe(1);
+      const posVal = positiveResult.kind === 'interval' ? positiveResult.value : positiveResult;
+      expect(posVal.lo).toBe(1);
+      expect(posVal.hi).toBe(1);
 
       const negativeResult = result.run!({ x: -3 }) as any;
-      expect(negativeResult.lo).toBe(-1);
-      expect(negativeResult.hi).toBe(-1);
+      const negVal = negativeResult.kind === 'interval' ? negativeResult.value : negativeResult;
+      expect(negVal.lo).toBe(-1);
+      expect(negVal.hi).toBe(-1);
     });
 
     it('should compile multi-branch Which to nested piecewise', () => {