diff --git a/src/__tests__/monzo.spec.ts b/src/__tests__/monzo.spec.ts index a6bad7b..40534d4 100644 --- a/src/__tests__/monzo.spec.ts +++ b/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); + } + }); +}); diff --git a/src/monzo.ts b/src/monzo.ts index 8a298f2..93a8196 100644 --- a/src/monzo.ts +++ b/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 { + 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 = 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; +}