Skip to content

Instantly share code, notes, and snippets.

@grondilu
Last active May 28, 2020 19:56
Show Gist options
  • Save grondilu/5916c9dbc3ec85da8372c8aa6315c951 to your computer and use it in GitHub Desktop.
Save grondilu/5916c9dbc3ec85da8372c8aa6315c951 to your computer and use it in GitHub Desktop.
Recursive matrix inversion
"use strict";
class Matrix {
constructor(a, b, c, d) {
if ([a,b,c,d].every(x => x instanceof Number || typeof(x) == "number")) {
this.n = 1;
} else if ([a,b,c,d].every(x => x instanceof Matrix && x.n == a.n)) {
this.n = a.n + 1;
} else throw new Error("wrong arguments type");
this.a = a;
this.b = b;
this.c = c;
this.d = d;
}
get inverse() {
let [a,b,c,d] = this.abcd;
if (this.n == 1) {
let det = a*d - b*c;
return det == 0 ? Infinity :
new Matrix(...[d,-b,-c,a].map(x => x/det));
} else {
/*
* a b
* c d
* x y = 1 0
* z w 0 1
*
* xa + yc = 1
* xb + yd = 0
* zb + wd = 1
* za + wc = 0
*
* x = 1/(a-b/d c)
* y = (1 - xa)/c
* z = 1/(b-a/c d)
* w = (1 - zb)/d
*
*/
let D = d.inverse;
if (D == Infinity) throw new Error("division by zero");
let x = a.add(b.mul(d.inverse).mul(c).neg).inverse,
y = identity(x.n).add(x.mul(a).neg).mul(c.inverse),
z = b.add(a.mul(c.inverse).mul(d).neg).inverse,
w = identity(z.n).add(z.mul(b).neg).mul(d.inverse);
return new Matrix(x, y, z, w);
}
}
get abcd() { return [this.a, this.b, this.c, this.d]; }
get neg() { return new Matrix(...this.abcd.map(x => this.n == 1 ? -x : x.neg)); }
add(that) {
if (that instanceof Number || typeof(that) == 'number')
return new Matrix(...this.abcd.map(x => x + that));
if (!(that instanceof Matrix))
throw new Error("wrong object");
if (this.n !== that.n)
throw new Error("size mistmatch");
return this.n == 1 ?
new Matrix(
this.a + that.a,
this.b + that.b,
this.c + that.c,
this.d + that.d
) :
new Matrix(
this.a.add(that.a),
this.b.add(that.b),
this.c.add(that.c),
this.d.add(that.d)
);
}
mul(that) {
if (!(that instanceof Matrix))
throw new Error("wrong object");
if (this.n !== that.n)
throw new Error("size mistmatch");
let [a,b,c,d] = this.abcd,
[x,y,z,w] = that.abcd;
return this.n == 1 ?
new Matrix(
a*x+b*z,a*y+b*w,
c*x+d*z,c*y+d*w
) :
new Matrix(
a.mul(x).add(b.mul(z)),a.mul(y).add(b.mul(w)),
c.mul(x).add(d.mul(z)),c.mul(y).add(d.mul(w))
);
}
get norm() {
return this.n == 1 ? this.abcd.map(x => x*x).reduce((a,b) => a+b) :
this.abcd.map(x => x.norm).reduce((a,b) => a+b);
}
}
function rand(n, r = Math.random) {
return n == 1 ?
new Matrix(...[1,2,3,4].map(r)) :
new Matrix(...[1,2,3,4].map(x => rand(n-1)));
}
function zeros(n) {
if (n == 1) return new Matrix(0, 0, 0, 0);
let z = zeros(n-1);
return new Matrix(z, z, z, z);
}
function identity(n) {
if (n == 1) return new Matrix(1, 0, 0, 1);
let i = identity(n-1),
z = zeros(n-1);
return new Matrix(i, z, z, i);
}
const n = 5;
let A = rand(n);
console.log(
A.mul(A.inverse).add(identity(n).neg).norm
);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment