From f6e03fc6c86268df132cfd15af884acf01090feb Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Wed, 24 Apr 2024 20:46:44 +0300 Subject: [PATCH] Implement a method for accurately measuring the sizes of small commas --- src/__tests__/monzo.spec.ts | 22 ++++ src/commas.ts | 35 +++++++ src/index.ts | 1 + src/monzo.ts | 96 ++++++++++++++++++ src/polyfils/sum-precise.ts | 197 ++++++++++++++++++++++++++++++++++++ 5 files changed, 351 insertions(+) create mode 100644 src/commas.ts create mode 100644 src/polyfils/sum-precise.ts diff --git a/src/__tests__/monzo.spec.ts b/src/__tests__/monzo.spec.ts index a6bad7b..d0e6674 100644 --- a/src/__tests__/monzo.spec.ts +++ b/src/__tests__/monzo.spec.ts @@ -2,6 +2,7 @@ import {describe, it, expect} from 'vitest'; import {Fraction} from '../fraction'; import { monzoToBigInt, + monzoToCents, monzoToFraction, primeLimit, toMonzo, @@ -372,3 +373,24 @@ describe('Prime limit calculator', () => { expect(primeLimit(2n ** 1025n)).toEqual(2); }); }); + +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 neutrino accurately', () => { + // XXX: Produces 1.6359101573382162e-10. Three digits of accuracy, smh... + expect(monzoToCents([1889, -2145, 138, 424])).toBeCloseTo( + 1.6361187484440885e-10, + 13 + ); + }); + + it('calculates the size of the demiquartervice comma accurately', () => { + expect(monzoToCents([-3, 2, -1, -1, 0, 0, -1, 0, 2])).toBeCloseTo( + 0.3636664332386927, + 14 + ); + }); +}); diff --git a/src/commas.ts b/src/commas.ts new file mode 100644 index 0000000..fe25ad2 --- /dev/null +++ b/src/commas.ts @@ -0,0 +1,35 @@ +// Mapping commas for measuring the sizes of small intervals + +export const MAPPING_COMMA_MONZOS = [ + [1], // Octave + [-1, 1], // Pure fifth + [554, -351, 1], // Quectisma + [-55, 30, 2, 1], // Nommisma + [-30, 27, -7, 0, 1], // Negative syntonoschisma / syntonisma + [9, 0, -1, 0, -3, 1], // Jacobin comma + [-2, 2, -1, -5, 0, 3, 1], // Aksial comma + [9, -3, -3, -2, 0, 0, 1, 1], // Decimillisma + [-1, -1, 0, -1, -2, 1, 1, 0, 1], // Broadviewsma +]; + +export const MAPPING_COMMA_CENTS = [ + 1200, 701.9550008653874, 0.10841011385118912, 0.1033601604170961, + -0.09303132362673557, 0.2601208102056527, 0.005150328654440726, + 0.010468503793319029, 0.34062647410756758, +]; + +// Auxiliary commas for chipping away at the 3-limit +export const SATANIC_COMMA_MONZO = [-1054, 665]; +export const SATANIC_COMMA_CENTS = 0.07557548263280008; + +export const MERCATOR_COMMA_MONZO = [-84, 53]; +export const MERCATOR_COMMA_CENTS = 3.61504586553314; + +export const PYTH_COMMA_MONZO = [-19, 12]; +export const PYTH_COMMA_CENTS = 23.46001038464901; + +export const PYTH_LIMMA_MONZO = [8, -5]; +export const PYTH_LIMMA_CENTS = 90.22499567306291; + +export const PYTH_TONE_MONZO = [-3, 2]; +export const PYTH_TONE_CENTS = 203.9100017307748; diff --git a/src/index.ts b/src/index.ts index 2e914f8..0100b6e 100644 --- a/src/index.ts +++ b/src/index.ts @@ -6,6 +6,7 @@ export * from './conversion'; export * from './combinations'; export * from './monzo'; export * from './approximation'; +export {sum} from './polyfils/sum-precise'; export interface AnyArray { [key: number]: any; diff --git a/src/monzo.ts b/src/monzo.ts index 3880a5d..118cfcc 100644 --- a/src/monzo.ts +++ b/src/monzo.ts @@ -1,5 +1,20 @@ import {Fraction, FractionValue} from './fraction'; import {BIG_INT_PRIMES, PRIMES} from './primes'; +import { + MAPPING_COMMA_CENTS, + MAPPING_COMMA_MONZOS, + MERCATOR_COMMA_CENTS, + MERCATOR_COMMA_MONZO, + PYTH_COMMA_CENTS, + PYTH_COMMA_MONZO, + PYTH_LIMMA_CENTS, + PYTH_LIMMA_MONZO, + PYTH_TONE_CENTS, + PYTH_TONE_MONZO, + SATANIC_COMMA_CENTS, + SATANIC_COMMA_MONZO, +} from './commas'; +import {sum} from './polyfils/sum-precise'; /** * Array of integers representing the exponents of prime numbers in the unique factorization of a rational number. @@ -593,3 +608,84 @@ function bigIntToMonzo7(n: bigint): [Monzo, bigint] { } return [result, n]; } + +export function monzoToCents(monzo: Monzo) { + monzo = [...monzo]; + while (monzo.length && !monzo[monzo.length - 1]) { + monzo.pop(); + } + const terms = []; + while (monzo.length > 2) { + const component = monzo.pop()!; + terms.push(component * MAPPING_COMMA_CENTS[monzo.length]); + const comma = MAPPING_COMMA_MONZOS[monzo.length]; + for (let i = 0; i < monzo.length; ++i) { + monzo[i] -= component * comma[i]; + } + } + if (monzo.length === 2) { + while (monzo[1] >= 665) { + terms.push(SATANIC_COMMA_CENTS); + monzo[0] -= SATANIC_COMMA_MONZO[0]; + monzo[1] -= SATANIC_COMMA_MONZO[1]; + } + while (monzo[1] <= -665) { + terms.push(-SATANIC_COMMA_CENTS); + monzo[0] += SATANIC_COMMA_MONZO[0]; + monzo[1] += SATANIC_COMMA_MONZO[1]; + } + + while (monzo[1] >= 53) { + terms.push(MERCATOR_COMMA_CENTS); + monzo[0] -= MERCATOR_COMMA_MONZO[0]; + monzo[1] -= MERCATOR_COMMA_MONZO[1]; + } + while (monzo[1] <= -53) { + terms.push(-MERCATOR_COMMA_CENTS); + monzo[0] += MERCATOR_COMMA_MONZO[0]; + monzo[1] += MERCATOR_COMMA_MONZO[1]; + } + + while (monzo[1] >= 12) { + terms.push(PYTH_COMMA_CENTS); + monzo[0] -= PYTH_COMMA_MONZO[0]; + monzo[1] -= PYTH_COMMA_MONZO[1]; + } + while (monzo[1] <= -12) { + terms.push(-PYTH_COMMA_CENTS); + monzo[0] += PYTH_COMMA_MONZO[0]; + monzo[1] += PYTH_COMMA_MONZO[1]; + } + + while (monzo[1] >= 5) { + terms.push(-PYTH_LIMMA_CENTS); + monzo[0] += PYTH_LIMMA_MONZO[0]; + monzo[1] += PYTH_LIMMA_MONZO[1]; + } + while (monzo[1] <= -5) { + terms.push(PYTH_LIMMA_CENTS); + monzo[0] -= PYTH_LIMMA_MONZO[0]; + monzo[1] -= PYTH_LIMMA_MONZO[1]; + } + + while (monzo[1] >= 2) { + terms.push(PYTH_TONE_CENTS); + monzo[0] -= PYTH_TONE_MONZO[0]; + monzo[1] -= PYTH_TONE_MONZO[1]; + } + while (monzo[1] <= -2) { + terms.push(-PYTH_TONE_CENTS); + monzo[0] += PYTH_TONE_MONZO[0]; + monzo[1] += PYTH_TONE_MONZO[1]; + } + } + while (monzo.length) { + const component = monzo.pop()!; + terms.push(component * MAPPING_COMMA_CENTS[monzo.length]); + const comma = MAPPING_COMMA_MONZOS[monzo.length]; + for (let i = 0; i < monzo.length; ++i) { + monzo[i] -= component * comma[i]; + } + } + return sum(terms); +} diff --git a/src/polyfils/sum-precise.ts b/src/polyfils/sum-precise.ts new file mode 100644 index 0000000..d799b1e --- /dev/null +++ b/src/polyfils/sum-precise.ts @@ -0,0 +1,197 @@ +// 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}; +} + +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; +}