Implement modular inverse

Add an option for non-strict behavior in the context of multiple interleaved cycles.

src/__tests__/index.spec.ts

@@ -22,6 +22,8 @@ import {
   monzoToCents,
   tenneyHeight,
   wilsonHeight,
+  modInv,
+  mmod,
 } from '../index';
 
 const FUZZ = 'FUZZ' in process.env;
@@ -621,3 +623,59 @@ describe('Wilson complexity measure', () => {
     }
   });
 });
+
+describe('Modular inverse calculator', () => {
+  it('finds modular inverses when they exist', () => {
+    for (let a = -30; a < 30; ++a) {
+      for (let b = 2; b < 20; ++b) {
+        const c = Math.abs(gcd(a, b));
+        if (c !== 1) {
+          expect(() => modInv(a, b)).toThrow();
+        } else {
+          expect(mmod(a * modInv(a, b), b)).toBe(1);
+        }
+      }
+    }
+  });
+
+  it("finds the next best alternative even if a true modular inverse doesn't exist", () => {
+    for (let a = -30; a < 30; ++a) {
+      for (let b = 1; b < 20; ++b) {
+        let c = Math.abs(gcd(a, b));
+        if (c === b) {
+          c = 0;
+        }
+        expect(mmod(a * modInv(a, b, false), b)).toBe(c);
+      }
+    }
+  });
+
+  it('can be used to break 15edo into 3 interleaved circles of fifths', () => {
+    const edo = 15;
+    const gen = 9;
+    const genInv = modInv(gen, edo, false);
+    const numCycles = mmod(gen * genInv, edo);
+    function f(step: number) {
+      const c = mmod(step, numCycles);
+      return mmod((step - c) * genInv, edo) + c;
+    }
+    expect([...Array(16).keys()].map(f)).toEqual([
+      0, // Root
+      1, // Root + 1
+      2, // Root + 2
+      6, // Gen * 2
+      7, // +1
+      8, // +2
+      12, // Gen * 4
+      13, // +1
+      14, // +2
+      3, // Gen
+      4, // Gen + 1
+      5, // Gen + 1
+      9, // Gen * 3
+      10, // +1
+      11, // +2
+      0, // Circle closes
+    ]);
+  });
+});

src/index.ts

@@ -128,6 +128,23 @@ export function iteratedEuclid(params: Iterable<number>) {
   return coefs;
 }
 
+/**
+ * Find modular inverse of a (mod b).
+ * @param a Number to find modular inverse of.
+ * @param b Modulus.
+ * @param strict Ensure that a * modInv(a, b) = 1 (mod b). If `strict = false` we have a * modInv(a, b) = gdc(a, b) (mod b) instead.
+ * @returns The modular inverse in the range {0, 1, ..., b - 1}.
+ */
+export function modInv(a: number, b: number, strict = true) {
+  const {gcd, coefA} = extendedEuclid(a, b);
+  if (strict && gcd !== 1) {
+    throw new Error(
+      `${a} does not have a modular inverse modulo ${b} since they're not coprime`
+    );
+  }
+  return mmod(coefA, b);
+}
+
 /**
  * Collection of unique fractions.
  */