Implement Tenney and Wilson complexity measures

Plug more holes in primeFactorize(). ref #34
xenharmonic-devs · Apr 30, 2024 · baf78ed · baf78ed
1 parent 2b16765
commit baf78ed
Show file tree

Hide file tree

Showing 4 changed files with 282 additions and 5 deletions.
diff --git a/src/__tests__/index.spec.ts b/src/__tests__/index.spec.ts
@@ -20,6 +20,8 @@ import {
   norm,
   valueToCents,
   monzoToCents,
+  tenneyHeight,
+  wilsonHeight,
 } from '../index';
 
 describe('Array equality tester', () => {
@@ -417,3 +419,203 @@ describe('Monzo size measure', () => {
     );
   });
 });
+
+describe('Tenney complexity measure', () => {
+  it('calculates the complexity of 88', () => {
+    expect(tenneyHeight(88)).toBeCloseTo(4.477);
+  });
+
+  it('calculates the complexity of 11/8', () => {
+    expect(tenneyHeight('11/8')).toBeCloseTo(4.477);
+  });
+
+  it('calculates the complexity of -11/8', () => {
+    expect(tenneyHeight('-11/8')).toBeCloseTo(4.477);
+  });
+
+  it('calculates the complexity of 8/11', () => {
+    expect(tenneyHeight([3, 0, 0, 0, -1])).toBeCloseTo(4.477);
+  });
+
+  it('has a value for zero', () => {
+    expect(tenneyHeight(0)).toBe(Infinity);
+  });
+
+  it.skip('fuzzes the fraction property', () => {
+    for (let i = 0; i < 100000; ++i) {
+      const n = Math.floor(Math.random() * Number.MAX_SAFE_INTEGER);
+      const d = Math.floor(Math.random() * (Number.MAX_SAFE_INTEGER - 1)) + 1;
+      const f = new Fraction(n, d);
+      expect(tenneyHeight(f.n) + tenneyHeight(f.d)).toBeCloseTo(
+        tenneyHeight(f)
+      );
+    }
+  });
+
+  it.skip('fuzzes the multiplicative property', () => {
+    for (let i = 0; i < 100000; ++i) {
+      const x = Math.floor(Math.random() * 94906265);
+      const y = Math.floor(Math.random() * 94906265);
+      expect(tenneyHeight(x) + tenneyHeight(y)).toBeCloseTo(
+        tenneyHeight(x * y)
+      );
+    }
+  });
+});
+
+describe('Wilson complexity measure', () => {
+  it('calculates the complexity of 88', () => {
+    expect(wilsonHeight(88)).toBe(17);
+  });
+
+  it('calculates the complexity of 11/8', () => {
+    expect(wilsonHeight('11/8')).toBe(17);
+  });
+
+  it('calculates the complexity of -11/8', () => {
+    expect(wilsonHeight('-11/8')).toBe(17);
+  });
+
+  it('calculates the complexity of 8/11', () => {
+    expect(wilsonHeight([3, 0, 0, 0, -1])).toBe(17);
+  });
+
+  it('calculates the complexity of -8/11', () => {
+    expect(
+      wilsonHeight(
+        new Map([
+          [-1, 1],
+          [2, 3],
+          [11, -1],
+        ])
+      )
+    ).toBe(17);
+  });
+
+  it('has a value for zero', () => {
+    expect(wilsonHeight(0)).toBe(Infinity);
+  });
+
+  it.each([
+    ['2/1', 2],
+    ['3/2', 5],
+    ['4/3', 7],
+    ['5/4', 9],
+    ['6/5', 10],
+    ['7/6', 12],
+    ['9/8', 12],
+    ['8/7', 13],
+    ['10/9', 13],
+    ['16/15', 16],
+    ['15/14', 17],
+    ['11/10', 18],
+    ['12/11', 18],
+    ['21/20', 19],
+    ['25/24', 19],
+    ['13/12', 20],
+    ['28/27', 20],
+    ['14/13', 22],
+    ['36/35', 22],
+    ['22/21', 23],
+    ['27/26', 24],
+    ['33/32', 24],
+    ['17/16', 25],
+    ['18/17', 25],
+    ['26/25', 25],
+    ['49/48', 25],
+    ['64/63', 25],
+    ['81/80', 25],
+    ['45/44', 26],
+    ['50/49', 26],
+    ['19/18', 27],
+    ['40/39', 27],
+    ['55/54', 27],
+    ['20/19', 28],
+    ['56/55', 29],
+    ['65/64', 30],
+    ['35/34', 31],
+    ['100/99', 31],
+    ['24/23', 32],
+    ['51/50', 32],
+    ['34/33', 33],
+    ['91/90', 33],
+    ['99/98', 33],
+    ['66/65', 34],
+    ['57/56', 35],
+    ['23/22', 36],
+    ['46/45', 36],
+    ['76/75', 36],
+    ['78/77', 36],
+    ['85/84', 36],
+    ['39/38', 37],
+    ['52/51', 37],
+    ['96/95', 37],
+    ['30/29', 39],
+    ['29/28', 40],
+    ['70/69', 40],
+    ['31/30', 41],
+    ['32/31', 41],
+    ['77/76', 41],
+    ['63/62', 46],
+    ['37/36', 47],
+    ['69/68', 47],
+    ['92/91', 47],
+    ['88/87', 49],
+    ['41/40', 52],
+    ['75/74', 52],
+    ['42/41', 53],
+    ['58/57', 53],
+    ['43/42', 55],
+    ['82/81', 55],
+    ['38/37', 58],
+    ['44/43', 58],
+    ['48/47', 58],
+    ['93/92', 61],
+    ['54/53', 64],
+    ['86/85', 67],
+    ['53/52', 70],
+    ['60/59', 71],
+    ['47/46', 72],
+    ['61/60', 73],
+    ['95/94', 73],
+    ['87/86', 77],
+    ['67/66', 83],
+    ['72/71', 83],
+    ['94/93', 83],
+    ['71/70', 85],
+    ['73/72', 85],
+    ['68/67', 88],
+    ['59/58', 90],
+    ['80/79', 92],
+    ['62/61', 94],
+    ['79/78', 97],
+    ['84/83', 97],
+    ['90/89', 102],
+    ['89/88', 106],
+    ['97/96', 110],
+    ['74/73', 112],
+    ['98/97', 113],
+    ['83/82', 126],
+  ])('Agrees with XenWiki on %s ~ %s', (fraction, height) => {
+    expect(wilsonHeight(fraction)).toBe(height);
+  });
+
+  it.skip('fuzzes the fraction property', () => {
+    for (let i = 0; i < 100; ++i) {
+      const n = Math.floor(Math.random() * 1073741823);
+      const d = Math.floor(Math.random() * 1073741822) + 1;
+      const f = new Fraction(n, d);
+      expect(wilsonHeight(f.n) + wilsonHeight(f.d)).toBe(wilsonHeight(f));
+    }
+  });
+
+  it.skip('fuzzes the multiplicative property', () => {
+    for (let i = 0; i < 10000; ++i) {
+      const x = Math.floor(Math.random() * 32768);
+      const y = Math.floor(Math.random() * 32768);
+      expect(wilsonHeight(x) + wilsonHeight(y), `Failed on ${x} * ${y}`).toBe(
+        wilsonHeight(x * y)
+      );
+    }
+  });
+});
diff --git a/src/__tests__/monzo.spec.ts b/src/__tests__/monzo.spec.ts
@@ -423,9 +423,25 @@ describe('Sparse monzos', () => {
   });
 
   it.each([
-    49305423, 4104956, 8375509, 27826943, 44852222, 22932439, 46933379,
-    59598447, 9693451, 54191546, 61834729, 26866018, 46410510, 47335837,
+    49305423,
+    4104956,
+    8375509,
+    27826943,
+    44852222,
+    22932439,
+    46933379,
+    59598447,
+    9693451,
+    54191546,
+    61834729,
+    26866018,
+    46410510,
+    47335837,
     43839566,
+    27470 * 5043,
+    17719 * 21909,
+    15812 * 27083,
+    21097 * 29468,
   ])('works on a tough case %s found during fuzzing', n => {
     const exponentByPrime = primeFactorize(n);
     let m = 1;

diff --git a/src/index.ts b/src/index.ts
@@ -1,6 +1,6 @@
-import {Fraction, mmod} from './fraction';
-import {Monzo, monzoToBigNumeratorDenominator} from './monzo';
-import {PRIME_CENTS} from './primes';
+import {Fraction, FractionValue, mmod} from './fraction';
+import {Monzo, monzoToBigNumeratorDenominator, primeFactorize} from './monzo';
+import {LOG_PRIMES, PRIMES, PRIME_CENTS} from './primes';
 import {sum} from './polyfills/sum-precise';
 
 export * from './fraction';
@@ -500,3 +500,50 @@ export function monzoToCents(monzo: Monzo) {
   }
   return Math.log1p(Number(delta) / Number(denominator)) * NATS_TO_CENTS;
 }
+
+/**
+ * Given fraction p/q calculate log(abs(p*q)).
+ * @param value Rational number or an array of its prime exponents.
+ * @returns The Tenney-height of the number.
+ */
+export function tenneyHeight(value: Monzo | FractionValue) {
+  if (Array.isArray(value)) {
+    return dotPrecise(
+      value.map(x => Math.abs(x)),
+      LOG_PRIMES
+    );
+  }
+  const {s, n, d} = new Fraction(value);
+  if (!s) {
+    return Infinity;
+  }
+  return Math.log(n) + Math.log(d);
+}
+
+/**
+ * Given fraction p/q calculate sopfr(p) + sopfr(q), ignoring sign.
+ * @param value Rational number, an array of its prime exponents or a `Map` of its prime exponents.
+ * @returns Sum of prime factors with repetition of p*q.
+ */
+export function wilsonHeight(
+  value: Monzo | FractionValue | Map<number, number>
+) {
+  if (Array.isArray(value)) {
+    return dot(
+      value.map(x => Math.abs(x)),
+      PRIMES
+    );
+  }
+  if (!(value instanceof Map)) {
+    value = primeFactorize(value);
+  }
+  if (value.has(0)) {
+    return Infinity;
+  }
+  value.delete(-1);
+  let result = 0;
+  for (const [prime, exponent] of value) {
+    result += Math.abs(exponent) * prime;
+  }
+  return result;
+}
diff --git a/src/monzo.ts b/src/monzo.ts
@@ -718,6 +718,18 @@ export function primeFactorize(value: FractionValue): Map<number, number> {
         } else if (subFactor === 1241) {
           result.set(17, (result.get(17) ?? 0) + 1);
           result.set(73, (result.get(73) ?? 0) + 1);
+        } else if (subFactor === 1681) {
+          result.set(41, (result.get(41) ?? 0) + 2);
+        } else if (subFactor === 3149) {
+          result.set(47, (result.get(47) ?? 0) + 1);
+          result.set(67, (result.get(67) ?? 0) + 1);
+        } else if (subFactor === 3869) {
+          result.set(53, (result.get(53) ?? 0) + 1);
+          result.set(73, (result.get(73) ?? 0) + 1);
+        } else if (subFactor === 65773) {
+          result.set(17, (result.get(17) ?? 0) + 1);
+          result.set(53, (result.get(53) ?? 0) + 1);
+          result.set(73, (result.get(73) ?? 0) + 1);
         } else {
           result.set(subFactor, (result.get(subFactor) ?? 0) + 1);
         }