From 598a58d5a25aeb29b963eee8eb9e3937c4eb7318 Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Mon, 25 Dec 2023 09:20:28 +0200 Subject: [PATCH] Move approximation methods to a separate file --- src/__tests__/approximation.spec.ts | 116 ++++++++++++ src/__tests__/index.spec.ts | 112 ------------ src/approximation.ts | 271 +++++++++++++++++++++++++++ src/index.ts | 273 +--------------------------- 4 files changed, 389 insertions(+), 383 deletions(-) create mode 100644 src/__tests__/approximation.spec.ts create mode 100644 src/approximation.ts diff --git a/src/__tests__/approximation.spec.ts b/src/__tests__/approximation.spec.ts new file mode 100644 index 0000000..2e75780 --- /dev/null +++ b/src/__tests__/approximation.spec.ts @@ -0,0 +1,116 @@ +import {describe, it, expect} from 'vitest'; +import { + approximateOddLimit, + approximatePrimeLimit, + approximatePrimeLimitWithErrors, + getConvergents, +} from '../approximation'; +import {valueToCents} from '../conversion'; +import {PRIMES} from '../primes'; + +describe('Odd limit approximator', () => { + it('can approximate tau in the 15-odd-limit', () => { + const approximation = approximateOddLimit(valueToCents(2 * Math.PI), 15)[0]; + expect(approximation.equals('44/7')).toBeTruthy(); + expect(approximation.valueOf()).toBeCloseTo(2 * Math.PI); + }); + + it('can approximate e in the 21-odd-limit', () => { + const approximations = approximateOddLimit(valueToCents(Math.E), 21); + expect(approximations[0].equals('19/7')).toBeTruthy(); + expect(approximations[0].valueOf()).toBeCloseTo(Math.E); + expect(approximations[7].equals('21/8')).toBeTruthy(); + expect(approximations[7].valueOf()).toBeCloseTo(Math.E, 0); + }); +}); + +describe('Prime limit approximator', () => { + it('can approximate pi in the 11-limit', () => { + const approximation = approximatePrimeLimit( + valueToCents(Math.PI), + PRIMES.indexOf(11), + 3 + )[0]; + expect(approximation.equals('12544/3993')).toBeTruthy(); + expect(approximation.valueOf()).toBeCloseTo(Math.PI); + }); + + it('can approximate pi in the 13-limit with a small sized result', () => { + const approximations = approximatePrimeLimit( + valueToCents(Math.PI), + PRIMES.indexOf(13), + 3, + 15, + 4 + ); + expect(approximations).toHaveLength(4); + }); + + it('can approximate the square root of two in the 7-limit within maximum error', () => { + const approximationsAndErrors = approximatePrimeLimitWithErrors( + 600, + PRIMES.indexOf(7), + 5, + 10 + ); + expect(approximationsAndErrors).toHaveLength(28); + approximationsAndErrors.forEach(([approximation, error]) => { + const cents = valueToCents(approximation.valueOf()); + const calculatedError = Math.abs(cents - 600); + expect(error).toBeCloseTo(calculatedError); + expect(calculatedError).toBeLessThanOrEqual(10); + }); + }); +}); + +describe('Convergent calculator', () => { + it('calculates the convergents of pi', () => { + const convergents = getConvergents(Math.PI, undefined, 10); + expect(convergents).toHaveLength(10); + expect(convergents[0].equals(3)).toBeTruthy(); + expect(convergents[1].equals('22/7')).toBeTruthy(); + expect(convergents[2].equals('333/106')).toBeTruthy(); + expect(convergents[3].equals('355/113')).toBeTruthy(); + expect(convergents[4].equals('103993/33102')).toBeTruthy(); + expect(convergents[5].equals('104348/33215')).toBeTruthy(); + expect(convergents[6].equals('208341/66317')).toBeTruthy(); + expect(convergents[7].equals('312689/99532')).toBeTruthy(); + expect(convergents[8].equals('833719/265381')).toBeTruthy(); + expect(convergents[9].equals('1146408/364913')).toBeTruthy(); + }); + + it('calculates the semiconvergents of pi', () => { + const semiconvergents = getConvergents(Math.PI, undefined, 13, true); + expect(semiconvergents).toHaveLength(13); + expect(semiconvergents[0].equals(3)).toBeTruthy(); + expect(semiconvergents[1].equals('13/4')).toBeTruthy(); + expect(semiconvergents[2].equals('16/5')).toBeTruthy(); + expect(semiconvergents[3].equals('19/6')).toBeTruthy(); + expect(semiconvergents[4].equals('22/7')).toBeTruthy(); + expect(semiconvergents[5].equals('179/57')).toBeTruthy(); + expect(semiconvergents[6].equals('201/64')).toBeTruthy(); + expect(semiconvergents[7].equals('223/71')).toBeTruthy(); + expect(semiconvergents[8].equals('245/78')).toBeTruthy(); + expect(semiconvergents[9].equals('267/85')).toBeTruthy(); + expect(semiconvergents[10].equals('289/92')).toBeTruthy(); + expect(semiconvergents[11].equals('311/99')).toBeTruthy(); + expect(semiconvergents[12].equals('333/106')).toBeTruthy(); + + let error = Infinity; + semiconvergents.forEach(semiconvergent => { + const newError = Math.abs(Math.PI - semiconvergent.valueOf()); + expect(newError).toBeLessThan(error); + error = newError; + }); + }); + + it('calculates the non-monotonic semiconvergents of pi', () => { + const semiconvergents = getConvergents(Math.PI, undefined, 5, true, true); + expect(semiconvergents).toHaveLength(5); + expect(semiconvergents[0].equals(3)).toBeTruthy(); + expect(semiconvergents[1].equals(4)).toBeTruthy(); + expect(semiconvergents[2].equals('7/2')).toBeTruthy(); + expect(semiconvergents[3].equals('10/3')).toBeTruthy(); + expect(semiconvergents[4].equals('13/4')).toBeTruthy(); + }); +}); diff --git a/src/__tests__/index.spec.ts b/src/__tests__/index.spec.ts index 6c8f02c..005b406 100644 --- a/src/__tests__/index.spec.ts +++ b/src/__tests__/index.spec.ts @@ -1,8 +1,5 @@ import {describe, it, expect} from 'vitest'; import { - approximateOddLimit, - approximatePrimeLimit, - approximatePrimeLimitWithErrors, arraysEqual, binomial, circleDifference, @@ -12,10 +9,8 @@ import { dot, extendedEuclid, gcd, - getConvergents, iteratedEuclid, norm, - PRIMES, valueToCents, } from '../index'; @@ -78,113 +73,6 @@ describe('binomial coefficient', () => { }); }); -describe('Odd limit approximator', () => { - it('can approximate tau in the 15-odd-limit', () => { - const approximation = approximateOddLimit(valueToCents(2 * Math.PI), 15)[0]; - expect(approximation.equals('44/7')).toBeTruthy(); - expect(approximation.valueOf()).toBeCloseTo(2 * Math.PI); - }); - - it('can approximate e in the 21-odd-limit', () => { - const approximations = approximateOddLimit(valueToCents(Math.E), 21); - expect(approximations[0].equals('19/7')).toBeTruthy(); - expect(approximations[0].valueOf()).toBeCloseTo(Math.E); - expect(approximations[7].equals('21/8')).toBeTruthy(); - expect(approximations[7].valueOf()).toBeCloseTo(Math.E, 0); - }); -}); - -describe('Prime limit approximator', () => { - it('can approximate pi in the 11-limit', () => { - const approximation = approximatePrimeLimit( - valueToCents(Math.PI), - PRIMES.indexOf(11), - 3 - )[0]; - expect(approximation.equals('12544/3993')).toBeTruthy(); - expect(approximation.valueOf()).toBeCloseTo(Math.PI); - }); - - it('can approximate pi in the 13-limit with a small sized result', () => { - const approximations = approximatePrimeLimit( - valueToCents(Math.PI), - PRIMES.indexOf(13), - 3, - 15, - 4 - ); - expect(approximations).toHaveLength(4); - }); - - it('can approximate the square root of two in the 7-limit within maximum error', () => { - const approximationsAndErrors = approximatePrimeLimitWithErrors( - 600, - PRIMES.indexOf(7), - 5, - 10 - ); - expect(approximationsAndErrors).toHaveLength(28); - approximationsAndErrors.forEach(([approximation, error]) => { - const cents = valueToCents(approximation.valueOf()); - const calculatedError = Math.abs(cents - 600); - expect(error).toBeCloseTo(calculatedError); - expect(calculatedError).toBeLessThanOrEqual(10); - }); - }); -}); - -describe('Convergent calculator', () => { - it('calculates the convergents of pi', () => { - const convergents = getConvergents(Math.PI, undefined, 10); - expect(convergents).toHaveLength(10); - expect(convergents[0].equals(3)).toBeTruthy(); - expect(convergents[1].equals('22/7')).toBeTruthy(); - expect(convergents[2].equals('333/106')).toBeTruthy(); - expect(convergents[3].equals('355/113')).toBeTruthy(); - expect(convergents[4].equals('103993/33102')).toBeTruthy(); - expect(convergents[5].equals('104348/33215')).toBeTruthy(); - expect(convergents[6].equals('208341/66317')).toBeTruthy(); - expect(convergents[7].equals('312689/99532')).toBeTruthy(); - expect(convergents[8].equals('833719/265381')).toBeTruthy(); - expect(convergents[9].equals('1146408/364913')).toBeTruthy(); - }); - - it('calculates the semiconvergents of pi', () => { - const semiconvergents = getConvergents(Math.PI, undefined, 13, true); - expect(semiconvergents).toHaveLength(13); - expect(semiconvergents[0].equals(3)).toBeTruthy(); - expect(semiconvergents[1].equals('13/4')).toBeTruthy(); - expect(semiconvergents[2].equals('16/5')).toBeTruthy(); - expect(semiconvergents[3].equals('19/6')).toBeTruthy(); - expect(semiconvergents[4].equals('22/7')).toBeTruthy(); - expect(semiconvergents[5].equals('179/57')).toBeTruthy(); - expect(semiconvergents[6].equals('201/64')).toBeTruthy(); - expect(semiconvergents[7].equals('223/71')).toBeTruthy(); - expect(semiconvergents[8].equals('245/78')).toBeTruthy(); - expect(semiconvergents[9].equals('267/85')).toBeTruthy(); - expect(semiconvergents[10].equals('289/92')).toBeTruthy(); - expect(semiconvergents[11].equals('311/99')).toBeTruthy(); - expect(semiconvergents[12].equals('333/106')).toBeTruthy(); - - let error = Infinity; - semiconvergents.forEach(semiconvergent => { - const newError = Math.abs(Math.PI - semiconvergent.valueOf()); - expect(newError).toBeLessThan(error); - error = newError; - }); - }); - - it('calculates the non-monotonic semiconvergents of pi', () => { - const semiconvergents = getConvergents(Math.PI, undefined, 5, true, true); - expect(semiconvergents).toHaveLength(5); - expect(semiconvergents[0].equals(3)).toBeTruthy(); - expect(semiconvergents[1].equals(4)).toBeTruthy(); - expect(semiconvergents[2].equals('7/2')).toBeTruthy(); - expect(semiconvergents[3].equals('10/3')).toBeTruthy(); - expect(semiconvergents[4].equals('13/4')).toBeTruthy(); - }); -}); - describe('Dot product', () => { it('can be used with all number arrays', () => { const a = new Float32Array([1, 2, 3, 4]); diff --git a/src/approximation.ts b/src/approximation.ts new file mode 100644 index 0000000..fea696c --- /dev/null +++ b/src/approximation.ts @@ -0,0 +1,271 @@ +import {valueToCents} from './conversion'; +import {Fraction, FractionValue, mmod} from './fraction'; +import {PRIMES, PRIME_CENTS} from './primes'; + +/** + * Calculate best rational approximations to a given fraction that are + * closer than any approximation with a smaller or equal denominator + * unless non-monotonic approximations are requested as well. + * @param value The fraction to simplify. + * @param maxDenominator Maximum denominator to include. + * @param maxLength Maximum length of the array of approximations. + * @param includeSemiconvergents Include semiconvergents. + * @param includeNonMonotonic Include non-monotonically improving approximations. + * @returns An array of (semi)convergents. + */ +export function getConvergents( + value: FractionValue, + maxDenominator?: number, + maxLength?: number, + includeSemiconvergents = false, + includeNonMonotonic = false +) { + const value_ = new Fraction(value); + /* + Glossary + cfDigit : the continued fraction digit + num : the convergent numerator + den : the convergent denominator + scnum : the semiconvergent numerator + scden : the semiconvergen denominator + cind : tracks indicies of convergents + */ + const result: Fraction[] = []; + const cf = value_.toContinued(); + const cind: number[] = []; + for (let d = 0; d < cf.length; d++) { + const cfDigit = cf[d]; + let num = cfDigit; + let den = 1; + // Calculate the convergent. + for (let i = d; i > 0; i--) { + [den, num] = [num, den]; + num += den * cf[i - 1]; + } + if (includeSemiconvergents && d > 0) { + const lowerBound = includeNonMonotonic ? 1 : Math.ceil(cfDigit / 2); + for (let i = lowerBound; i < cfDigit; i++) { + const scnum = num - (cfDigit - i) * result[cind[d - 1]].n; + const scden = den - (cfDigit - i) * result[cind[d - 1]].d; + if (scden > maxDenominator!) break; + const convergent = new Fraction(scnum, scden); + if (includeNonMonotonic) { + result.push(convergent); + } else { + // See https://en.wikipedia.org/wiki/Continued_fraction#Semiconvergents + // for the origin of this half-rule + if (2 * i > cfDigit) { + result.push(convergent); + } else if ( + convergent + .sub(value_) + .abs() + .compare(result[result.length - 1].sub(value_).abs()) < 0 + ) { + result.push(convergent); + } + } + if (result.length >= maxLength!) { + return result; + } + } + } + if (den > maxDenominator!) break; + cind.push(result.length); + result.push(new Fraction(num, den)); + if (result.length >= maxLength!) { + return result; + } + } + return result; +} + +// Cache of odd limit fractions. Expanded as necessary. +const ODD_FRACTIONS = [new Fraction(1), new Fraction(1, 3), new Fraction(3)]; +const ODD_CENTS = [0, -PRIME_CENTS[1], PRIME_CENTS[1]]; +const ODD_BREAKPOINTS = [1, 3]; +const TWO = new Fraction(2); + +/** + * Approximate a musical interval by ratios of which neither the numerator or denominator + * exceeds a specified limit, once all powers of 2 are removed. + * @param cents Size of the musical interval measured in cents. + * @param limit Maximum odd limit. + * @returns All odd limit fractions within 600 cents of the input value sorted by closeness with cent offsets attached. + */ +export function approximateOddLimitWithErrors(cents: number, limit: number) { + const breakpointIndex = (limit - 1) / 2; + // Expand cache. + while (ODD_BREAKPOINTS.length <= breakpointIndex) { + const newLimit = ODD_BREAKPOINTS.length * 2 + 1; + for (let numerator = 1; numerator <= newLimit; numerator += 2) { + for (let denominator = 1; denominator <= newLimit; denominator += 2) { + const fraction = new Fraction(numerator, denominator); + let novel = true; + for (let i = 0; i < ODD_FRACTIONS.length; ++i) { + if (fraction.equals(ODD_FRACTIONS[i])) { + novel = false; + break; + } + } + if (novel) { + ODD_FRACTIONS.push(fraction); + ODD_CENTS.push(valueToCents(fraction.valueOf())); + } + } + } + ODD_BREAKPOINTS.push(ODD_FRACTIONS.length); + } + + // Find closest odd limit fractions modulo octaves. + const results: [Fraction, number][] = []; + for (let i = 0; i < ODD_BREAKPOINTS[breakpointIndex]; ++i) { + const oddCents = ODD_CENTS[i]; + const remainder = mmod(cents - oddCents, 1200); + // Overshot + if (remainder <= 600) { + // Rounding done to eliminate floating point jitter. + const exponent = Math.round((cents - oddCents - remainder) / 1200); + const error = remainder; + // Exponentiate to add the required number of octaves. + results.push([ODD_FRACTIONS[i].mul(TWO.pow(exponent)!), error]); + } + // Undershot + else { + const exponent = Math.round((cents - oddCents - remainder) / 1200) + 1; + const error = 1200 - remainder; + results.push([ODD_FRACTIONS[i].mul(TWO.pow(exponent)!), error]); + } + } + + results.sort((a, b) => a[1] - b[1]); + + return results; +} + +/** + * Approximate a musical interval by ratios of which neither the numerator or denominator + * exceeds a specified limit, once all powers of 2 are removed. + * @param cents Size of the musical interval measured in cents. + * @param limit Maximum odd limit. + * @returns All odd limit fractions within 600 cents of the input value sorted by closeness. + */ +export function approximateOddLimit(cents: number, limit: number) { + return approximateOddLimitWithErrors(cents, limit).map(result => result[0]); +} + +/** + * Approximate a musical interval by ratios of which are within a prime limit with + * exponents that do not exceed the maximimum, exponent of 2 ignored. + * @param cents Size of the musical interval measured in cents. + * @param limitIndex The ordinal of the prime of the limit. + * @param maxError Maximum error from the interval for inclusion in the result. + * @param maxLength Maximum number of approximations to return. + * @returns All valid fractions within `maxError` cents of the input value sorted by closeness with cent offsets attached. + */ +export function approximatePrimeLimitWithErrors( + cents: number, + limitIndex: number, + maxExponent: number, + maxError = 600, + maxLength?: number +) { + if (maxError > 600) { + throw new Error('Maximum search distance is 600 cents'); + } + const results: [Fraction, number][] = []; + + function push(error: number, generator: () => Fraction) { + if (maxLength === undefined || results.length < maxLength) { + results.push([generator(), error]); + if (results.length === maxLength) { + results.sort((a, b) => a[1] - b[1]); + } + } else { + if (error > results[results.length - 1][1]) { + return; + } + for (let index = results.length - 1; ; index--) { + if (error <= results[index][1]) { + results.splice(index, 0, [generator(), error]); + break; + } + } + results.pop(); + } + } + + function accumulate( + approximation: Fraction, + approximationCents: number, + index: number + ) { + if (index > limitIndex) { + // Procedure is the same as in approximateOddLimitWithErrors + const remainder = mmod(cents - approximationCents, 1200); + if (remainder <= 600) { + const error = remainder; + if (error > maxError) { + return; + } + push(error, () => + approximation.mul( + TWO.pow( + Math.round((cents - approximationCents - remainder) / 1200) + )! + ) + ); + } else { + const error = 1200 - remainder; + if (error > maxError) { + return; + } + push(error, () => + approximation.mul( + TWO.pow( + Math.round((cents - approximationCents - remainder) / 1200) + 1 + )! + ) + ); + } + return; + } + approximation = approximation.div(PRIMES[index] ** maxExponent); + approximationCents -= PRIME_CENTS[index] * maxExponent; + for (let i = -maxExponent; i <= maxExponent; ++i) { + accumulate(approximation, approximationCents, index + 1); + approximation = approximation.mul(PRIMES[index]); + approximationCents += PRIME_CENTS[index]; + } + } + accumulate(new Fraction(1), 0, 1); + + results.sort((a, b) => a[1] - b[1]); + + return results; +} + +/** + * Approximate a musical interval by ratios of which are within a prime limit with + * exponents that do not exceed the maximimum, exponent of 2 ignored. + * @param cents Size of the musical interval measured in cents. + * @param limitIndex The ordinal of the prime of the limit. + * @param maxError Maximum error from the interval for inclusion in the result. + * @param maxLength Maximum number of approximations to return. + * @returns All valid fractions within `maxError` cents of the input value sorted by closenesss. + */ +export function approximatePrimeLimit( + cents: number, + limitIndex: number, + maxExponent: number, + maxError = 600, + maxLength?: number +) { + return approximatePrimeLimitWithErrors( + cents, + limitIndex, + maxExponent, + maxError, + maxLength + ).map(result => result[0]); +} diff --git a/src/index.ts b/src/index.ts index a9eaabc..794f255 100644 --- a/src/index.ts +++ b/src/index.ts @@ -1,12 +1,11 @@ -import {valueToCents} from './conversion'; -import {Fraction, FractionValue, mmod} from './fraction'; -import {PRIMES, PRIME_CENTS} from './primes'; +import {Fraction, mmod} from './fraction'; export * from './fraction'; export * from './primes'; export * from './conversion'; export * from './combinations'; export * from './monzo'; +export * from './approximation'; export interface AnyArray { [key: number]: any; @@ -125,84 +124,6 @@ export function iteratedEuclid(params: Iterable) { return coefs; } -/** - * Calculate best rational approximations to a given fraction that are - * closer than any approximation with a smaller or equal denominator - * unless non-monotonic approximations are requested as well. - * @param value The fraction to simplify. - * @param maxDenominator Maximum denominator to include. - * @param maxLength Maximum length of the array of approximations. - * @param includeSemiconvergents Include semiconvergents. - * @param includeNonMonotonic Include non-monotonically improving approximations. - * @returns An array of (semi)convergents. - */ -export function getConvergents( - value: FractionValue, - maxDenominator?: number, - maxLength?: number, - includeSemiconvergents = false, - includeNonMonotonic = false -) { - const value_ = new Fraction(value); - /* - Glossary - cfDigit : the continued fraction digit - num : the convergent numerator - den : the convergent denominator - scnum : the semiconvergent numerator - scden : the semiconvergen denominator - cind : tracks indicies of convergents - */ - const result: Fraction[] = []; - const cf = value_.toContinued(); - const cind: number[] = []; - for (let d = 0; d < cf.length; d++) { - const cfDigit = cf[d]; - let num = cfDigit; - let den = 1; - // Calculate the convergent. - for (let i = d; i > 0; i--) { - [den, num] = [num, den]; - num += den * cf[i - 1]; - } - if (includeSemiconvergents && d > 0) { - const lowerBound = includeNonMonotonic ? 1 : Math.ceil(cfDigit / 2); - for (let i = lowerBound; i < cfDigit; i++) { - const scnum = num - (cfDigit - i) * result[cind[d - 1]].n; - const scden = den - (cfDigit - i) * result[cind[d - 1]].d; - if (scden > maxDenominator!) break; - const convergent = new Fraction(scnum, scden); - if (includeNonMonotonic) { - result.push(convergent); - } else { - // See https://en.wikipedia.org/wiki/Continued_fraction#Semiconvergents - // for the origin of this half-rule - if (2 * i > cfDigit) { - result.push(convergent); - } else if ( - convergent - .sub(value_) - .abs() - .compare(result[result.length - 1].sub(value_).abs()) < 0 - ) { - result.push(convergent); - } - } - if (result.length >= maxLength!) { - return result; - } - } - } - if (den > maxDenominator!) break; - cind.push(result.length); - result.push(new Fraction(num, den)); - if (result.length >= maxLength!) { - return result; - } - } - return result; -} - /** * Collection of unique fractions. */ @@ -301,196 +222,6 @@ export function clamp(minValue: number, maxValue: number, value: number) { return value; } -// Cache of odd limit fractions. Expanded as necessary. -const ODD_FRACTIONS = [new Fraction(1), new Fraction(1, 3), new Fraction(3)]; -const ODD_CENTS = [0, -PRIME_CENTS[1], PRIME_CENTS[1]]; -const ODD_BREAKPOINTS = [1, 3]; -const TWO = new Fraction(2); - -/** - * Approximate a musical interval by ratios of which neither the numerator or denominator - * exceeds a specified limit, once all powers of 2 are removed. - * @param cents Size of the musical interval measured in cents. - * @param limit Maximum odd limit. - * @returns All odd limit fractions within 600 cents of the input value sorted by closeness with cent offsets attached. - */ -export function approximateOddLimitWithErrors(cents: number, limit: number) { - const breakpointIndex = (limit - 1) / 2; - // Expand cache. - while (ODD_BREAKPOINTS.length <= breakpointIndex) { - const newLimit = ODD_BREAKPOINTS.length * 2 + 1; - for (let numerator = 1; numerator <= newLimit; numerator += 2) { - for (let denominator = 1; denominator <= newLimit; denominator += 2) { - const fraction = new Fraction(numerator, denominator); - let novel = true; - for (let i = 0; i < ODD_FRACTIONS.length; ++i) { - if (fraction.equals(ODD_FRACTIONS[i])) { - novel = false; - break; - } - } - if (novel) { - ODD_FRACTIONS.push(fraction); - ODD_CENTS.push(valueToCents(fraction.valueOf())); - } - } - } - ODD_BREAKPOINTS.push(ODD_FRACTIONS.length); - } - - // Find closest odd limit fractions modulo octaves. - const results: [Fraction, number][] = []; - for (let i = 0; i < ODD_BREAKPOINTS[breakpointIndex]; ++i) { - const oddCents = ODD_CENTS[i]; - const remainder = mmod(cents - oddCents, 1200); - // Overshot - if (remainder <= 600) { - // Rounding done to eliminate floating point jitter. - const exponent = Math.round((cents - oddCents - remainder) / 1200); - const error = remainder; - // Exponentiate to add the required number of octaves. - results.push([ODD_FRACTIONS[i].mul(TWO.pow(exponent)!), error]); - } - // Undershot - else { - const exponent = Math.round((cents - oddCents - remainder) / 1200) + 1; - const error = 1200 - remainder; - results.push([ODD_FRACTIONS[i].mul(TWO.pow(exponent)!), error]); - } - } - - results.sort((a, b) => a[1] - b[1]); - - return results; -} - -/** - * Approximate a musical interval by ratios of which neither the numerator or denominator - * exceeds a specified limit, once all powers of 2 are removed. - * @param cents Size of the musical interval measured in cents. - * @param limit Maximum odd limit. - * @returns All odd limit fractions within 600 cents of the input value sorted by closeness. - */ -export function approximateOddLimit(cents: number, limit: number) { - return approximateOddLimitWithErrors(cents, limit).map(result => result[0]); -} - -/** - * Approximate a musical interval by ratios of which are within a prime limit with - * exponents that do not exceed the maximimum, exponent of 2 ignored. - * @param cents Size of the musical interval measured in cents. - * @param limitIndex The ordinal of the prime of the limit. - * @param maxError Maximum error from the interval for inclusion in the result. - * @param maxLength Maximum number of approximations to return. - * @returns All valid fractions within `maxError` cents of the input value sorted by closeness with cent offsets attached. - */ -export function approximatePrimeLimitWithErrors( - cents: number, - limitIndex: number, - maxExponent: number, - maxError = 600, - maxLength?: number -) { - if (maxError > 600) { - throw new Error('Maximum search distance is 600 cents'); - } - const results: [Fraction, number][] = []; - - function push(error: number, generator: () => Fraction) { - if (maxLength === undefined || results.length < maxLength) { - results.push([generator(), error]); - if (results.length === maxLength) { - results.sort((a, b) => a[1] - b[1]); - } - } else { - if (error > results[results.length - 1][1]) { - return; - } - for (let index = results.length - 1; ; index--) { - if (error <= results[index][1]) { - results.splice(index, 0, [generator(), error]); - break; - } - } - results.pop(); - } - } - - function accumulate( - approximation: Fraction, - approximationCents: number, - index: number - ) { - if (index > limitIndex) { - // Procedure is the same as in approximateOddLimitWithErrors - const remainder = mmod(cents - approximationCents, 1200); - if (remainder <= 600) { - const error = remainder; - if (error > maxError) { - return; - } - push(error, () => - approximation.mul( - TWO.pow( - Math.round((cents - approximationCents - remainder) / 1200) - )! - ) - ); - } else { - const error = 1200 - remainder; - if (error > maxError) { - return; - } - push(error, () => - approximation.mul( - TWO.pow( - Math.round((cents - approximationCents - remainder) / 1200) + 1 - )! - ) - ); - } - return; - } - approximation = approximation.div(PRIMES[index] ** maxExponent); - approximationCents -= PRIME_CENTS[index] * maxExponent; - for (let i = -maxExponent; i <= maxExponent; ++i) { - accumulate(approximation, approximationCents, index + 1); - approximation = approximation.mul(PRIMES[index]); - approximationCents += PRIME_CENTS[index]; - } - } - accumulate(new Fraction(1), 0, 1); - - results.sort((a, b) => a[1] - b[1]); - - return results; -} - -/** - * Approximate a musical interval by ratios of which are within a prime limit with - * exponents that do not exceed the maximimum, exponent of 2 ignored. - * @param cents Size of the musical interval measured in cents. - * @param limitIndex The ordinal of the prime of the limit. - * @param maxError Maximum error from the interval for inclusion in the result. - * @param maxLength Maximum number of approximations to return. - * @returns All valid fractions within `maxError` cents of the input value sorted by closenesss. - */ -export function approximatePrimeLimit( - cents: number, - limitIndex: number, - maxExponent: number, - maxError = 600, - maxLength?: number -) { - return approximatePrimeLimitWithErrors( - cents, - limitIndex, - maxExponent, - maxError, - maxLength - ).map(result => result[0]); -} - /** * Calculate the inner (dot) product of two arrays of real numbers. * @param a The first array of numbers.