Implement a faster fraction to monzo + residual converter

.eslintrc.json

@@ -6,6 +6,7 @@
     "@typescript-eslint/no-explicit-any": 0,
     "@typescript-eslint/ban-ts-comment": 0,
     "node/no-unpublished-import": ["error", {"allowModules": ["vitest"]}],
+    "prefer-const": ["error", {"destructuring": "all"}],
     "no-restricted-syntax": ["error", "SequenceExpression"]
   }
 }

src/__benchmarks__/monzo.bench.ts

@@ -1,13 +1,21 @@
 import {describe, bench} from 'vitest';
 
 // It's important to use the distributed versions for a realistic comparison
-import {toMonzoLegacy, primeLimitLegacy} from '../../legacy/legacy';
-import {toMonzo, primeLimit} from '../../dist';
+import {
+  toMonzoLegacy,
+  primeLimitLegacy,
+  toMonzoAndResidualLegacy,
+} from '../../legacy/legacy';
+import {toMonzo, primeLimit, toMonzoAndResidual} from '../../dist';
 
 function randInt() {
   return Math.ceil(Math.random() * 1000000000);
 }
 
+function randNumComponents() {
+  return 2 + Math.floor(Math.random() * 10);
+}
+
 describe('Number to prime exponent vector conversion', () => {
   bench('Old implementation', () => {
     try {
@@ -31,3 +39,23 @@ describe('Prime limit calculator', () => {
     primeLimit(randInt());
   });
 });
+
+describe('Monzo with residual', () => {
+  bench('Current implementation', () => {
+    toMonzoAndResidual(randInt(), randNumComponents());
+  });
+
+  bench('Old implementation', () => {
+    toMonzoAndResidualLegacy(randInt(), randNumComponents());
+  });
+});
+
+describe('Monzo with residual (bigint)', () => {
+  bench('Current implementation', () => {
+    toMonzoAndResidual(BigInt(randInt()), randNumComponents());
+  });
+
+  bench('Old implementation', () => {
+    toMonzoAndResidualLegacy(BigInt(randInt()), randNumComponents());
+  });
+});

src/__tests__/monzo.spec.ts

@@ -8,6 +8,34 @@ import {
   toMonzoAndResidual,
 } from '../monzo';
 
+function toMonzoAndResidual11(n: number): [number[], number] {
+  const result = [0, 0, 0, 0, 0];
+  if (!n) {
+    return [result, n];
+  }
+  while (n % 2 === 0) {
+    n /= 2;
+    result[0]++;
+  }
+  while (n % 3 === 0) {
+    n /= 3;
+    result[1]++;
+  }
+  while (n % 5 === 0) {
+    n /= 5;
+    result[2]++;
+  }
+  while (n % 7 === 0) {
+    n /= 7;
+    result[3]++;
+  }
+  while (n % 11 === 0) {
+    n /= 11;
+    result[4]++;
+  }
+  return [result, n];
+}
+
 describe('Monzo converter', () => {
   it('can break down an integer to its prime components', () => {
     const monzo = toMonzo(360);
@@ -145,6 +173,62 @@ describe('Fraction to monzo converter', () => {
     expect(monzo[0]).toBe(1);
     expect(monzo[1]).toBe(0);
   });
+
+  it('leaves negative residual for integers', () => {
+    const [monzo, residual] = toMonzoAndResidual(-10, 2);
+    expect(residual.toFraction()).toBe('-5');
+    expect(monzo).toHaveLength(2);
+    expect(monzo[0]).toBe(1);
+    expect(monzo[1]).toBe(0);
+  });
+
+  it('works just below the int32 boundary', () => {
+    const [monzo, residual] = toMonzoAndResidual(2 ** 30, 1);
+    expect(monzo).toEqual([30]);
+    expect(residual.isUnity()).toBe(true);
+  });
+
+  it('works at the int32 boundary', () => {
+    const [monzo, residual] = toMonzoAndResidual(2 ** 31, 1);
+    expect(monzo).toEqual([31]);
+    expect(residual.isUnity()).toBe(true);
+  });
+
+  it('works just above the int32 boundary', () => {
+    const [monzo, residual] = toMonzoAndResidual(2 ** 32, 1);
+    expect(monzo).toEqual([32]);
+    expect(residual.isUnity()).toBe(true);
+  });
+
+  it('works just below the IEEE limit', () => {
+    const [monzo, residual] = toMonzoAndResidual(2n ** 1023n, 1);
+    expect(monzo).toEqual([1023]);
+    expect(residual).toBe(1n);
+  });
+
+  it('works at the IEEE limit', () => {
+    const [monzo, residual] = toMonzoAndResidual(2n ** 1024n, 1);
+    expect(monzo).toEqual([1024]);
+    expect(residual).toBe(1n);
+  });
+
+  it('works just above the IEEE limit', () => {
+    const [monzo, residual] = toMonzoAndResidual(2n ** 1025n, 1);
+    expect(monzo).toEqual([1025]);
+    expect(residual).toBe(1n);
+  });
+
+  it('agrees with the reference implementation', () => {
+    for (let n = -1000; n <= 1000; ++n) {
+      const [monzo, residual] = toMonzoAndResidual(n, 5);
+      const [bigMonzo, bigResidual] = toMonzoAndResidual(BigInt(n), 5);
+      const [reference, refResidual] = toMonzoAndResidual11(n);
+      expect(monzo).toEqual(reference);
+      expect(bigMonzo).toEqual(reference);
+      expect(residual.equals(refResidual)).toBe(true);
+      expect(bigResidual).toBe(BigInt(refResidual));
+    }
+  });
 });
 
 describe('Monzo to fraction converter', () => {

src/legacy.ts

@@ -162,3 +162,86 @@ function bigIntPrimeLimit(
     }
   }
 }
+
+/**
+ * Extract the exponents of the prime factors of a rational number.
+ * @param n Rational number to convert to a monzo.
+ * @param numberOfComponents Number of components in the result.
+ * @returns The monzo representing `n` and a multiplicative residue that cannot be represented in the given limit.
+ */
+export function toMonzoAndResidualLegacy(
+  n: bigint,
+  numberOfComponents: number
+): [Monzo, bigint];
+export function toMonzoAndResidualLegacy(
+  n: FractionValue,
+  numberOfComponents: number
+): [Monzo, Fraction];
+export function toMonzoAndResidualLegacy(
+  n: FractionValue | bigint,
+  numberOfComponents: number
+): [Monzo, Fraction] | [Monzo, bigint] {
+  if (typeof n === 'bigint') {
+    return bigIntToMonzoAndResidualLegacy(n, numberOfComponents);
+  }
+  n = new Fraction(n);
+  const numerator = n.n;
+  const denominator = n.d;
+
+  if (!n.n) {
+    return [Array(numberOfComponents).fill(0), new Fraction(0)];
+  }
+
+  let nProbe = 1;
+  let dProbe = 1;
+
+  const result = Array(numberOfComponents).fill(-1);
+  for (let i = 0; i < numberOfComponents; ++i) {
+    let lastProbe;
+    do {
+      lastProbe = nProbe;
+      nProbe *= PRIMES[i];
+      result[i]++;
+    } while (numerator % nProbe === 0);
+    nProbe = lastProbe;
+
+    // The fraction is in lowest terms so we know that positive components exclude negative components.
+    if (result[i]) {
+      continue;
+    }
+
+    result[i] = 1;
+    do {
+      lastProbe = dProbe;
+      dProbe *= PRIMES[i];
+      result[i]--;
+    } while (denominator % dProbe === 0);
+    dProbe = lastProbe;
+  }
+
+  return [result, (n as Fraction).div(new Fraction(nProbe, dProbe))];
+}
+
+function bigIntToMonzoAndResidualLegacy(
+  n: bigint,
+  numberOfComponents: number
+): [Monzo, bigint] {
+  if (!n) {
+    return [Array(numberOfComponents).fill(0), 0n];
+  }
+
+  let probe = 1n;
+
+  const result = Array(numberOfComponents).fill(-1);
+  for (let i = 0; i < numberOfComponents; ++i) {
+    let lastProbe;
+    do {
+      lastProbe = probe;
+      probe *= BIG_INT_PRIMES[i];
+      result[i]++;
+    } while (n % probe === 0n);
+    probe = lastProbe;
+  }
+
+  return [result, n / probe];
+}