From a7bc711d066d5215672684db8d1ef79588da285b Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Sat, 29 Jun 2024 17:44:54 +0300 Subject: [PATCH] Implement Hermite normal form and associated methods ref #43 --- src/__tests__/basis.spec.ts | 97 ++++++++++- src/__tests__/hnf.spec.ts | 97 +++++++++++ src/basis.ts | 159 +++++++++++++++--- src/hnf.ts | 321 ++++++++++++++++++++++++++++++++++++ src/index.ts | 1 + 5 files changed, 645 insertions(+), 30 deletions(-) create mode 100644 src/__tests__/hnf.spec.ts create mode 100644 src/hnf.ts diff --git a/src/__tests__/basis.spec.ts b/src/__tests__/basis.spec.ts index aeec2d7..ef33ee7 100644 --- a/src/__tests__/basis.spec.ts +++ b/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], + ]); + }); +}); diff --git a/src/__tests__/hnf.spec.ts b/src/__tests__/hnf.spec.ts new file mode 100644 index 0000000..e0b00cd --- /dev/null +++ b/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], + ]); + }); +}); diff --git a/src/basis.ts b/src/basis.ts index fd045a5..d4bff75 100644 --- a/src/basis.ts +++ b/src/basis.ts @@ -1,16 +1,19 @@ import {Fraction, FractionValue} from './fraction'; import { FractionalMonzo, + Monzo, ProtoFractionalMonzo, add, fractionalAdd, fractionalDot, fractionalScale, fractionalSub, + monzosEqual, scale, sub, } from './monzo'; import {dot} from './number-array'; +import {transpose, hnf, integerDet, antitranspose, padMatrix} from './hnf'; /** * Result of Gram–Schmidt process without normalization. @@ -694,27 +697,6 @@ export function fractionalDet(matrix: ProtoFractionalMonzo[]) { return result; } -/** - * Transpose a 2-D matrix (swap rows and columns). - * @param matrix Matrix to transpose. - * @returns The transpose. - */ -export function transpose(matrix: number[][]) { - let width = 0; - for (const row of matrix) { - width = Math.max(row.length, width); - } - const result: number[][] = []; - for (let i = 0; i < width; ++i) { - const row: number[] = []; - for (let j = 0; j < matrix.length; ++j) { - row.push(matrix[j][i] ?? 0); - } - result.push(row); - } - return result; -} - /** * Transpose a 2-D matrix with rational entries (swap rows and columns). * @param matrix Matrix to transpose. @@ -840,3 +822,138 @@ export function fractionalMatsub( } return result; } + +// Finds the Hermite normal form and 'defactors' it. +// Defactoring is also known as saturation. +// This removes torsion from the map. +// Algorithm as described by: +// +// Clément Pernet and William Stein. +// Fast Computation of HNF of Random Integer Matrices. +// Journal of Number Theory. +// https://doi.org/10.1016/j.jnt.2010.01.017 +// See section 8. + +/** + * Compute the Hermite normal form with torsion removed. + * @param M The input matrix. + * @returns The defactored Hermite normal form. + */ +export function defactoredHnf(M: number[][]): number[][] { + // Need to convert to bigint so that intermediate results don't blow up. + const bigM = M.map(row => row.map(BigInt)); + const K = hnf(transpose(bigM)); + while (K.length > M.length) { + K.pop(); + } + const determinant = integerDet(K); + if (determinant === 1n) { + return hnf(bigM).map(row => row.map(Number)); + } + const S = inv(transpose(K).map(row => row.map(Number))); + const D = matmul(S, M).map(row => row.map(Math.round)); + return hnf(D); +} + +/** + * Compute the canonical form of the input. + * @param M Input maps. + * @returns Defactored Hermite normal form or the antitranspose sandwich for commas bases. + */ +export function canonical(M: number[][]): number[][] { + for (const row of M) { + if (row.length < M.length) { + // Comma basis + return antitranspose(defactoredHnf(antitranspose(M))); + } + } + return defactoredHnf(M); +} + +// Babai's nearest plane algorithm for solving approximate CVP +// `basis` should be LLL reduced first + +/** + * Solve approximate CVP using Babai's nearest plane algorithm. + * @param v Vector to reduce. + * @param basis LLL basis to reduce with. + * @param dual Optional precalculated geometric duals of the orthogonalized basis. + * @returns The reduced vector. + */ +export function nearestPlane( + v: number[], + basis: number[][], + dual?: number[][] +) { + // Body moved to respell to save a sub() call. + return sub(v, respell(v, basis, dual)); +} + +/** + * Respell a monzo represting a rational number to a simpler form. + * @param monzo Array of prime exponents to simplify. + * @param commas Monzos representing near-unisons to simplify by. Should be LLL reduced to work properly. + * @param commaOrthoDuals Optional precalculated geometric duals of the orthogonalized comma basis. + * @returns An array of prime exponents representing a simpler rational number. + */ +export function respell( + monzo: Monzo, + commas: Monzo[], + commaOrthoDuals?: Monzo[] +) { + if (commaOrthoDuals === undefined) { + commaOrthoDuals = gram(commas).dual; + } + monzo = [...monzo]; + + for (let i = commaOrthoDuals.length - 1; i >= 0; --i) { + const mu = dot(monzo, commaOrthoDuals[i]); + monzo = sub(monzo, scale(commas[i], Math.round(mu))); + } + return monzo; +} + +/** + * Solve Ax = b in the integers. + * @param A Coefficients of unknowns. + * @param b Target vector or a matrix of column vectors. + * @returns A vector mapped to the target vector by A or a matrix of row vectors. + */ +export function solveDiophantine(A: number[][], b: number[][]): number[][]; +export function solveDiophantine(A: number[][], b: number[]): number[]; +export function solveDiophantine( + A: number[][], + b: number[] | number[][] +): typeof b { + const hasMultiple = Array.isArray(b[0]); + + const B: number[][] = hasMultiple + ? padMatrix(b as number[][]).M + : (b as number[]).map(c => [c]); + + // Need to convert to bigint so that intermediate results don't blow up. + const {width, height, M} = padMatrix(A.map(row => row.map(BigInt))); + for (let i = 0; i < height; ++i) { + M[i].push(...B[i].map(BigInt)); + } + const H = hnf(M); + while (H.length > width) { + H.pop(); + } + const c: number[][] = []; + for (const row of H) { + c.push(row.splice(width, row.length - width).map(Number)); + } + const S = inv(H.map(row => row.map(Number))); + const sol = matmul(S, c).map(row => row.map(Math.round)); + + // Check solution(s). + const BS = matmul(A, sol); + for (let i = 0; i < B.length; ++i) { + if (!monzosEqual(B[i], BS[i])) { + throw new Error('Could not solve system'); + } + } + + return hasMultiple ? sol : sol.map(row => row[0]); +} diff --git a/src/hnf.ts b/src/hnf.ts new file mode 100644 index 0000000..4e4a686 --- /dev/null +++ b/src/hnf.ts @@ -0,0 +1,321 @@ +/** + * Algorithm adapted from https://github.com/lan496/hsnf + * Guaranteed to not overflow with BigInt matrices + */ +import {bigAbs} from './monzo'; + +function abs(x: T): T { + if (typeof x === 'number') { + return Math.abs(x) as T; + } + return bigAbs(x) as T; +} + +function floorDiv(x: T, y: T): T { + if (typeof x === 'number') { + return Math.floor(x / (y as typeof x)) as T; + } + if (x >= 0n !== y >= 0n && x % (y as typeof x)) { + // @ts-ignore + return x / y - 1n; + } + return (x / y) as T; +} + +function getPivot(A: T[][], i1: number, j: number) { + let idx: number | undefined; + let valmin: number | bigint | undefined; + + for (let i = i1; i < A.length; ++i) { + if (!A[i][j]) { + continue; + } + if (valmin === undefined || abs(A[i][j]) < valmin) { + idx = i; + valmin = abs(A[i][j]); + } + } + return idx; +} + +/** + * Fix a 2-D matrix to have full rows (pad with zeros). + * @param M Input matrix. + * @returns Height, width, the padded matrix and the corresponding 0 or 0n and 1 or 1n. + */ +export function padMatrix(M: T[][]) { + const height = M.length; + if (!height) { + return { + height, + width: 0, + zero: 0 as T, + one: 1 as T, + M: [], + }; + } + let width = 0; + let zero: T | undefined; + for (const row of M) { + width = Math.max(width, row.length); + if (row.length) { + // @ts-ignore + zero = typeof row[0] === 'number' ? 0 : 0n; + } + } + M = M.map(row => [...row]); + if (zero === undefined) { + return { + height, + width, + zero: 0 as T, + one: 1 as T, + M, + }; + } + const one = typeof zero === 'number' ? 1 : 1n; + for (const row of M) { + while (row.length < width) { + row.push(zero); + } + } + return { + height, + width, + zero, + one, + M, + }; +} + +/** + * Compute the Hermite normal form of a matrix of integers. + * @param A The input matrix. + * @returns The input in Hermite normal form. + */ +export function hnf(A: T[][]): T[][] { + const {width, height, zero, one, M} = padMatrix(A); + let si = 0; + let sj = 0; + + // @ts-ignore + const negOne: T = -one; + + while (true) { + if (si === height || sj === width) { + return M; + } + + // choose a pivot + const row = getPivot(M, si, sj); + + if (row === undefined) { + // if there does not remain non-zero elements, go to a next column + sj = sj + 1; + continue; + } + // swap + [M[si], M[row]] = [M[row], M[si]]; + + // eliminate the s-th column entries + for (let i = si + 1; i < height; ++i) { + if (M[i][sj]) { + const k = floorDiv(M[i][sj], M[si][sj]); + for (let j = 0; j < width; ++j) { + // @ts-ignore + M[i][j] -= k * M[si][j]; + } + } + } + + // if there does not remain non-zero element in s-th column, find a next entry + let rowDone = true; + for (let i = si + 1; i < height; ++i) { + if (M[i][sj]) { + rowDone = false; + } + } + if (rowDone) { + if (M[si][sj] < zero) { + for (let j = 0; j < width; ++j) { + // @ts-ignore + M[si][j] *= negOne; + } + } + + if (M[si][sj]) { + for (let i = 0; i < si; ++i) { + const k = floorDiv(M[i][sj], M[si][sj]); + if (k) { + for (let j = 0; j < width; ++j) { + // @ts-ignore + M[i][j] -= k * M[si][j]; + } + } + } + } + + si += 1; + sj += 1; + } + } +} + +/** + * Prune rows filled with falsy values from a 2-D matrix. + * @param A: Matrix to prune in-place. + */ +export function pruneZeroRows(A: any[][]) { + for (let i = 0; i < A.length; ++i) { + if (!A[i].some(Boolean)) { + A.splice(i, 1); + i--; + } + } +} + +// exact integer determinant using Bareiss algorithm +// modified slightly from: +// https://stackoverflow.com/questions/66192894/precise-determinant-of-integer-nxn-matrix + +/** + * Compute the determinant of a matrix of integers. + * @param A The input matrix. + * @returns The determinant. + */ +export function integerDet(A: T[][]): T { + const {width, height, zero, one, M} = padMatrix(A); + if (!width || !height || width !== height) { + return zero; + } + let sign = one; + let prev = sign; + for (let i = 0; i < width - 1; ++i) { + if (!M[i][i]) { + // swap with another row having nonzero i's elem + let swapto: number | undefined; + for (let j = i + 1; j < height; ++j) { + if (M[j][i]) { + swapto = j; + break; + } + } + if (swapto === undefined) { + return zero; // all M[*][i] are zero => zero determinant + } + // swap rows + [M[i], M[swapto]] = [M[swapto], M[i]]; + sign = -sign; + } + for (let j = i + 1; j < height; ++j) { + for (let k = i + 1; k < width; ++k) { + // assert (M[j, k] * M[i, i] - M[j, i] * M[i, k]) % prev == 0 + // @ts-ignore + M[j][k] = floorDiv(M[j][k] * M[i][i] - M[j][i] * M[i][k], prev); + } + } + prev = M[i][i]; + } + // @ts-ignore + return sign * M.pop()!.pop()!; +} + +/** + * Anti-transpose a 2-D matrix (swap rows and columns along the other diagonal). + * @param matrix Matrix to antitranspose. + * @returns The antitranspose. + */ +export function antitranspose(matrix: T[][]): T[][] { + const {width, height, M} = padMatrix(matrix); + const result: T[][] = []; + for (let i = width - 1; i >= 0; --i) { + const row: T[] = []; + for (let j = height - 1; j >= 0; --j) { + row.push(M[j][i]); + } + result.push(row); + } + return result; +} + +/** + * Transpose a 2-D matrix (swap rows and columns). + * @param matrix Matrix to transpose. + * @returns The transpose. + */ +export function transpose(matrix: T[][]): T[][] { + const {width, height, M} = padMatrix(matrix); + const result: T[][] = []; + for (let i = 0; i < width; ++i) { + const row: T[] = []; + for (let j = 0; j < height; ++j) { + row.push(M[j][i]); + } + result.push(row); + } + return result; +} + +// Find the left kernel (nullspace) of M. +// Adjoin an identity block matrix and solve for HNF. +// This is equivalent to the highschool maths method, +// but using HNF instead of Gaussian elimination. + +/** + * Find the left kernel (nullspace) of the input matrix. + * @param A The input matrix. + * @returns The kernel matrix. + */ +export function kernel(A: T[][]): T[][] { + const {width, height, zero, one, M} = padMatrix(A); + for (let i = 0; i < width; ++i) { + const row = Array(width).fill(zero); + row[i] = one; + M.push(row); + } + const K = transpose(hnf(transpose(M))); + for (let i = 0; i < height; ++i) { + K.shift(); + } + for (const row of K) { + for (let i = 0; i < height; ++i) { + row.shift(); + } + } + return K; +} + +/** + * Find the right kernel (nullspace) of the input matrix. + * @param A The input matrix. + * @returns The kernel matrix. + */ +export function cokernel(A: T[][]): T[][] { + return transpose(kernel(transpose(A))); +} + +/** + * Find the preimage X of A such that AX = I. + * @param A The input matrix. + * @returns The preimage. + */ +export function preimage(A: T[][]): T[][] { + const {width, height, zero, one, M} = padMatrix(A); + const B = transpose(M); + for (let i = 0; i < width; ++i) { + for (let j = 0; j < width; ++j) { + // @ts-ignore + B[i].push(i === j ? one : zero); + } + } + const H = hnf(B); + while (H.length > height) { + H.pop(); + } + for (const row of H) { + for (let i = 0; i < height; ++i) { + row.shift(); + } + } + return transpose(H); +} diff --git a/src/index.ts b/src/index.ts index 29ddb40..eaa5df3 100644 --- a/src/index.ts +++ b/src/index.ts @@ -11,6 +11,7 @@ export * from './monzo'; export * from './approximation'; export * from './number-array'; export * from './basis'; +export * from './hnf'; export {sum} from './polyfills/sum-precise'; export interface AnyArray {