Implement a method for accurately measuring the sizes of small commas

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

src/index.ts

@@ -1,11 +1,15 @@
 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';
 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;
+}

src/monzo.ts

@@ -278,6 +278,27 @@ export function monzoToBigInt(monzo: Iterable<number>) {
   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<number>) {
+  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.

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