Implement Hermite normal form and associated methods

src/__tests__/basis.spec.ts

@@ -20,7 +20,9 @@ import {
   matsub,
   matmul,
   minor,
-  transpose,
+  canonical,
+  respell,
+  solveDiophantine,
 } from '../basis';
 import {
   FractionalMonzo,
@@ -35,6 +37,7 @@ import {
 import {dot} from '../number-array';
 import {LOG_PRIMES, PRIMES} from '../primes';
 import {Fraction} from '../fraction';
+import {cokernel, kernel, preimage, transpose} from '../hnf';
 
 const FUZZ = 'FUZZ' in process.env;
 
@@ -747,14 +750,6 @@ describe('Determinant', () => {
 });
 
 describe('Transpose', () => {
-  it('transposes a 3x2 matrix', () => {
-    const mat = [[1, 2], [3], [4, 5]];
-    expect(transpose(mat)).toEqual([
-      [1, 3, 4],
-      [2, 0, 5],
-    ]);
-  });
-
   it('transposes a 3x2 matrix with rational entries', () => {
     const mat = [[1, 0.5], [3], ['2/7', 5]];
     expect(
@@ -826,3 +821,87 @@ describe('Auxiliary matrix methods', () => {
     ]);
   });
 });
+
+describe("Sin-tel's example", () => {
+  it('agrees with temper/example.py', () => {
+    // 31 & 19
+    const tMap = [
+      [19, 30, 44, 53],
+      [31, 49, 72, 87],
+    ];
+    // Canonical form
+    const c = [
+      [1, 0, -4, -13],
+      [0, 1, 4, 10],
+    ];
+    const can = canonical(tMap);
+    expect(can).toEqual(c);
+
+    let commas = transpose(kernel(can));
+    expect(commas.map(comma => monzoToFraction(comma).toFraction())).toEqual([
+      '126/125',
+      '1275989841/1220703125',
+    ]);
+
+    commas = lenstraLenstraLovasz(commas).basis;
+    expect(commas.map(comma => monzoToFraction(comma).toFraction())).toEqual([
+      '126/125',
+      '81/80',
+    ]);
+
+    const gens = transpose(preimage(can));
+
+    expect(gens.map(gen => monzoToFraction(gen).toFraction())).toEqual([
+      '20253807/9765625',
+      '3',
+    ]);
+
+    const simpleGens = gens.map(gen => respell(gen, commas));
+
+    expect(simpleGens.map(gen => monzoToFraction(gen).toFraction())).toEqual([
+      '2',
+      '3',
+    ]);
+  });
+
+  it('agrees with temper/example_commas.py', () => {
+    const commas = ['81/80', '225/224'].map(toMonzo);
+
+    expect(commas).toEqual([
+      [-4, 4, -1],
+      [-5, 2, 2, -1],
+    ]);
+
+    const tMap = cokernel(transpose(commas));
+    expect(tMap).toEqual([
+      [1, 0, -4, -13],
+      [0, 1, 4, 10],
+    ]);
+  });
+});
+
+describe('Diophantine equation solver', () => {
+  it.skip('solves a diophantine equation', () => {
+    const A = [
+      [1, 2, 3],
+      [-4, 5, 6],
+      [7, -8, 9],
+    ];
+    const b = [4, 28, -28];
+    const solution = solveDiophantine(A, b);
+    expect(solution).toEqual([-3, 2, 1]);
+  });
+
+  it('converts standard basis commas to subgroup basis monzos (pinkan)', () => {
+    const subgroup = '2.3.13/5.19/5'.split('.').map(toMonzo);
+    const commas = ['676/675', '1216/1215'].map(toMonzo);
+    const subgroupMonzos = solveDiophantine(
+      transpose(subgroup),
+      transpose(commas)
+    );
+    expect(transpose(subgroupMonzos)).toEqual([
+      [2, -3, 2, 0],
+      [6, -5, 0, 1],
+    ]);
+  });
+});

src/__tests__/hnf.spec.ts

@@ -0,0 +1,97 @@
+import {describe, it, expect} from 'vitest';
+import {antitranspose, hnf, integerDet, preimage, transpose} from '../hnf';
+
+describe('Hermite normal form', () => {
+  it('works on a small matrix', () => {
+    const A = [
+      [2, 3, 6, 2],
+      [5, 6, 1, 6],
+      [8, 3, 1, 1],
+    ];
+
+    const H = hnf(A);
+
+    expect(H).toEqual([
+      [1, 0, 50, -11],
+      [0, 3, 28, -2],
+      [0, 0, 61, -13],
+    ]);
+  });
+
+  it('works on a small matrix (bigint)', () => {
+    const A = [
+      [2n, 3n, 6n, 2n],
+      [5n, 6n, 1n, 6n],
+      [8n, 3n, 1n, 1n],
+    ];
+
+    const H = hnf(A);
+
+    expect(H).toEqual([
+      [1n, 0n, 50n, -11n],
+      [0n, 3n, 28n, -2n],
+      [0n, 0n, 61n, -13n],
+    ]);
+  });
+});
+
+describe('Integer determinant', () => {
+  it('works on a small matrix', () => {
+    const mat = [
+      [-2, -1, 2],
+      [2, 1, 4],
+      [-3, 3, -1],
+    ];
+    const determinant = integerDet(mat);
+    expect(determinant).toBe(54);
+  });
+
+  it('works on a small matrix (bigint)', () => {
+    const mat = [
+      [-2n, -1n, 2n],
+      [2n, 1n, 4n],
+      [-3n, 3n, -1n],
+    ];
+    const determinant = integerDet(mat);
+    expect(determinant).toBe(54n);
+  });
+});
+
+describe('Transpose', () => {
+  it('transposes a 3x2 matrix', () => {
+    const mat = [[1, 2], [3], [4, 5]];
+    expect(transpose(mat)).toEqual([
+      [1, 3, 4],
+      [2, 0, 5],
+    ]);
+  });
+});
+
+describe('Anti-transpose', () => {
+  it('anti-transposes a 3x2 matrix', () => {
+    const mat = [
+      [1, 2],
+      [3, 4],
+      [5, 6],
+    ];
+    expect(antitranspose(mat)).toEqual([
+      [6, 4, 2],
+      [5, 3, 1],
+    ]);
+  });
+});
+
+describe('Preimage', () => {
+  it('obtains the preimage associated with 5-limit meantone', () => {
+    const map = [
+      [1, 0, -4],
+      [0, 1, 4],
+    ];
+    const gensT = preimage(map);
+    expect(gensT).toEqual([
+      [1, 0],
+      [0, 1],
+      [0, 0],
+    ]);
+  });
+});