fix: improve numeric correctness in arithmetic operations and add regression tests

Author: arnog

CHANGELOG.md

@@ -1,5 +1,22 @@
 ### [Unreleased]
 
+- **Numeric correctness fixes** — several arithmetic operations produced
+  silently-wrong results:
+  - Complex `pow` (machine precision) used `arg ** n` instead of `arg · n` in
+    De Moivre's formula, so e.g. `i²` gave `−0.78 + 0.62i` instead of `−1`. A
+    negative integer exponent on a complex base also dropped the imaginary part
+    (`(1+i)⁻² → −0.5i` now).
+  - Complex reciprocal/`inv` (machine and arbitrary precision) divided the
+    conjugate by `|z|` instead of `|z|²` (`1/(2i) → −0.5i` now).
+  - Exact `x^(1/n)` took the n-th root of the numerator (always 1) instead of
+    the denominator, returning the base unchanged (`8^(1/3) → 2` now).
+  - Exact `floor`/`ceil`/`round` routed integers and rationals through a float,
+    losing digits beyond 2⁵³; they now compute exactly with bigints.
+  - `BigDecimal.mod` gave a wrong remainder when the quotient exceeded the
+    working precision (e.g. `1e60 mod 3`); the truncated quotient is now exact.
+  - Skewness and kurtosis used incorrect central-moment formulas; a symmetric
+    sample now has skewness 0 and `[1,2,3,4,5]` has (non-excess) kurtosis 1.7.
+
 - **Matrix (tensor) linear-algebra fixes** — operations were broken beyond the
   hardcoded 2×2 cases:
   - `Determinant` no longer throws for 4×4 and larger matrices, and now returns

src/big-decimal/big-decimal.ts

@@ -629,9 +629,19 @@ export class BigDecimal {
     if (Number.isFinite(thisExp) && Number.isFinite(otherExp)) {
       if (other.significand === 0n) return BigDecimal.NAN; // x mod 0 → NaN
       if (this.significand === 0n) return fromRaw(0n, 0); // 0 mod x → 0
-      return this.sub(this.div(other).trunc().mul(other)).toPrecision(
-        BigDecimal.precision
-      );
+
+      // Compute trunc(this / other) EXACTLY with bigints. Using the
+      // precision-bounded `div` here would round the quotient before
+      // truncating, producing a wrong remainder whenever |this / other|
+      // exceeds the working precision (e.g. `1e60 mod 3`).
+      //   this / other = (s₁ · 10^(e₁−e₂)) / s₂
+      const ediff = thisExp - otherExp;
+      const num =
+        ediff >= 0 ? this.significand * pow10(ediff) : this.significand;
+      const den =
+        ediff >= 0 ? other.significand : other.significand * pow10(-ediff);
+      const q = num / den; // bigint division truncates toward zero
+      return this.sub(fromRaw(q, 0).mul(other));
     }
 
     // Slow path: NaN or Infinity

src/compute-engine/numeric-value/big-numeric-value.ts

@@ -212,11 +212,9 @@ export class BigNumericValue extends NumericValue {
     if (this.isNegativeOne) return this;
     if (this.im === 0) return this.clone(this.decimal.inv());
 
-    const d = Math.hypot(this.re, this.im);
-    const bigD = this.decimal
-      .mul(this.decimal)
-      .add(this.im * this.im)
-      .sqrt();
+    // 1/z = conj(z) / |z|²  (not / |z|).
+    const d = this.re * this.re + this.im * this.im;
+    const bigD = this.decimal.mul(this.decimal).add(this.im * this.im);
     return this.clone({ re: this.decimal.div(bigD), im: -this.im / d });
   }
 

src/compute-engine/numeric-value/exact-numeric-value.ts

@@ -499,8 +499,10 @@ export class ExactNumericValue extends NumericValue {
         exponent = { re: exponent.re, im: exponent.im };
       } else {
         if (exponent instanceof ExactNumericValue) {
+          // Exponent 1/n (numerator rational[0] === 1) ⇒ n-th root, where n is
+          // the denominator rational[1] (not the numerator).
           if (exponent.radical === 1 && exponent.rational[0] == 1)
-            return this.root(exponent.rational[0]);
+            return this.root(Number(exponent.rational[1]));
         }
         exponent = exponent.re;
       }
@@ -709,21 +711,56 @@ export class ExactNumericValue extends NumericValue {
     return this.factory(this.bignumRe).exp();
   }
 
+  /**
+   * Floor/ceil/round of a pure rational (`radical === 1`) computed exactly with
+   * bigints. Routing through `this.re` (a float) would lose digits for
+   * integers/rationals larger than 2^53.
+   */
+  private _integerPart(mode: 'floor' | 'ceil' | 'round'): ExactNumericValue {
+    let n = BigInt(this.rational[0]);
+    let d = BigInt(this.rational[1]);
+    if (d < 0n) {
+      n = -n;
+      d = -d;
+    }
+    let q: bigint;
+    if (mode === 'round') {
+      // Round half toward +∞ (matches JS `Math.round`): floor((2n + d) / (2d)).
+      const m = 2n * n + d;
+      const dd = 2n * d;
+      q = m / dd;
+      if (m % dd !== 0n && m < 0n) q -= 1n;
+    } else {
+      q = n / d; // bigint division truncates toward zero
+      const r = n % d;
+      if (r !== 0n) {
+        if (mode === 'floor' && n < 0n) q -= 1n;
+        if (mode === 'ceil' && n > 0n) q += 1n;
+      }
+    }
+    return this.clone({ rational: [q, BigInt(1)], radical: 1 });
+  }
+
+  // An exact value is an integer iff it has no radical part and a unit
+  // denominator. (`this.type` returns `'finite_integer'`, never `'integer'`.)
   floor(): NumericValue {
     if (this.isNaN) return this.clone(NaN);
-    if (this.type === 'integer') return this;
+    if (this.radical === 1 && isInteger(this.rational)) return this;
+    if (this.radical === 1) return this._integerPart('floor');
     return this.clone(Math.floor(this.re));
   }
 
   ceil(): NumericValue {
     if (this.isNaN) return this.clone(NaN);
-    if (this.type === 'integer') return this;
+    if (this.radical === 1 && isInteger(this.rational)) return this;
+    if (this.radical === 1) return this._integerPart('ceil');
     return this.clone(Math.ceil(this.re));
   }
 
   round(): NumericValue {
     if (this.isNaN) return this.clone(NaN);
-    if (this.type === 'integer') return this;
+    if (this.radical === 1 && isInteger(this.rational)) return this;
+    if (this.radical === 1) return this._integerPart('round');
     return this.clone(Math.round(this.re));
   }
 

src/compute-engine/numeric-value/machine-numeric-value.ts

@@ -179,7 +179,8 @@ export class MachineNumericValue extends NumericValue {
     if (this.isNegativeOne) return this;
     if (this.im === 0) return this.clone(1 / this.decimal);
 
-    const d = Math.hypot(this.re, this.im);
+    // 1/z = conj(z) / |z|²  (not / |z|).
+    const d = this.re * this.re + this.im * this.im;
     return this.clone({ re: this.decimal / d, im: -this.im / d });
   }
 
@@ -373,7 +374,11 @@ export class MachineNumericValue extends NumericValue {
       if (exponent < 0) return this.clone({ im: Infinity }); // Complex/unsigned infinity
     }
 
-    if (exponent < 0) return this.clone(1 / this.decimal ** -exponent);
+    // Real base: 1/xⁿ. (Complex bases fall through to the De Moivre branch
+    // below, which handles negative exponents too — using only `this.decimal`
+    // here would drop the imaginary part.)
+    if (exponent < 0 && this.im === 0)
+      return this.clone(1 / this.decimal ** -exponent);
 
     if (this.im === 0) return this.clone(this.decimal ** exponent);
 
@@ -382,7 +387,9 @@ export class MachineNumericValue extends NumericValue {
     const modulus = Math.sqrt(a * a + b * b);
     const argument = Math.atan2(b, a);
     const newModulus = modulus ** exponent;
-    const newArgument = argument ** exponent;
+    // De Moivre: zⁿ = |z|ⁿ · (cos(n·arg) + i·sin(n·arg)). The new argument is
+    // n·arg, not argⁿ.
+    const newArgument = argument * exponent;
     return this.clone({
       re: newModulus * Math.cos(newArgument),
       im: newModulus * Math.sin(newArgument),

src/compute-engine/numerics/statistics.ts

@@ -118,57 +118,61 @@ export function bigPopulationStandardDeviation(
 export function kurtosis(values: Iterable<number>): number {
   let sum = 0;
   let sum2 = 0;
+  let sum3 = 0;
   let sum4 = 0;
   let count = 0;
   for (const op of values) {
     const v = op;
     if (!Number.isFinite(v)) return NaN;
     sum += v;
     sum2 += v * v;
+    sum3 += v * v * v;
     sum4 += v * v * v * v;
     count++;
   }
   if (count === 0) return NaN;
-  const s2 = (sum2 - (sum * sum) / count) / (count - 1);
-  return (
+  const n = count;
+  const m = sum / n;
+  // Central moments: m2 = (1/n)Σ(x−m)², m4 = (1/n)Σ(x−m)⁴.
+  const m2 = (sum2 - (sum * sum) / n) / n;
+  const m4 =
     (sum4 -
-      (4 * sum * sum2) / count +
-      (6 * sum * sum * sum) / count / count -
-      (3 * sum * sum * sum * sum) / count / count / count) /
-    (s2 * s2)
-  );
+      4 * m * sum3 +
+      6 * m * m * sum2 -
+      4 * m * m * m * sum +
+      n * m * m * m * m) /
+    n;
+  // Non-excess kurtosis β₂ = m4 / m2² (a normal distribution gives 3).
+  return m4 / (m2 * m2);
 }
 
 export function bigKurtosis(values: Iterable<BigDecimal>): BigDecimal {
   let sum = BigDecimal.ZERO;
   let sum2 = BigDecimal.ZERO;
+  let sum3 = BigDecimal.ZERO;
   let sum4 = BigDecimal.ZERO;
   let count = 0;
   for (const op of values) {
     const v = op;
     if (!v.isFinite()) return BigDecimal.NAN;
     sum = sum.add(v);
     sum2 = sum2.add(v.mul(v));
+    sum3 = sum3.add(v.mul(v).mul(v));
     sum4 = sum4.add(v.mul(v).mul(v).mul(v));
     count++;
   }
   if (count === 0) return BigDecimal.NAN;
-  const bdCount = new BigDecimal(count);
-  const s2 = sum2.sub(sum.mul(sum).div(bdCount)).div(new BigDecimal(count - 1));
-  return sum4
-    .sub(sum.mul(sum2).mul(new BigDecimal(4)).div(bdCount))
-    .add(sum.mul(sum).mul(sum).mul(new BigDecimal(6)).div(bdCount).div(bdCount))
-    .sub(
-      sum
-        .mul(sum)
-        .mul(sum)
-        .mul(sum)
-        .mul(new BigDecimal(3))
-        .div(bdCount)
-        .div(bdCount)
-        .div(bdCount)
-    )
-    .div(s2.mul(s2));
+  const m = sum.div(count); // mean
+  // Central moments: m2 = (1/n)Σ(x−m)², m4 = (1/n)Σ(x−m)⁴.
+  const m2 = sum2.sub(sum.mul(sum).div(count)).div(count);
+  const m4 = sum4
+    .sub(m.mul(sum3).mul(4))
+    .add(m.mul(m).mul(sum2).mul(6))
+    .sub(m.mul(m).mul(m).mul(sum).mul(4))
+    .add(m.mul(m).mul(m).mul(m).mul(count))
+    .div(count);
+  // Non-excess kurtosis β₂ = m4 / m2² (a normal distribution gives 3).
+  return m4.div(m2.mul(m2));
 }
 
 export function skewness(values: Iterable<number>): number {
@@ -185,9 +189,13 @@ export function skewness(values: Iterable<number>): number {
     count++;
   }
   if (count === 0) return NaN;
-  const s2 = (sum2 - (sum * sum) / count) / (count - 1);
-  const s3 = (sum3 - (sum2 * sum) / count) / (count - 1);
-  return (s3 / Math.pow(s2, 3 / 2)) * Math.sqrt(count * 1);
+  const n = count;
+  const m = sum / n;
+  // Central moments: m2 = (1/n)Σ(x−m)², m3 = (1/n)Σ(x−m)³.
+  const m2 = (sum2 - (sum * sum) / n) / n;
+  const m3 = (sum3 - 3 * m * sum2 + 3 * m * m * sum - n * m * m * m) / n;
+  // Moment coefficient of skewness g₁ = m3 / m2^(3/2).
+  return m3 / Math.pow(m2, 3 / 2);
 }
 
 export function bigSkewness(values: Iterable<BigDecimal>): BigDecimal {
@@ -204,12 +212,16 @@ export function bigSkewness(values: Iterable<BigDecimal>): BigDecimal {
     count++;
   }
   if (count === 0) return BigDecimal.NAN;
-  const bdCount = new BigDecimal(count);
-  const s2 = sum2.sub(sum.mul(sum).div(bdCount)).div(new BigDecimal(count - 1));
-  const s3 = sum3
-    .sub(sum2.mul(sum).div(bdCount))
-    .div(new BigDecimal(count - 1));
-  return s3.div(s2.pow(new BigDecimal(1.5))).mul(bdCount.sqrt());
+  const m = sum.div(count); // mean
+  // Central moments: m2 = (1/n)Σ(x−m)², m3 = (1/n)Σ(x−m)³.
+  const m2 = sum2.sub(sum.mul(sum).div(count)).div(count);
+  const m3 = sum3
+    .sub(m.mul(sum2).mul(3))
+    .add(m.mul(m).mul(sum).mul(3))
+    .sub(m.mul(m).mul(m).mul(count))
+    .div(count);
+  // Moment coefficient of skewness g₁ = m3 / m2^(3/2) = m3 / (m2·√m2).
+  return m3.div(m2.mul(m2.sqrt()));
 }
 
 export function mode(values: Iterable<number>): number {

test/big-decimal/big-decimal.test.ts

@@ -1592,6 +1592,17 @@ describe('mod()', () => {
     expect(result.sub(new BigDecimal('1.5')).abs().lt(new BigDecimal('1e-30'))).toBe(true);
   });
 
+  // REVIEW.md D5: when the quotient exceeds the working precision the old code
+  // rounded `this / other` before truncating, giving a wrong remainder. The
+  // quotient here (~3.3e59) far exceeds the precision of 50 digits.
+  test('1e60 mod 3 = 1 (quotient exceeds precision)', () => {
+    expect(new BigDecimal('1e60').mod(new BigDecimal('3')).eq(1)).toBe(true);
+  });
+
+  test('2e60 mod 3 = 2 (quotient exceeds precision)', () => {
+    expect(new BigDecimal('2e60').mod(new BigDecimal('3')).eq(2)).toBe(true);
+  });
+
   test('-10 mod 3 = -1', () => {
     const result = new BigDecimal('-10').mod(new BigDecimal('3'));
     expect(result.sub(new BigDecimal('-1')).abs().lt(new BigDecimal('1e-30'))).toBe(true);

test/compute-engine/internals/numeric-value.test.ts

@@ -1,4 +1,5 @@
 import { engine as ce } from '../../utils';
+import { ComputeEngine } from '../../../src/compute-engine';
 
 describe('constructor', () => {
   it('should create from an integer', () => {
@@ -180,3 +181,74 @@ describe('power', () => {
     expect(b.toString()).toMatchInlineSnapshot(`192`);
   });
 });
+
+// Regressions for the numeric-value bugs reported in REVIEW.md (D1–D4, D10).
+describe('Numeric-value correctness (REVIEW.md D1–D4)', () => {
+  // `BigDecimal.precision` is a global static, so a machine-precision engine is
+  // created in beforeAll (not at collection time) and the global precision is
+  // restored in afterAll — otherwise it leaks into other tests.
+  let machine: ComputeEngine;
+  let savedPrecision: number;
+  beforeAll(() => {
+    savedPrecision = ce.precision;
+    machine = new ComputeEngine({ precision: 'machine' });
+  });
+  afterAll(() => {
+    ce.precision = savedPrecision;
+  });
+
+  // D1: complex pow used De Moivre's `argument ** exponent` instead of
+  // `argument * exponent` (machine precision).
+  it('D1: machine complex pow (i^2 = -1)', () => {
+    const r = machine._numericValue(machine.complex(0, 1)).pow(2);
+    expect(r.re).toBeCloseTo(-1, 12);
+    expect(r.im).toBeCloseTo(0, 12);
+  });
+
+  // D10: negative exponent on a complex base dropped the imaginary part.
+  it('D10: machine complex negative pow ((1+i)^-2 = -0.5i)', () => {
+    const r = machine._numericValue(machine.complex(1, 1)).pow(-2);
+    expect(r.re).toBeCloseTo(0, 12);
+    expect(r.im).toBeCloseTo(-0.5, 12);
+  });
+
+  // D2: complex inv divided the conjugate by |z| instead of |z|².
+  it('D2: machine complex inv (1/(2i) = -0.5i)', () => {
+    const r = machine._numericValue(machine.complex(0, 2)).inv();
+    expect(r.re).toBeCloseTo(0, 12);
+    expect(r.im).toBeCloseTo(-0.5, 12);
+  });
+
+  it('D2: bignum complex inv (1/(2i) = -0.5i)', () => {
+    const r = ce._numericValue(ce.complex(0, 2)).inv();
+    expect(r.re).toBeCloseTo(0, 12);
+    expect(r.im).toBeCloseTo(-0.5, 12);
+  });
+
+  // D3: exact pow with a 1/n exponent took the n-th root of the numerator
+  // (always 1) instead of the denominator, so it returned the base unchanged.
+  it('D3: exact n-th root (8^(1/3) = 2, 27^(1/3) = 3)', () => {
+    expect(ce._numericValue(8).pow(ce._numericValue([1, 3])).toString()).toEqual(
+      '2'
+    );
+    expect(
+      ce._numericValue(27).pow(ce._numericValue([1, 3])).toString()
+    ).toEqual('3');
+  });
+
+  // D4: floor/ceil/round routed exact integers/rationals through a float,
+  // losing digits beyond 2^53. They now compute exactly with bigints.
+  it('D4: exact floor of a large integer keeps every digit', () => {
+    const big = 123456789012345678901234567890n;
+    expect(ce._numericValue(big).floor().toString()).toEqual(big.toString());
+  });
+
+  it('D4: exact floor/ceil/round of rationals', () => {
+    expect(ce._numericValue([7, 2]).floor().toString()).toEqual('3');
+    expect(ce._numericValue([-7, 2]).floor().toString()).toEqual('-4');
+    expect(ce._numericValue([7, 2]).ceil().toString()).toEqual('4');
+    expect(ce._numericValue([-7, 2]).ceil().toString()).toEqual('-3');
+    expect(ce._numericValue([5, 2]).round().toString()).toEqual('3');
+    expect(ce._numericValue([-5, 2]).round().toString()).toEqual('-2');
+  });
+});