From 6b69940896a6582b98b7c2c37e9fc46b324da468 Mon Sep 17 00:00:00 2001
From: Lumi Pakkanen <lumi.pakkanen@gmail.com>
Date: Thu, 2 May 2024 07:53:14 +0300
Subject: [PATCH] Remove residual restrictions from monzo dot product

---
 src/__tests__/monzo.spec.ts | 25 +++++++++++++++++++++++
 src/monzo.ts                | 40 ++++++++++++++++++++++++++++++++++++-
 2 files changed, 64 insertions(+), 1 deletion(-)

diff --git a/src/__tests__/monzo.spec.ts b/src/__tests__/monzo.spec.ts
index fa1056b7..38b9bc4e 100644
--- a/src/__tests__/monzo.spec.ts
+++ b/src/__tests__/monzo.spec.ts
@@ -425,4 +425,29 @@ describe('Extended monzo with time', () => {
       90, 120, 180, 360,
     ]);
   });
+
+  it('can calculate the dot product between syntonic and porcupine commas (prime exponents)', () => {
+    const neg27 = TimeMonzo.fromFraction('81/80').dot(
+      TimeMonzo.fromFraction('250/243')
+    );
+    expect(neg27.toFraction()).toBe('-27');
+  });
+
+  it('can calculate the dot product between syntonic and porcupine commas (residuals)', () => {
+    const neg27 = new TimeMonzo(new Fraction(0), [], new Fraction('81/80')).dot(
+      new TimeMonzo(new Fraction(0), [], new Fraction('250/243'))
+    );
+    expect(neg27.toFraction()).toBe('-27');
+  });
+
+  it('can calculate the dot product between syntonic and porcupine commas (mixed)', () => {
+    const neg27 = new TimeMonzo(
+      new Fraction(0),
+      [new Fraction(-4)],
+      new Fraction('81/5')
+    ).dot(
+      new TimeMonzo(new Fraction(0), [new Fraction(1)], new Fraction('125/243'))
+    );
+    expect(neg27.toFraction()).toBe('-27');
+  });
 });
diff --git a/src/monzo.ts b/src/monzo.ts
index 4805fa1c..e98e319e 100644
--- a/src/monzo.ts
+++ b/src/monzo.ts
@@ -13,6 +13,7 @@ import {
   BIG_INT_PRIMES,
   monzoToCents,
   sum,
+  primeFactorize,
 } from 'xen-dev-utils';
 
 import {
@@ -137,6 +138,30 @@ function min(a: Fraction, b: Fraction) {
   return a;
 }
 
+function intDot(a: number, b: number) {
+  const factor = gcd(a, b);
+  if (factor === 1) {
+    return 0;
+  }
+  a /= factor;
+  b /= factor;
+  let total = 0;
+  for (const [prime, exponent] of primeFactorize(factor)) {
+    let aCount = exponent;
+    while (!(a % prime)) {
+      ++aCount;
+      a /= prime;
+    }
+    let bCount = exponent;
+    while (!(b % prime)) {
+      ++bCount;
+      b /= prime;
+    }
+    total += aCount * bCount;
+  }
+  return total;
+}
+
 /**
  * Arbitrary (but inaccurate) value measured in time-related units (usually Hz).
  *
@@ -1887,7 +1912,20 @@ export class TimeMonzo {
       }
       return result;
     }
-    throw new Error('Residuals prevent calculating the dot product.');
+    if (!this.residual.s && !other.residual.s) {
+      throw new Error('Dot product of 0 is ambiguous.');
+    }
+    if (this.residual.s < 0 && other.residual.s < 0) {
+      throw new Error('Dot product of -1 is ambiguous.');
+    }
+    const {n, d} = this.residual;
+    const {n: p, d: q} = other.residual;
+    const resDot = intDot(n, p) - intDot(n, q) - intDot(d, p) + intDot(d, q);
+    let result = this.timeExponent.mul(other.timeExponent).add(resDot);
+    for (let i = 0; i < this.primeExponents.length; ++i) {
+      result = result.add(this.primeExponents[i].mul(other.primeExponents[i]));
+    }
+    return result;
   }
 
   /**