From 380f5e4faa6bed3b70a897a2f915d8ad931c56d2 Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Fri, 21 Jun 2024 14:02:14 +0300 Subject: [PATCH] Implement helpers related to lattice bases Variants for add, sub, dot and scale. Unnormalized Gram-Schmidt orthogonalization. LLL basis reduction. --- src/__tests__/basis.spec.ts | 283 ++++++++++++++++++++++++++++++++++++ src/__tests__/monzo.spec.ts | 31 ++++ src/basis.ts | 216 +++++++++++++++++++++++++++ src/index.ts | 67 +-------- src/monzo.ts | 167 ++++++++++++++++++++- src/number-array.ts | 64 ++++++++ 6 files changed, 757 insertions(+), 71 deletions(-) create mode 100644 src/__tests__/basis.spec.ts create mode 100644 src/basis.ts create mode 100644 src/number-array.ts diff --git a/src/__tests__/basis.spec.ts b/src/__tests__/basis.spec.ts new file mode 100644 index 0000000..fd524e9 --- /dev/null +++ b/src/__tests__/basis.spec.ts @@ -0,0 +1,283 @@ +import {describe, expect, it} from 'vitest'; +import { + fractionalGram, + fractionalLenstraLenstraLovasz, + gram, + lenstraLenstraLovasz, +} from '../basis'; +import { + applyWeights, + fractionalDot, + fractionalMonzosEqual, + monzoToFraction, + monzosEqual, + toMonzo, + unapplyWeights, +} from '../monzo'; +import {dot} from '../number-array'; +import {LOG_PRIMES} from '../primes'; + +const FUZZ = 'FUZZ' in process.env; + +describe('Gram process', () => { + it('orthogonalizes a basis', () => { + const basis = [ + [1, 2, 3], + [4, -5, 6], + [-7, 8, 9], + ]; + const {ortho, dual} = gram(basis); + // Leading orientation + expect(monzosEqual(basis[0], ortho[0])); + // Geometric duals + expect(dot(ortho[0], dual[0])).toBeCloseTo(1); + expect(dot(ortho[1], dual[1])).toBeCloseTo(1); + expect(dot(ortho[2], dual[2])).toBeCloseTo(1); + // Orthogonality + expect(dot(ortho[0], ortho[1])).toBeCloseTo(0); + expect(dot(ortho[0], ortho[2])).toBeCloseTo(0); + expect(dot(ortho[1], ortho[2])).toBeCloseTo(0); + // Value + expect(ortho.map(o => o.map(c => c.toFixed(2)))).toEqual([ + ['1.00', '2.00', '3.00'], + ['3.14', '-6.71', '3.43'], + ['-7.46', '-1.66', '3.59'], + ]); + }); + + it('handles non-basis', () => { + const basis = [ + [1, 2, 3, 4], + [0, 1, 1, 0], + [1, 1, 2, 4], + [-1, 2, 0, 0], + ]; + const {ortho} = gram(basis); + // Pseudo-orthogonality + expect(dot(ortho[0], ortho[1])).toBeCloseTo(0); + expect(dot(ortho[0], ortho[2])).toBeCloseTo(0); + expect(dot(ortho[0], ortho[3])).toBeCloseTo(0); + expect(dot(ortho[1], ortho[2])).toBeCloseTo(0); + expect(dot(ortho[1], ortho[3])).toBeCloseTo(0); + expect(dot(ortho[2], ortho[3])).toBeCloseTo(0); + }); + + it.runIf(FUZZ)('Fuzzes for random bases', () => { + for (let k = 0; k < 100000; ++k) { + const basis: number[][] = []; + for (let i = Math.random() * 10; i > 0; --i) { + const row: number[] = []; + for (let j = Math.random() * 10; j > 0; --j) { + row.push(Math.random() * 100 - 50); + } + basis.push(row); + } + const {ortho} = gram(basis); + // Pseudo-orthogonality + for (let i = 0; i < basis.length; ++i) { + for (let j = 0; j < i; ++j) { + expect(dot(ortho[i], ortho[j])).toBeCloseTo(0); + } + } + } + }); +}); + +describe('Gram process for arrays of fractions', () => { + it('orthogonalizes a basis', () => { + const basis = [ + [1, 2, 3], + [4, -5, 6], + [-7, 8, 9], + ]; + const {ortho, dual} = fractionalGram(basis); + // Leading orientation + expect(fractionalMonzosEqual(basis[0], ortho[0])); + // Geometric duals + expect(fractionalDot(ortho[0], dual[0]).toFraction()).toBe('1'); + expect(fractionalDot(ortho[1], dual[1]).toFraction()).toBe('1'); + expect(fractionalDot(ortho[2], dual[2]).toFraction()).toBe('1'); + // Orthogonality + expect(fractionalDot(ortho[0], ortho[1]).n).toBe(0); + expect(fractionalDot(ortho[0], ortho[2]).n).toBe(0); + expect(fractionalDot(ortho[1], ortho[2]).n).toBe(0); + // Value + expect(ortho.map(o => o.map(c => c.toFraction()))).toEqual([ + ['1', '2', '3'], + ['22/7', '-47/7', '24/7'], + ['-3483/467', '-774/467', '1677/467'], + ]); + }); + + it('handles non-basis', () => { + const basis = [ + [1, 2, 3, 4], + [0, 1, 1, 0], + [1, 1, 2, 4], + [-1, 2, 0, 0], + ]; + const {ortho} = fractionalGram(basis); + // Pseudo-orthogonality + expect(fractionalDot(ortho[0], ortho[1]).n).toBeCloseTo(0); + expect(fractionalDot(ortho[0], ortho[2]).n).toBeCloseTo(0); + expect(fractionalDot(ortho[0], ortho[3]).n).toBeCloseTo(0); + expect(fractionalDot(ortho[1], ortho[2]).n).toBeCloseTo(0); + expect(fractionalDot(ortho[1], ortho[3]).n).toBeCloseTo(0); + expect(fractionalDot(ortho[2], ortho[3]).n).toBeCloseTo(0); + }); +}); + +describe('LLL basis reduction', () => { + it('can LLL reduce', () => { + const basis = [ + [1, 1, 1], + [-1, 0, 2], + [3, 5, 6], + ]; + const lll = lenstraLenstraLovasz(basis); + // Size-reduction + for (let i = 0; i < 3; ++i) { + for (let j = 0; j < i; ++j) { + expect( + Math.abs(dot(lll.basis[i], lll.gram.dual[j])) + ).toBeLessThanOrEqual(0.5); + } + } + // Lovász condition + for (let k = 1; k < 3; ++k) { + const ok = lll.gram.ortho[k]; + const ok1 = lll.gram.ortho[k - 1]; + const mu = dot(lll.basis[k], lll.gram.dual[k - 1]); + const n1 = dot(ok1, ok1); + expect((n1 * 3) / 4).toBeLessThanOrEqual(dot(ok, ok) + n1 * mu * mu); + } + + expect(lll.basis).toEqual([ + [0, 1, 0], + [1, 0, 1], + [-1, 0, 2], + ]); + }); + + it('handles non-basis', () => { + const basis = [ + [1, 2, 3, 4], + [0, 1, 1, 0], + [1, 1, 2, 4], + [-1, 2, 0, 0], + ]; + const lll = lenstraLenstraLovasz(basis); + expect(lll.basis).toEqual([ + [0, 0, 0, 0], + [0, 1, 1, 0], + [-1, 1, -1, 0], + [0, 0, -1, 4], + ]); + }); + + it('can mess up the basis of miracle with naïve weights', () => { + const basis = ['225/224', '1029/1024'].map(toMonzo); + const lll = lenstraLenstraLovasz(basis); + expect(lll.basis.map(m => monzoToFraction(m).toFraction())).toEqual([ + '225/224', + '2401/2400', + ]); + }); + + it('can fix the basis of miracle with Tenney weights', () => { + const basis = ['225/224', '2401/2400'] + .map(toMonzo) + .map(m => applyWeights(m, LOG_PRIMES)); + const lll = lenstraLenstraLovasz(basis); + const commas = lll.basis + .map(m => unapplyWeights(m, LOG_PRIMES).map(Math.round)) + .map(m => monzoToFraction(m).toFraction()); + expect(commas).toEqual(['225/224', '1029/1024']); + }); + + it.runIf(FUZZ)('Fuzzes for random bases', () => { + for (let k = 0; k < 1000; ++k) { + let basis: number[][] = []; + for (let i = Math.random() * 10; i > 0; --i) { + const row: number[] = []; + for (let j = Math.random() * 10; j > 0; --j) { + row.push(Math.random() * 100 - 50); + } + basis.push(row); + } + if (Math.random() < 0.5) { + basis = basis.map(row => row.map(Math.round)); + } + const lll = lenstraLenstraLovasz(basis); + if (lll.gram.squaredLengths.every(l => l)) { + // Size-reduction + for (let i = 0; i < basis.length; ++i) { + for (let j = 0; j < i; ++j) { + expect( + Math.abs(dot(lll.basis[i], lll.gram.dual[j])) + ).toBeLessThanOrEqual(0.5); + } + } + // Lovász condition + for (let k = 1; k < basis.length; ++k) { + const ok = lll.gram.ortho[k]; + const ok1 = lll.gram.ortho[k - 1]; + const mu = dot(lll.basis[k], lll.gram.dual[k - 1]); + const n1 = dot(ok1, ok1); + expect((n1 * 3) / 4).toBeLessThanOrEqual(dot(ok, ok) + n1 * mu * mu); + } + } + } + }); +}); + +describe('Precise LLL basis reduction', () => { + it('can LLL reduce', () => { + const basis = [ + [1, 1, 1], + [-1, 0, 2], + [3, 5, 6], + ]; + const lll = fractionalLenstraLenstraLovasz(basis); + // Size-reduction + for (let i = 0; i < 3; ++i) { + for (let j = 0; j < i; ++j) { + expect( + fractionalDot(lll.basis[i], lll.gram.dual[j]).compare(0.5) + ).toBeLessThanOrEqual(0); + } + } + // Lovász condition + for (let k = 1; k < 3; ++k) { + const ok = lll.gram.ortho[k]; + const ok1 = lll.gram.ortho[k - 1]; + const mu = fractionalDot(lll.basis[k], lll.gram.dual[k - 1]); + const n1 = fractionalDot(ok1, ok1); + expect( + n1.mul('3/4').compare(fractionalDot(ok, ok).add(n1.mul(mu.mul(mu)))) + ).toBeLessThanOrEqual(0); + } + + expect(lll.basis.map(row => row.map(f => f.valueOf()))).toEqual([ + [0, 1, 0], + [1, 0, 1], + [-1, 0, 2], + ]); + }); + + it('handles non-basis', () => { + const basis = [ + [1, 2, 3, 4], + [0, 1, 1, 0], + [1, 1, 2, 4], + [-1, 2, 0, 0], + ]; + const lll = fractionalLenstraLenstraLovasz(basis); + expect(lll.basis.map(row => row.map(f => f.valueOf()))).toEqual([ + [0, 0, 0, 0], + [0, 1, 1, 0], + [-1, 1, -1, 0], + [0, 0, -1, 4], + ]); + }); +}); diff --git a/src/__tests__/monzo.spec.ts b/src/__tests__/monzo.spec.ts index 6ab2b2f..87b710a 100644 --- a/src/__tests__/monzo.spec.ts +++ b/src/__tests__/monzo.spec.ts @@ -1,6 +1,9 @@ import {describe, it, expect} from 'vitest'; import {Fraction} from '../fraction'; import { + fractionalAdd, + fractionalMonzosEqual, + fractionalNorm, monzoToBigInt, monzoToFraction, primeFactorize, @@ -468,3 +471,31 @@ describe('Sparse monzos', () => { } }); }); + +describe('Fractional monzo methods', () => { + it('test for equality between two monzos (equal)', () => { + const yes = fractionalMonzosEqual( + ['1/2', '7/9'], + [0.5, new Fraction(14, 18), 0] + ); + expect(yes).toBe(true); + }); + + it('test for equality between two monzos (not equal)', () => { + const no = fractionalMonzosEqual( + ['1/2', '7/9'], + [0.75, new Fraction(7, 9)] + ); + expect(no).toBe(false); + }); + + it('adds two fractional monzos', () => { + const result = fractionalAdd(['1/2', '2/3'], [new Fraction(1), 0.75]); + expect(result.map(f => f.toFraction())).toEqual(['3/2', '17/12']); + }); + + it('measures the naïve squared length of a fractional monzo', () => { + const l2 = fractionalNorm([0.5, '1/3', '5/7']); + expect(l2.toFraction()).toBe('1537/1764'); + }); +}); diff --git a/src/basis.ts b/src/basis.ts new file mode 100644 index 0000000..f6088c5 --- /dev/null +++ b/src/basis.ts @@ -0,0 +1,216 @@ +import {Fraction, FractionValue} from './fraction'; +import { + FractionalMonzo, + ProtoFractionalMonzo, + fractionalDot, + fractionalScale, + fractionalSub, + scale, + sub, +} from './monzo'; +import {dot} from './number-array'; + +/** + * Result of Gram–Schmidt process without normalization. + */ +export type GramResult = { + /** Orthogonal basis with the leading basis element intact. */ + ortho: number[][]; + /** Squared lengths of the orthogonal basis. */ + squaredLengths: number[]; + /** Geometric duals of the orthogonal basis. */ + dual: number[][]; +}; + +/** + * Result of Gram–Schmidt process without normalization. + */ +export type FractionalGramResult = { + /** Orthogonal basis with the leading basis element intact. */ + ortho: FractionalMonzo[]; + /** Squared lengths of the orthogonal basis. */ + squaredLengths: Fraction[]; + /** Geometric duals of the orthogonal basis. */ + dual: FractionalMonzo[]; +}; + +/** + * Result of Lenstra-Lenstra-Lovász basis reduction. + */ +export type LLLResult = { + /** Basis that's short and nearly orthogonal. */ + basis: number[][]; + /** Gram-Schmidt process results. */ + gram: GramResult; +}; + +/** + * Result of Lenstra-Lenstra-Lovász basis reduction. + */ +export type FractionalLLLResult = { + /** Basis that's short and nearly orthogonal. */ + basis: FractionalMonzo[]; + /** Gram-Schmidt process results. */ + gram: FractionalGramResult; +}; + +/** + * Perform Gram–Schmidt process without normalization. + * @param basis Array of basis elements. + * @param epsilon Threshold for zero. + * @returns The orthogonalized basis and its dual (duals of near-zero basis elements are coerced to zero). + */ +export function gram(basis: number[][], epsilon = 1e-12): GramResult { + const ortho: number[][] = []; + const squaredLengths: number[] = []; + const dual: number[][] = []; + for (let i = 0; i < basis.length; ++i) { + ortho.push([...basis[i]]); + for (let j = 0; j < i; ++j) { + ortho[i] = sub(ortho[i], scale(ortho[j], dot(ortho[i], dual[j]))); + } + squaredLengths.push(dot(ortho[i], ortho[i])); + if (squaredLengths[i] > epsilon) { + dual.push(scale(ortho[i], 1 / squaredLengths[i])); + } else { + dual.push(scale(ortho[i], 0)); + } + } + return {ortho, squaredLengths, dual}; +} + +/** + * Perform Gram–Schmidt process without normalization. + * @param basis Array of rational basis elements. + * @returns The orthogonalized basis and its dual as arrays of fractions (duals of zero basis elements are coerced to zero). + */ +export function fractionalGram( + basis: ProtoFractionalMonzo[] +): FractionalGramResult { + const ortho: Fraction[][] = []; + const squaredLengths: Fraction[] = []; + const dual: Fraction[][] = []; + for (let i = 0; i < basis.length; ++i) { + ortho.push(basis[i].map(f => new Fraction(f))); + for (let j = 0; j < i; ++j) { + ortho[i] = fractionalSub( + ortho[i], + fractionalScale(ortho[j], fractionalDot(ortho[i], dual[j])) + ); + } + squaredLengths.push(fractionalDot(ortho[i], ortho[i])); + if (squaredLengths[i].n) { + dual.push(fractionalScale(ortho[i], squaredLengths[i].inverse())); + } else { + dual.push(fractionalScale(ortho[i], squaredLengths[i])); + } + } + return {ortho, squaredLengths, dual}; +} + +/** + * Preform Lenstra-Lenstra-Lovász basis reduction. + * @param basis Array of basis elements. + * @param delta Lovász coefficient. + * @param epsilon Threshold for zero. + * @param maxIterations Maximum number of iterations to perform. + * @returns The basis processed to be short and nearly orthogonal alongside Gram-Schmidt coefficients. + */ +export function lenstraLenstraLovasz( + basis: number[][], + delta = 0.75, + epsilon = 1e-12, + maxIterations = 10000 +): LLLResult { + // https://en.wikipedia.org/wiki/Lenstra%E2%80%93Lenstra%E2%80%93Lov%C3%A1sz_lattice_basis_reduction_algorithm#LLL_algorithm_pseudocode + basis = basis.map(row => [...row]); + let {ortho, squaredLengths, dual} = gram(basis, epsilon); + let k = 1; + while (k < basis.length && maxIterations--) { + for (let j = k - 1; j >= 0; --j) { + const mu = dot(basis[k], dual[j]); + if (Math.abs(mu) > 0.5) { + basis[k] = sub(basis[k], scale(basis[j], Math.round(mu))); + + ({ortho, squaredLengths, dual} = gram(basis, epsilon)); + } + } + const mu = dot(basis[k], dual[k - 1]); + if (squaredLengths[k] > (delta - mu * mu) * squaredLengths[k - 1]) { + k++; + } else { + const bk = basis[k]; + basis[k] = basis[k - 1]; + basis[k - 1] = bk; + + ({ortho, squaredLengths, dual} = gram(basis, epsilon)); + + k = Math.max(k - 1, 1); + } + } + return { + basis, + gram: { + ortho, + squaredLengths, + dual, + }, + }; +} + +const HALF = new Fraction(1, 2); + +/** + * Preform Lenstra-Lenstra-Lovász basis reduction using rational numbers. + * @param basis Array of rational basis elements. + * @param delta Lovász coefficient. + * @param maxIterations Maximum number of iterations to perform. + * @returns The basis processed to be short and nearly orthogonal alongside Gram-Schmidt coefficients. + */ +export function fractionalLenstraLenstraLovasz( + basis: ProtoFractionalMonzo[], + delta: FractionValue = '3/4', + maxIterations = 10000 +): FractionalLLLResult { + const result = basis.map(row => row.map(f => new Fraction(f))); + const delta_ = new Fraction(delta); + let {ortho, squaredLengths, dual} = fractionalGram(result); + let k = 1; + while (k < result.length && maxIterations--) { + for (let j = k - 1; j >= 0; --j) { + const mu = fractionalDot(result[k], dual[j]); + if (mu.abs().compare(HALF) > 0) { + result[k] = fractionalSub( + result[k], + fractionalScale(result[j], mu.round()) + ); + + ({ortho, squaredLengths, dual} = fractionalGram(result)); + } + } + const mu = fractionalDot(result[k], dual[k - 1]); + if ( + squaredLengths[k].compare( + delta_.sub(mu.mul(mu)).mul(squaredLengths[k - 1]) + ) > 0 + ) { + k++; + } else { + const bk = result[k]; + result[k] = result[k - 1]; + result[k - 1] = bk; + + ({ortho, squaredLengths, dual} = fractionalGram(result)); + + k = Math.max(k - 1, 1); + } + } + return { + basis: result, + gram: { + ortho, + squaredLengths, + dual, + }, + }; +} diff --git a/src/index.ts b/src/index.ts index df58177..29ddb40 100644 --- a/src/index.ts +++ b/src/index.ts @@ -1,7 +1,7 @@ import {Fraction, FractionValue, mmod} from './fraction'; import {Monzo, monzoToBigNumeratorDenominator, primeFactorize} from './monzo'; +import {dot, dotPrecise} from './number-array'; import {LOG_PRIMES, PRIMES, PRIME_CENTS} from './primes'; -import {sum} from './polyfills/sum-precise'; export * from './fraction'; export * from './primes'; @@ -9,6 +9,8 @@ export * from './conversion'; export * from './combinations'; export * from './monzo'; export * from './approximation'; +export * from './number-array'; +export * from './basis'; export {sum} from './polyfills/sum-precise'; export interface AnyArray { @@ -16,11 +18,6 @@ export interface AnyArray { length: number; } -export interface NumberArray { - [key: number]: number; - length: number; -} - /** * Check if the contents of two arrays are equal using '==='. * @param a The first array. @@ -243,64 +240,6 @@ export function clamp(minValue: number, maxValue: number, value: number) { return value; } -/** - * Calculate the inner (dot) product of two arrays of real numbers. - * @param a The first array of numbers. - * @param b The second array of numbers. - * @returns The dot product. - */ -export function dot(a: NumberArray, b: NumberArray): number { - const numComponents = Math.min(a.length, b.length); - let result = 0; - for (let i = 0; i < numComponents; ++i) { - result += a[i] * b[i]; - } - return result; -} - -/** - * Calculate the inner (dot) product of two arrays of real numbers. - * The resulting terms are summed accurately using Shewchuk's algorithm. - * @param a The first array of numbers. - * @param b The second array of numbers. - * @returns The dot product. - */ -export function dotPrecise(a: NumberArray, b: NumberArray): number { - const numComponents = Math.min(a.length, b.length); - function* terms() { - for (let i = 0; i < numComponents; ++i) { - yield a[i] * b[i]; - } - } - return sum(terms()); -} - -/** - * Calculate the norm (vector length) of an array of real numbers. - * @param array The array to measure. - * @param type Type of measurement. - * @returns The length of the vector. - */ -export function norm( - array: NumberArray, - type: 'euclidean' | 'taxicab' | 'maximum' = 'euclidean' -) { - let result = 0; - for (let i = 0; i < array.length; ++i) { - if (type === 'taxicab') { - result += Math.abs(array[i]); - } else if (type === 'maximum') { - result = Math.max(result, Math.abs(array[i])); - } else { - result += array[i] * array[i]; - } - } - if (type === 'euclidean') { - return Math.sqrt(result); - } - return result; -} - /** * Calculate the difference between two cents values such that equave equivalence is taken into account. * @param a The first pitch measured in cents. diff --git a/src/monzo.ts b/src/monzo.ts index 175a186..6501a3a 100644 --- a/src/monzo.ts +++ b/src/monzo.ts @@ -6,6 +6,16 @@ import {BIG_INT_PRIMES, PRIMES} from './primes'; */ export type Monzo = number[]; +/** + * Array of rational numbers representing the exponents of prime numbers in the unique factorization of a radical number i.e. n-th root. + */ +export type FractionalMonzo = Fraction[]; + +/** + * Array of rational-looking numbers suitable as arguments for methods that work with fractional monzos. + */ +export type ProtoFractionalMonzo = FractionValue[]; + // The limit at which ((n ^ (n-1)) & n) is no longer equal to the two's factor. const BIT_MAGIC_LIMIT = 2 ** 31; @@ -28,7 +38,7 @@ export function monzosEqual(a: Monzo, b: Monzo) { if (a === b) { return true; } - for (let i = 0; i < Math.min(a.length, b.length); ++i) { + for (let i = Math.min(a.length, b.length) - 1; i >= 0; --i) { if (a[i] !== b[i]) { return false; } @@ -50,6 +60,41 @@ export function monzosEqual(a: Monzo, b: Monzo) { return true; } +/** + * Check if two fractional monzos are equal. + * @param a The first monzo. + * @param b The second monzo. + * @returns `true` if the two values are equal. + */ +export function fractionalMonzosEqual( + a: ProtoFractionalMonzo, + b: ProtoFractionalMonzo +): boolean { + if (a === b) { + return true; + } + for (let i = Math.min(a.length, b.length) - 1; i >= 0; --i) { + if (!new Fraction(a[i]).equals(b[i])) { + return false; + } + } + if (a.length > b.length) { + for (let i = b.length; i < a.length; ++i) { + if (new Fraction(a[i]).n) { + return false; + } + } + } + if (b.length > a.length) { + for (let i = a.length; i < b.length; ++i) { + if (new Fraction(b[i]).n) { + return false; + } + } + } + return true; +} + /** * Add two monzos. * @param a The first monzo. @@ -57,13 +102,31 @@ export function monzosEqual(a: Monzo, b: Monzo) { * @returns A monzo that represents the product of the two numbers represented by `a` and `b`. */ export function add(a: Monzo, b: Monzo): Monzo { - if (a.length < b.length) { - return add(b, a); - } const result = [...a]; - for (let i = 0; i < b.length; ++i) { + for (let i = Math.min(a.length, b.length) - 1; i >= 0; --i) { result[i] += b[i]; } + result.push(...b.slice(result.length)); + return result; +} + +/** + * Add two fractional monzos. + * @param a The first monzo. + * @param b The second monzo. + * @returns A monzo that represents the product of the two numbers represented by `a` and `b`. + */ +export function fractionalAdd( + a: ProtoFractionalMonzo, + b: ProtoFractionalMonzo +): FractionalMonzo { + const result: FractionalMonzo = a.map(f => new Fraction(f)); + for (let i = Math.min(a.length, b.length) - 1; i >= 0; --i) { + result[i] = result[i].add(b[i]); + } + while (result.length < b.length) { + result.push(new Fraction(b[result.length])); + } return result; } @@ -75,7 +138,7 @@ export function add(a: Monzo, b: Monzo): Monzo { */ export function sub(a: Monzo, b: Monzo): Monzo { const result = [...a]; - for (let i = 0; i < Math.min(a.length, b.length); ++i) { + for (let i = Math.min(a.length, b.length) - 1; i >= 0; --i) { result[i] -= b[i]; } while (result.length < b.length) { @@ -84,6 +147,26 @@ export function sub(a: Monzo, b: Monzo): Monzo { return result; } +/** + * Subtract two fractional monzos. + * @param a The first monzo. + * @param b The second monzo. + * @returns A monzo that represents the division of the two numbers represented by `a` and `b`. + */ +export function fractionalSub( + a: ProtoFractionalMonzo, + b: ProtoFractionalMonzo +): FractionalMonzo { + const result: FractionalMonzo = a.map(f => new Fraction(f)); + for (let i = Math.min(a.length, b.length) - 1; i >= 0; --i) { + result[i] = result[i].sub(b[i]); + } + while (result.length < b.length) { + result.push(new Fraction(b[result.length]).neg()); + } + return result; +} + /** * Scale a monzo by a scalar. * @param monzo The monzo to scale. @@ -94,6 +177,62 @@ export function scale(monzo: Monzo, amount: number) { return monzo.map(component => component * amount); } +/** + * Scale a monzo by a scalar. + * @param monzo The monzo to scale. + * @param amount The amount to scale by. + * @returns The scalar multiple. + */ +export function fractionalScale( + monzo: ProtoFractionalMonzo, + amount: FractionValue +): FractionalMonzo { + return monzo.map(component => new Fraction(component).mul(amount)); +} + +/** + * Calculate the inner (dot) product of two arrays of rational numbers. + * @param a The first array of numbers. + * @param b The second array of numbers. + * @returns The dot product. + */ +export function fractionalDot( + a: ProtoFractionalMonzo, + b: ProtoFractionalMonzo +): Fraction { + let result = new Fraction(0); + for (let i = Math.min(a.length, b.length) - 1; i >= 0; --i) { + result = result.add(new Fraction(a[i]).mul(b[i])); + } + return result; +} + +/** + * Calculate the norm (vector length) of an array of rational numbers. + * @param array The array to measure. + * @param type Type of measurement. (Euclidean norm can be obtained using L2 and calling .sqrt() on the result.) + * @returns The length of the vector. + */ +export function fractionalNorm( + array: ProtoFractionalMonzo, + type: 'L2' | 'taxicab' | 'maximum' = 'L2' +): Fraction { + let result = new Fraction(0); + for (let i = 0; i < array.length; ++i) { + const a = new Fraction(array[i]).abs(); + if (type === 'taxicab') { + result = result.add(a); + } else if (type === 'maximum') { + if (result.compare(a) < 0) { + result = a; + } + } else { + result = result.add(a.mul(a)); + } + } + return result; +} + /** * Multiply two monzos component-wise. * @param monzo The first monzo. @@ -102,12 +241,26 @@ export function scale(monzo: Monzo, amount: number) { */ export function applyWeights(monzo: Monzo, weights: Monzo) { const result = [...monzo]; - for (let i = 0; i < Math.min(monzo.length, weights.length); ++i) { + for (let i = Math.min(monzo.length, weights.length) - 1; i >= 0; --i) { result[i] *= weights[i]; } return result; } +/** + * Divide two monzos component-wise. + * @param monzo The first monzo. + * @param weights The second monzo. Missing values interpreted as 1 (no change). + * @returns The first monzo unweighted by the second. + */ +export function unapplyWeights(monzo: Monzo, weights: Monzo) { + const result = [...monzo]; + for (let i = Math.min(monzo.length, weights.length) - 1; i >= 0; --i) { + result[i] /= weights[i]; + } + return result; +} + /** * Accumulate a monzo into the first one. * @param target The monzo to accumulate into. diff --git a/src/number-array.ts b/src/number-array.ts new file mode 100644 index 0000000..4afa9a5 --- /dev/null +++ b/src/number-array.ts @@ -0,0 +1,64 @@ +import {sum} from './polyfills/sum-precise'; + +export interface NumberArray { + [key: number]: number; + length: number; +} + +/** + * Calculate the inner (dot) product of two arrays of real numbers. + * @param a The first array of numbers. + * @param b The second array of numbers. + * @returns The dot product. + */ +export function dot(a: NumberArray, b: NumberArray): number { + const numComponents = Math.min(a.length, b.length); + let result = 0; + for (let i = 0; i < numComponents; ++i) { + result += a[i] * b[i]; + } + return result; +} + +/** + * Calculate the inner (dot) product of two arrays of real numbers. + * The resulting terms are summed accurately using Shewchuk's algorithm. + * @param a The first array of numbers. + * @param b The second array of numbers. + * @returns The dot product. + */ +export function dotPrecise(a: NumberArray, b: NumberArray): number { + const numComponents = Math.min(a.length, b.length); + function* terms() { + for (let i = 0; i < numComponents; ++i) { + yield a[i] * b[i]; + } + } + return sum(terms()); +} + +/** + * Calculate the norm (vector length) of an array of real numbers. + * @param array The array to measure. + * @param type Type of measurement. + * @returns The length of the vector. + */ +export function norm( + array: NumberArray, + type: 'euclidean' | 'L2' | 'taxicab' | 'maximum' = 'euclidean' +) { + let result = 0; + for (let i = 0; i < array.length; ++i) { + if (type === 'taxicab') { + result += Math.abs(array[i]); + } else if (type === 'maximum') { + result = Math.max(result, Math.abs(array[i])); + } else { + result += array[i] * array[i]; + } + } + if (type === 'euclidean') { + return Math.sqrt(result); + } + return result; +}