Function partialFraction

src/compute-engine/boxed-expression/factor.ts:972–1219

(
  expr: Expression,
  variable: string
)

      for (let j = 0; j < cols; j++) {
        m[row][j] = m[row][j] * pivotVal - factor * m[currentRow][j];
      }
    }

    currentRow++;
  }

  // Consistency check: after elimination, a row whose coefficient entries are
  // all zero but whose RHS is nonzero encodes 0 = c ≠ 0 — the system has no
  // solution. Without this the back-substitution below silently ignores such
  // rows and returns a spurious (often all-zero) "solution", which made
  // partialFraction collapse e.g. 1/(x·(x²+x)) to 0 (the factors x and x²+x
  // share the root 0, so the irreducible-factor template is inconsistent).
  for (let row = 0; row < rows; row++) {
    let allZero = true;
    for (let col = 0; col < numVars; col++) {
      if (m[row][col] !== 0n) {
        allZero = false;
        break;
      }
    }
    if (allZero && m[row][cols - 1] !== 0n) return null; // inconsistent
  }

  // Back substitution: extract solutions as reduced [numerator, denominator]
  const solution: [bigint, bigint][] = new Array(numVars);
  for (let col = 0; col < numVars; col++) {
    const pr = pivotRow[col];
    if (pr === -1) {
      // Free variable — set to 0
      solution[col] = [0n, 1n];
      continue;
    }

    let num = m[pr][cols - 1];
    let den = m[pr][col];
    if (den === 0n) return null; // Inconsistent

    // Normalize sign onto the numerator and reduce the fraction
    if (den < 0n) {
      num = -num;
      den = -den;
    }
    const g = gcd(num < 0n ? -num : num, den);
    solution[col] = [num / g, den / g];
  }

  return solution;
}

/** GCD of two non-negative bigints (1 for gcd(0, 0), to avoid /0). */
function gcd(a: bigint, b: bigint): bigint {
  while (b) {
    [a, b] = [b, a % b];
  }
  return a || 1n;
}