Implement a method for accurately measuring the sizes of small commas

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
+    );
+  });
+});

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;

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;

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);
+}

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<number>) {
+  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;
+}