feat: add factorial and binomial simplification rules (#289)

Add simplification rules for factorial quotients (n!/k! via partial products),
binomial detection (n!/(k!(n-k)!) → Binomial(n,k)), binomial identities
(C(n,0)→1, C(n,1)→n, C(n,n)→1), and factorial sum factoring
(n!-(n-1)! → (n-1)!*(n-1)). Works for both concrete integers and symbolic
expressions.

https://claude.ai/code/session_0114soGikCy942KAFJU8xn79

Co-authored-by: Claude <noreply@anthropic.com>

Author: arnog

src/compute-engine/boxed-expression/simplify.ts

@@ -433,6 +433,8 @@ function simplifyNonCommutativeFunction(
   const isAbsRule = because?.startsWith('|');
   // Quotient-power distribution: a/(b/c)^d -> a*(c/b)^d eliminates nested fractions
   const isQuotientPowerRule = because === 'a / (b/c)^d -> a * (c/b)^d';
+  // Factorial factoring: n! - (n-1)! -> (n-1)! * (n-1) is structurally preferred
+  const isFactorialFactoring = because === 'factor common factorial';
   // Expand may produce more nodes but enables term cancellation
   // Accept when expansion reduces terms or eliminates Power-of-Add patterns
   const isExpandWithSimplification =
@@ -474,6 +476,7 @@ function simplifyNonCommutativeFunction(
     !isRootSignRule &&
     !isAbsRule &&
     !isQuotientPowerRule &&
+    !isFactorialFactoring &&
     !isExpandWithSimplification
   )
     return steps;

src/compute-engine/symbolic/simplify-divide.ts

@@ -1,17 +1,20 @@
 import type { Expression, RuleStep } from '../global-types';
 import { isFunction, isNumber } from '../boxed-expression/type-guards';
 import { asRational } from '../boxed-expression/numerics';
+import { baseOffset } from './simplify-factorial';
 
 /**
  * Division simplification rules consolidated from simplify-rules.ts.
- * Handles ~5 patterns for simplifying Divide expressions.
  *
  * Patterns:
  * - a/a -> 1 (when a ≠ 0)
  * - 1/(1/a) -> a (when a ≠ 0)
  * - a/(1/b) -> a*b (when b ≠ 0)
  * - a/(b/c) -> a*c/b (when c ≠ 0)
  * - 0/a -> 0 (when a ≠ 0)
+ * - n!/k! -> partial product (concrete integers)
+ * - n!/k! -> (k+1)(k+2)...n (symbolic, small constant diff)
+ * - n!/(k!(n-k)!) -> Binomial(n, k)
  *
  * IMPORTANT: Do not call .simplify() on results to avoid infinite recursion.
  */
@@ -125,5 +128,119 @@ export function simplifyDivide(x: Expression): RuleStep | undefined {
     }
   }
 
+  // ── Factorial quotient: Factorial(a) / Factorial(b) ──
+  if (
+    num.operator === 'Factorial' &&
+    denom.operator === 'Factorial' &&
+    isFunction(num) &&
+    isFunction(denom)
+  ) {
+    const a = num.op1;
+    const b = denom.op1;
+
+    // Concrete integer case: compute partial product directly
+    if (
+      isNumber(a) &&
+      isNumber(b) &&
+      a.isInteger &&
+      b.isInteger &&
+      a.isNonNegative &&
+      b.isNonNegative
+    ) {
+      const aVal = BigInt(a.re);
+      const bVal = BigInt(b.re);
+      if (aVal >= bVal) {
+        // n!/k! = (k+1)(k+2)...n
+        let result = 1n;
+        for (let i = bVal + 1n; i <= aVal; i++) result *= i;
+        return { value: ce.number(result), because: 'n!/k! partial product' };
+      } else {
+        // n < k: n!/k! = 1/((n+1)(n+2)...k)
+        let result = 1n;
+        for (let i = aVal + 1n; i <= bVal; i++) result *= i;
+        return {
+          value: ce.number([1, result]),
+          because: 'n!/k! -> 1/(partial product)',
+        };
+      }
+    }
+
+    // Symbolic case: check if a and b differ by a small integer constant
+    const aBO = baseOffset(a);
+    const bBO = baseOffset(b);
+    if (aBO && bBO && aBO.base.isSame(bBO.base)) {
+      const d = aBO.offset - bBO.offset;
+      if (Number.isInteger(d) && d >= 1 && d <= 8) {
+        // a!/b! = (b+1)(b+2)...a
+        let product: Expression = b.add(ce.One);
+        for (let i = 2; i <= d; i++) {
+          product = product.mul(b.add(ce.number(i)));
+        }
+        return { value: product, because: 'n!/k! -> (k+1)..n' };
+      }
+      if (Number.isInteger(d) && d <= -1 && d >= -8) {
+        // a < b: a!/b! = 1/((a+1)(a+2)...b)
+        let product: Expression = a.add(ce.One);
+        for (let i = 2; i <= -d; i++) {
+          product = product.mul(a.add(ce.number(i)));
+        }
+        return {
+          value: ce.One.div(product),
+          because: 'n!/k! -> 1/((n+1)..k)',
+        };
+      }
+    }
+  }
+
+  // ── Binomial detection: n! / (k! * (n-k)!) → Binomial(n, k) ──
+  if (
+    num.operator === 'Factorial' &&
+    isFunction(num) &&
+    denom.operator === 'Multiply' &&
+    isFunction(denom)
+  ) {
+    const n = num.op1;
+    const factorialOps = denom.ops.filter(
+      (op) => op.operator === 'Factorial' && isFunction(op)
+    );
+    const otherOps = denom.ops.filter(
+      (op) => !(op.operator === 'Factorial' && isFunction(op))
+    );
+
+    if (
+      factorialOps.length === 2 &&
+      otherOps.length === 0 &&
+      isFunction(factorialOps[0]) &&
+      isFunction(factorialOps[1])
+    ) {
+      const k = factorialOps[0].op1;
+      const m = factorialOps[1].op1;
+
+      // Check if k + m = n (numeric)
+      if (
+        isNumber(n) &&
+        isNumber(k) &&
+        isNumber(m) &&
+        k.re + m.re === n.re
+      ) {
+        // Use the smaller of k, m for efficiency
+        const smaller = k.re <= m.re ? k : m;
+        return {
+          value: ce._fn('Binomial', [n, smaller]),
+          because: 'n!/(k!(n-k)!) -> Binomial',
+        };
+      }
+
+      // Symbolic: check if k + m structurally equals n
+      const sum = k.add(m);
+      if (sum.isSame(n)) {
+        return {
+          value: ce._fn('Binomial', [n, k]),
+          because: 'n!/(k!(n-k)!) -> Binomial',
+        };
+      }
+    }
+  }
+
   return undefined;
 }

src/compute-engine/symbolic/simplify-factorial.ts

@@ -0,0 +1,221 @@
+import type { Expression, RuleStep } from '../global-types';
+import { isFunction, isNumber } from '../boxed-expression/type-guards';
+
+/**
+ * Extracts base + integer offset from an expression.
+ * - Symbol `n` → { base: n, offset: 0 }
+ * - Add(n, 3) → { base: n, offset: 3 }
+ * - Add(n, -2) → { base: n, offset: -2 }
+ * - Multiply(2, x) → { base: Multiply(2, x), offset: 0 }
+ * - Number → null (pure numeric, no symbolic base)
+ */
+export function baseOffset(
+  expr: Expression
+): { base: Expression; offset: number } | null {
+  if (isNumber(expr)) return null;
+
+  if (expr.operator === 'Add' && isFunction(expr) && expr.nops === 2) {
+    const [op1, op2] = [expr.op1, expr.op2];
+    if (isNumber(op2) && Number.isInteger(op2.re) && !isNumber(op1))
+      return { base: op1, offset: op2.re };
+    if (isNumber(op1) && Number.isInteger(op1.re) && !isNumber(op2))
+      return { base: op2, offset: op1.re };
+  }
+
+  return { base: expr, offset: 0 };
+}
+
+/**
+ * Simplification rules for Binomial and Choose expressions.
+ *
+ * Patterns:
+ * - C(n, 0) → 1
+ * - C(n, 1) → n
+ * - C(n, n) → 1
+ * - C(n, n-1) → n
+ */
+export function simplifyBinomial(x: Expression): RuleStep | undefined {
+  if (x.operator !== 'Binomial' && x.operator !== 'Choose') return undefined;
+  if (!isFunction(x)) return undefined;
+
+  const n = x.op1;
+  const k = x.op2;
+  if (!n || !k) return undefined;
+
+  const ce = x.engine;
+
+  // C(n, 0) → 1
+  if (k.is(0)) return { value: ce.One, because: 'C(n,0) -> 1' };
+
+  // C(n, 1) → n
+  if (k.is(1)) return { value: n, because: 'C(n,1) -> n' };
+
+  // C(n, n) → 1
+  if (k.isSame(n)) return { value: ce.One, because: 'C(n,n) -> 1' };
+
+  // C(n, n-1) → n (structural check via baseOffset)
+  const nBO = baseOffset(n);
+  const kBO = baseOffset(k);
+  if (
+    nBO &&
+    kBO &&
+    nBO.base.isSame(kBO.base) &&
+    nBO.offset - kBO.offset === 1
+  )
+    return { value: n, because: 'C(n,n-1) -> n' };
+
+  return undefined;
+}
+
+/**
+ * Extracts factorial information from a term in an Add expression.
+ * Handles: Factorial(n), Negate(Factorial(n)), Multiply(c, Factorial(n))
+ */
+function extractFactorialTerm(
+  term: Expression
+): { coeff: number; factArg: Expression } | null {
+  // Direct: Factorial(n)
+  if (term.operator === 'Factorial' && isFunction(term))
+    return { coeff: 1, factArg: term.op1 };
+
+  // Negate(Factorial(n))
+  if (
+    term.operator === 'Negate' &&
+    isFunction(term) &&
+    term.op1.operator === 'Factorial' &&
+    isFunction(term.op1)
+  )
+    return { coeff: -1, factArg: term.op1.op1 };
+
+  // Multiply with integer coefficient and Factorial
+  if (term.operator === 'Multiply' && isFunction(term)) {
+    const ops = term.ops;
+    let factIdx = -1;
+    for (let i = 0; i < ops.length; i++) {
+      if (ops[i].operator === 'Factorial' && isFunction(ops[i])) {
+        if (factIdx >= 0) return null; // Multiple factorials
+        factIdx = i;
+      }
+    }
+    if (factIdx < 0) return null;
+
+    const factOp = ops[factIdx];
+    if (!isFunction(factOp)) return null;
+    const factArg = factOp.op1;
+    const others = ops.filter((_, i) => i !== factIdx);
+
+    if (
+      others.length === 1 &&
+      isNumber(others[0]) &&
+      Number.isInteger(others[0].re)
+    )
+      return { coeff: others[0].re, factArg };
+
+    return null;
+  }
+
+  return null;
+}
+
+/**
+ * Simplification for Add expressions containing factorial terms.
+ * Factors out the common (smallest) factorial for symbolic expressions.
+ *
+ * Examples:
+ * - n! - (n-1)! → (n-1)! * (n - 1)
+ * - (n+1)! - n! → n! * n
+ * - (n+1)! + n! → n! * (n + 2)
+ */
+export function simplifyFactorialAdd(x: Expression): RuleStep | undefined {
+  if (x.operator !== 'Add' || !isFunction(x)) return undefined;
+
+  const ops = x.ops;
+  if (ops.length < 2) return undefined;
+
+  // Extract factorial info from each operand
+  const factTerms: Array<{
+    coeff: number;
+    factArg: Expression;
+    index: number;
+  }> = [];
+
+  for (let i = 0; i < ops.length; i++) {
+    const info = extractFactorialTerm(ops[i]);
+    if (info) factTerms.push({ ...info, index: i });
+  }
+
+  // Need at least 2 factorial terms
+  if (factTerms.length < 2) return undefined;
+
+  const ce = x.engine;
+
+  // Get base+offset for each factorial argument
+  const bos = factTerms.map((f) => ({ ...f, bo: baseOffset(f.factArg) }));
+  const validBOs = bos.filter((b) => b.bo !== null);
+
+  if (validBOs.length < 2) return undefined;
+
+  // Check if all share the same symbolic base
+  const refBase = validBOs[0].bo!.base;
+  if (!validBOs.every((b) => b.bo!.base.isSame(refBase))) return undefined;
+
+  // Find the minimum offset (smallest factorial)
+  let minOffset = Infinity;
+  for (const b of validBOs) {
+    if (b.bo!.offset < minOffset) minOffset = b.bo!.offset;
+  }
+
+  // Find the argument expression for the minimum factorial
+  const minFactArg = validBOs.find((v) => v.bo!.offset === minOffset)!.factArg;
+
+  // Build the inner terms: for each factorial, compute the partial product
+  const innerTerms: Expression[] = [];
+
+  for (const b of validBOs) {
+    const d = b.bo!.offset - minOffset;
+
+    if (d === 0) {
+      // This term IS the minimum factorial
+      innerTerms.push(ce.number(b.coeff));
+    } else if (d > 0 && d <= 8) {
+      // Build product (minArg+1)(minArg+2)...(minArg+d)
+      let product: Expression = minFactArg.add(ce.One);
+      for (let i = 2; i <= d; i++) {
+        product = product.mul(minFactArg.add(ce.number(i)));
+      }
+      if (b.coeff === 1) {
+        innerTerms.push(product);
+      } else if (b.coeff === -1) {
+        innerTerms.push(product.neg());
+      } else {
+        innerTerms.push(ce.number(b.coeff).mul(product));
+      }
+    } else {
+      // Difference too large or negative (shouldn't happen since we found min)
+      return undefined;
+    }
+  }
+
+  // Sum the inner terms
+  const innerSum =
+    innerTerms.length === 1
+      ? innerTerms[0]
+      : ce.function('Add', innerTerms);
+
+  // Result: Factorial(minFactArg) * innerSum
+  // Use _fn to avoid canonicalization that might re-distribute the product
+  const factorialExpr = ce._fn('Factorial', [minFactArg]);
+  const factored = ce._fn('Multiply', [factorialExpr, innerSum]);
+
+  // Include any non-factorial terms from the original Add
+  const factIndices = new Set(factTerms.map((f) => f.index));
+  const nonFactTerms = ops.filter((_, i) => !factIndices.has(i));
+
+  if (nonFactTerms.length > 0)
+    return {
+      value: ce._fn('Add', [factored, ...nonFactTerms]),
+      because: 'factor common factorial',
+    };
+
+  return { value: factored, because: 'factor common factorial' };
+}