From c133dd135b0057a3cad0043ddb3bbf3a07d786d0 Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Thu, 25 Apr 2024 08:43:46 +0300 Subject: [PATCH] Implement a method for accurately measuring the sizes of small commas Add a polyfill for Array.sumPrecise. Implement monzoToBigNumeratorDenominator. Implement dotPrecise. --- src/__tests__/index.spec.ts | 35 ++++++ src/index.ts | 53 ++++++++- src/monzo.ts | 21 ++++ src/polyfills/sum-precise.ts | 202 +++++++++++++++++++++++++++++++++++ 4 files changed, 310 insertions(+), 1 deletion(-) create mode 100644 src/polyfills/sum-precise.ts diff --git a/src/__tests__/index.spec.ts b/src/__tests__/index.spec.ts index 1ebb389..07d3614 100644 --- a/src/__tests__/index.spec.ts +++ b/src/__tests__/index.spec.ts @@ -19,6 +19,7 @@ import { iteratedEuclid, norm, valueToCents, + monzoToCents, } from '../index'; describe('Array equality tester', () => { @@ -382,3 +383,37 @@ describe('Constant structure checker with a margin of equivalence', () => { expect(hasMarginConstantStructure([1199, 1200], 2)).toBe(false); }); }); + +describe('Monzo size measure', () => { + it('calculates the size of the perfect fourth accurately', () => { + expect(monzoToCents([2, -1])).toBeCloseTo(498.0449991346125, 12); + }); + + it('calculates the size of the rascal accurately', () => { + expect(monzoToCents([-7470, 2791, 1312])).toBeCloseTo( + 5.959563411893381e-6, + 24 + ); + }); + + it('calculates the size of the neutrino accurately', () => { + expect(monzoToCents([1889, -2145, 138, 424])).toBeCloseTo( + 1.6361187484440885e-10, + 24 + ); + }); + + it('calculates the size of the demiquartervice comma accurately', () => { + expect(monzoToCents([-3, 2, -1, -1, 0, 0, -1, 0, 2])).toBeCloseTo( + 0.3636664332386927, + 15 + ); + }); + + it('calculates the size of the negative junebug comma accurately', () => { + expect(monzoToCents([-1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1])).toBeCloseTo( + -6.104006661651758, + 15 + ); + }); +}); diff --git a/src/index.ts b/src/index.ts index 2e914f8..5d22e48 100644 --- a/src/index.ts +++ b/src/index.ts @@ -1,4 +1,7 @@ import {Fraction, mmod} from './fraction'; +import {Monzo, monzoToBigNumeratorDenominator} from './monzo'; +import {PRIME_CENTS} from './primes'; +import {sum} from './polyfills/sum-precise'; export * from './fraction'; export * from './primes'; @@ -6,6 +9,7 @@ export * from './conversion'; export * from './combinations'; export * from './monzo'; export * from './approximation'; +export {sum} from './polyfills/sum-precise'; export interface AnyArray { [key: number]: any; @@ -229,13 +233,31 @@ export function clamp(minValue: number, maxValue: number, value: number) { * @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 < Math.min(a.length, b.length); ++i) { + 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. @@ -449,3 +471,32 @@ export function hasMarginConstantStructure( } return true; } + +const NATS_TO_CENTS = 1200 / Math.LN2; +const IEEE_LIMIT = 2n ** 1024n; + +/** + * Measure the size of a monzo in cents. + * Monzos representing small rational numbers (commas) are measured accurately. + * @param monzo Array or prime exponents, possibly fractional. + * @returns The size of the represented number in cents (1200ths of an octave). + */ +export function monzoToCents(monzo: Monzo) { + const result = dotPrecise(monzo, PRIME_CENTS); + if (Math.abs(result) > 10) { + return result; + } + for (const component of monzo) { + if (!Number.isInteger(component)) { + return result; + } + } + let {numerator, denominator} = monzoToBigNumeratorDenominator(monzo); + let delta = numerator - denominator; + // The answer is smaller than 10 cents so no need to check delta here or worry about its sign + while (denominator >= IEEE_LIMIT) { + delta >>= 1n; + denominator >>= 1n; + } + return Math.log1p(Number(delta) / Number(denominator)) * NATS_TO_CENTS; +} diff --git a/src/monzo.ts b/src/monzo.ts index 3880a5d..8a298f2 100644 --- a/src/monzo.ts +++ b/src/monzo.ts @@ -278,6 +278,27 @@ export function monzoToBigInt(monzo: Iterable) { return result; } +/** + * Convert a monzo to the BigInt fraction it represents. + * @param monzo Iterable of prime exponents. + * @returns Record with keys 'numerator' and 'denominator containing BigInts. + */ +export function monzoToBigNumeratorDenominator(monzo: Iterable) { + let numerator = 1n; + let denominator = 1n; + let index = 0; + for (const component of monzo) { + if (component > 0) { + numerator *= BIG_INT_PRIMES[index] ** BigInt(component); + } + if (component < 0) { + denominator *= BIG_INT_PRIMES[index] ** BigInt(-component); + } + index++; + } + return {numerator, denominator}; +} + /** * Calculate the prime limit of an integer or a fraction. * @param n Integer or fraction to calculate prime limit for. diff --git a/src/polyfills/sum-precise.ts b/src/polyfills/sum-precise.ts new file mode 100644 index 0000000..ef55507 --- /dev/null +++ b/src/polyfills/sum-precise.ts @@ -0,0 +1,202 @@ +// Stolen from: https://github.com/tc39/proposal-math-sum/blob/main/polyfill/polyfill.mjs +// Linted and type-annotated by Lumi Pakkanen. + +/* +https://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps +Shewchuk's algorithm for exactly floating point addition +as implemented in Python's fsum: https://github.com/python/cpython/blob/48dfd74a9db9d4aa9c6f23b4a67b461e5d977173/Modules/mathmodule.c#L1359-L1474 +adapted to handle overflow via an additional "biased" partial, representing 2**1024 times its actual value +*/ + +// exponent 11111111110, significand all 1s +const MAX_DOUBLE = 1.79769313486231570815e308; // i.e. (2**1024 - 2**(1023 - 52)) + +// exponent 11111111110, significand all 1s except final 0 +const PENULTIMATE_DOUBLE = 1.79769313486231550856e308; // i.e. (2**1024 - 2 * 2**(1023 - 52)) + +// exponent 11111001010, significand all 0s +const MAX_ULP = MAX_DOUBLE - PENULTIMATE_DOUBLE; // 1.99584030953471981166e+292, i.e. 2**(1023 - 52) + +// prerequisite: Math.abs(x) >= Math.abs(y) +function twosum(x: number, y: number) { + const hi = x + y; + const lo = y - (hi - x); + return {hi, lo}; +} + +/** + * Accurately add up elements from an iterable using Shewchuk's algorithm. + * @param iterable Numbers to sum together. + * @returns The sum of the elements. + */ +export function sum(iterable: Iterable) { + const partials: number[] = []; + + let overflow = 0; // conceptually 2**1024 times this value; the final partial + + // for purposes of the polyfill we're going to ignore closing the iterator, sorry + const iterator = iterable[Symbol.iterator](); + const next = iterator.next.bind(iterator); + + // in C this would be done using a goto + function drainNonFiniteValue(current: number) { + while (true) { + const {done, value} = next(); + if (done) { + return current; + } + if (!Number.isFinite(value)) { + // summing any distinct two of the three non-finite values gives NaN + // summing any one of them with itself gives itself + if (!Object.is(value, current)) { + current = NaN; + } + } + } + } + + // handle list of -0 special case + while (true) { + const {done, value} = next(); + if (done) { + return -0; + } + if (!Object.is(value, -0)) { + if (!Number.isFinite(value)) { + return drainNonFiniteValue(value); + } + partials.push(value); + break; + } + } + + // main loop + while (true) { + const {done, value} = next(); + if (done) { + break; + } + let x = +value; + if (!Number.isFinite(x)) { + return drainNonFiniteValue(x); + } + + // we're updating partials in place, but it is maybe easier to understand if you think of it as making a new copy + let actuallyUsedPartials = 0; + // let newPartials = []; + for (let y of partials) { + if (Math.abs(x) < Math.abs(y)) { + [x, y] = [y, x]; + } + let {hi, lo} = twosum(x, y); + if (Math.abs(hi) === Infinity) { + const sign = hi === Infinity ? 1 : -1; + overflow += sign; + if (Math.abs(overflow) >= 2 ** 53) { + throw new RangeError('overflow'); + } + + x = x - sign * 2 ** 1023 - sign * 2 ** 1023; + if (Math.abs(x) < Math.abs(y)) { + [x, y] = [y, x]; + } + ({hi, lo} = twosum(x, y)); + } + if (lo !== 0) { + partials[actuallyUsedPartials] = lo; + ++actuallyUsedPartials; + // newPartials.push(lo); + } + x = hi; + } + partials.length = actuallyUsedPartials; + // assert.deepStrictEqual(partials, newPartials) + // partials = newPartials + + if (x !== 0) { + partials.push(x); + } + } + + // compute the exact sum of partials, stopping once we lose precision + let n = partials.length - 1; + let hi = 0; + let lo = 0; + + if (overflow !== 0) { + const next = n >= 0 ? partials[n] : 0; + --n; + if ( + Math.abs(overflow) > 1 || + (overflow > 0 && next > 0) || + (overflow < 0 && next < 0) + ) { + return overflow > 0 ? Infinity : -Infinity; + } + // here we actually have to do the arithmetic + // drop a factor of 2 so we can do it without overflow + // assert(Math.abs(overflow) === 1) + ({hi, lo} = twosum(overflow * 2 ** 1023, next / 2)); + lo *= 2; + if (Math.abs(2 * hi) === Infinity) { + // stupid edge case: rounding to the maximum value + // MAX_DOUBLE has a 1 in the last place of its significand, so if we subtract exactly half a ULP from 2**1024, the result rounds away from it (i.e. to infinity) under ties-to-even + // but if the next partial has the opposite sign of the current value, we need to round towards MAX_DOUBLE instead + // this is the same as the "handle rounding" case below, but there's only one potentially-finite case we need to worry about, so we just hardcode that one + if (hi > 0) { + if ( + hi === 2 ** 1023 && + lo === -(MAX_ULP / 2) && + n >= 0 && + partials[n] < 0 + ) { + return MAX_DOUBLE; + } + return Infinity; + } else { + if ( + hi === -(2 ** 1023) && + lo === MAX_ULP / 2 && + n >= 0 && + partials[n] > 0 + ) { + return -MAX_DOUBLE; + } + return -Infinity; + } + } + if (lo !== 0) { + partials[n + 1] = lo; + ++n; + lo = 0; + } + hi *= 2; + } + + while (n >= 0) { + const x = hi; + const y = partials[n]; + --n; + // assert: Math.abs(x) > Math.abs(y) + ({hi, lo} = twosum(x, y)); + if (lo !== 0) { + break; + } + } + + // handle rounding + // when the roundoff error is exactly half of the ULP for the result, we need to check one more partial to know which way to round + if ( + n >= 0 && + ((lo < 0.0 && partials[n] < 0.0) || (lo > 0.0 && partials[n] > 0.0)) + ) { + const y = lo * 2.0; + const x = hi + y; + const yr = x - hi; + if (y === yr) { + hi = x; + } + } + + return hi; +}