Implement integer factorization into primes

ref #34

src/__tests__/monzo.spec.ts

@@ -3,10 +3,12 @@ import {Fraction} from '../fraction';
 import {
   monzoToBigInt,
   monzoToFraction,
+  primeFactorize,
   primeLimit,
   toMonzo,
   toMonzoAndResidual,
 } from '../monzo';
+import {isPrime} from '../primes';
 
 function toMonzoAndResidual11(n: number): [number[], number] {
   const result = [0, 0, 0, 0, 0];
@@ -130,6 +132,12 @@ describe('Monzo converter', () => {
       }
     }
   });
+
+  it('refuses to factor a negative fraction', () => {
+    expect(() => toMonzo('-1/2')).toThrow(
+      'Cannot convert fraction -1/2 to monzo'
+    );
+  });
 });
 
 describe('Fraction to monzo converter', () => {
@@ -372,3 +380,72 @@ describe('Prime limit calculator', () => {
     expect(primeLimit(2n ** 1025n)).toEqual(2);
   });
 });
+
+describe('Sparse monzos', () => {
+  it('factorizes 12', () => {
+    const exponentByPrime = primeFactorize(12);
+    expect(exponentByPrime).toHaveLength(2);
+    expect(exponentByPrime.get(2)).toBe(2);
+    expect(exponentByPrime.get(3)).toBe(1);
+  });
+
+  it('factorizes 0', () => {
+    const factors = primeFactorize(0);
+    expect(factors).toHaveLength(1);
+    expect(factors.get(0)).toBe(1);
+  });
+
+  it('factorizes -35', () => {
+    const factors = primeFactorize(-35);
+    expect(factors).toHaveLength(3);
+    expect(factors.get(-1)).toBe(1);
+    expect(factors.get(5)).toBe(1);
+    expect(factors.get(7)).toBe(1);
+  });
+
+  it('factorizes 81/80', () => {
+    const factors = primeFactorize('81/80');
+    expect(factors).toHaveLength(3);
+    expect(factors.get(2)).toBe(-4);
+    expect(factors.get(3)).toBe(4);
+    expect(factors.get(5)).toBe(-1);
+  });
+
+  it('factorizes 1073741823', () => {
+    const factors = primeFactorize(1073741823);
+    expect(factors).toHaveLength(6);
+    expect(factors.get(3)).toBe(2);
+    expect(factors.get(7)).toBe(1);
+    expect(factors.get(11)).toBe(1);
+    expect(factors.get(31)).toBe(1);
+    expect(factors.get(151)).toBe(1);
+    expect(factors.get(331)).toBe(1);
+  });
+
+  it.each([
+    49305423, 4104956, 8375509, 27826943, 44852222, 22932439, 46933379,
+    59598447, 9693451, 54191546, 61834729, 26866018, 46410510, 47335837,
+    43839566,
+  ])('works on a tough case %s found during fuzzing', n => {
+    const exponentByPrime = primeFactorize(n);
+    let m = 1;
+    for (const prime of exponentByPrime.keys()) {
+      expect(isPrime(prime)).toBe(true);
+      m *= prime ** exponentByPrime.get(prime)!;
+    }
+    expect(m).toBe(n);
+  });
+
+  it.skip('fuzzes for more broken cases', () => {
+    for (let i = 0; i < 100; ++i) {
+      const n = Math.floor(Math.random() * 62837328) + 1;
+      const exponentByPrime = primeFactorize(n);
+      let m = 1;
+      for (const prime of exponentByPrime.keys()) {
+        expect(isPrime(prime)).toBe(true);
+        m *= prime ** exponentByPrime.get(prime)!;
+      }
+      expect(m).toBe(n);
+    }
+  });
+});

src/monzo.ts

@@ -1,4 +1,4 @@
-import {Fraction, FractionValue} from './fraction';
+import {Fraction, FractionValue, gcd} from './fraction';
 import {BIG_INT_PRIMES, PRIMES} from './primes';
 
 /**
@@ -158,6 +158,11 @@ export function toMonzo(n: FractionValue | bigint): Monzo {
   }
   if (typeof n !== 'number') {
     n = new Fraction(n);
+    if ((n as Fraction).s !== 1) {
+      throw new Error(
+        `Cannot convert fraction ${(n as Fraction).toFraction()} to monzo`
+      );
+    }
     return sub(toMonzo(n.n), toMonzo(n.d));
   }
   if (n < 1 || Math.round(n) !== n) {
@@ -614,3 +619,113 @@ function bigIntToMonzo7(n: bigint): [Monzo, bigint] {
   }
   return [result, n];
 }
+
+// Condition: m, n < 2**30
+function modMul(m: number, n: number, modulus: number) {
+  let result = 0;
+  let current = m;
+  while (n) {
+    if (n & 1) {
+      result = (result + current) % modulus;
+    }
+    current = (current << 1) % modulus;
+    n >>= 1;
+  }
+  return result;
+}
+
+function pollardRhoStep(n: number, p: number, seed = 1) {
+  // n² + s mod p
+  return modMul(n, n, p) + seed;
+}
+
+// Condition: n !== 1
+function pollardRhoFactor(n: number, seed = 2, stepSeed = 1) {
+  let x = seed;
+  let y = x;
+  let d = 1;
+  while (d === 1) {
+    x = pollardRhoStep(x, n, stepSeed);
+    y = pollardRhoStep(pollardRhoStep(y, n, stepSeed), n, stepSeed);
+    d = gcd(Math.abs(x - y), n);
+  }
+  return d;
+}
+
+// Condition n ∈ 7-limit residue
+// Checked up to 20000000.
+function rhoCascade(n: number) {
+  let factor = pollardRhoFactor(n);
+  if (factor === n) {
+    factor = pollardRhoFactor(n, 3);
+    if (factor === n) {
+      factor = pollardRhoFactor(n, 9);
+      if (factor === n) {
+        factor = pollardRhoFactor(n, 4, 2);
+        if (factor === n) {
+          return pollardRhoFactor(n, 1, 2);
+        }
+      }
+    }
+  }
+  return factor;
+}
+
+/**
+ * Factorize a number into a `Map` instace with prime numbers as keys and their multiplicity as values.
+ * @param value Rational number to factorize.
+ * @returns A sparse monzo.
+ */
+export function primeFactorize(value: FractionValue): Map<number, number> {
+  if (typeof value !== 'number' || !Number.isInteger(value)) {
+    const {s, n, d} = new Fraction(value);
+    const nResult = primeFactorize(s * n);
+    const dResult = primeFactorize(d);
+    for (const [prime, exponent] of dResult) {
+      nResult.set(prime, -exponent);
+    }
+    return nResult;
+  }
+  const result: Map<number, number> = new Map();
+  if (value === 0) {
+    result.set(0, 1);
+    return result;
+  } else if (value < 0) {
+    result.set(-1, 1);
+    value = -value;
+  }
+  if (value > 1073741823) {
+    throw new Error('Factorization not implemented above 1073741823.');
+  }
+  let [monzo, residual] = intToMonzo7(value);
+  for (let i = 0; i < monzo.length; ++i) {
+    if (monzo[i]) {
+      result.set(PRIMES[i], monzo[i]);
+    }
+  }
+  // This is entirely ad. hoc. with holes patched as they came up during fuzzing.
+  while (residual !== 1) {
+    let factor = rhoCascade(residual);
+    residual /= factor;
+    let subFactor = rhoCascade(factor);
+    if (subFactor === factor) {
+      result.set(factor, (result.get(factor) ?? 0) + 1);
+    } else {
+      while (factor !== subFactor) {
+        if (subFactor === 901) {
+          result.set(17, (result.get(17) ?? 0) + 1);
+          result.set(53, (result.get(53) ?? 0) + 1);
+        } else if (subFactor === 1241) {
+          result.set(17, (result.get(17) ?? 0) + 1);
+          result.set(73, (result.get(73) ?? 0) + 1);
+        } else {
+          result.set(subFactor, (result.get(subFactor) ?? 0) + 1);
+        }
+        factor /= subFactor;
+        subFactor = rhoCascade(factor);
+      }
+      result.set(factor, (result.get(factor) ?? 0) + 1);
+    }
+  }
+  return result;
+}