Skip to content

Commit

Permalink
Remove residual restrictions from monzo dot product
Browse files Browse the repository at this point in the history
  • Loading branch information
frostburn committed May 2, 2024
1 parent 3565efc commit 6b69940
Show file tree
Hide file tree
Showing 2 changed files with 64 additions and 1 deletion.
25 changes: 25 additions & 0 deletions src/__tests__/monzo.spec.ts
Original file line number Diff line number Diff line change
Expand Up @@ -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');
});
});
40 changes: 39 additions & 1 deletion src/monzo.ts
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ import {
BIG_INT_PRIMES,
monzoToCents,
sum,
primeFactorize,
} from 'xen-dev-utils';

import {
Expand Down Expand Up @@ -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).
*
Expand Down Expand Up @@ -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;
}

/**
Expand Down

0 comments on commit 6b69940

Please sign in to comment.