feat: Implement Vector type (#28)

The Vector type is simply a generic type for built-in JavaScript Arrays. The parameter type is the type of all of the entries of the Array. The Vector type also supports inhomogeneous arrays by using the special type `Unknown` as the argument type, but note that when computing with inhomogeneous arrays, method dispatch must be performed separately for every entry in a calculation, making all operations considerably slower than on homogeneous Vector instances.

Note also that arithmetic operations on nested Vectors, e.g. `Vector(Vector(NumberT))` are defined so as to interpret such entities as ordinary matrices, represented in row-major format (i.e., the component `Vector(NumberT)` items of such an entity are the _rows_ of the matrix.

Reviewed-on: #28
Co-authored-by: Glen Whitney <glen@studioinfinity.org>
Co-committed-by: Glen Whitney <glen@studioinfinity.org>

src/vector/Vector.js

@@ -0,0 +1,48 @@
+import {Type, Unknown} from '#core/Type.js'
+
+const isVector = v => Array.isArray(v)
+
+export const Vector = new Type(isVector, {
+   specialize(CompType) {
+      const compTest = CompType.test
+      const specTest = v => isVector(v) && v.every(compTest)
+      const typeName = `Vector(${CompType})`
+      if (CompType.concrete) {
+         const typeOptions = {typeName}
+         if ('zero' in CompType) typeOptions.zero = [CompType.zero]
+         const vectorCompType = new Type(specTest, typeOptions)
+         vectorCompType.Component = CompType
+         vectorCompType.vectorDepth = (CompType.vectorDepth ?? 0) + 1
+         return vectorCompType
+      }
+      // Wrapping a generic type in Vector
+      return new Type(specTest, {
+         typeName,
+         specialize: (...args) => this.specialize(CompType.specialize(...args)),
+         specializesTo: VT => this.specializesTo(VT)
+            && CompType.specializesTo(VT.Component),
+         refine: (v, typer) => {
+            if (!v.length) {
+               throw new RangeError(
+                  `no way to refine ${typeName} on an empty vector`)
+            }
+            const eltTypes = v.map(elt => CompType.refine(elt, typer))
+            const newCompType = eltTypes[0]
+            if (eltTypes.some(T => T !== newCompType)) {
+               throw new TypeError(
+                  `can't refine ${typeName} to ${v}; `
+                  + `not all entries have type ${newCompType}`)
+            }
+            return this.specialize(newCompType)
+         }
+      })
+   },
+   specializesTo: VT => 'vectorDepth' in VT && VT.vectorDepth > 0,
+   refine(v, typer) {
+      if (!v.length) return this.specialize(Unknown) // what else could we do?
+      const eltTypes = v.map(elt => typer(elt))
+      let compType = eltTypes[0]
+      if (eltTypes.some(T => T !== compType)) compType = Unknown
+      return this.specialize(compType)
+   }
+})

src/vector/__test__/Vector.spec.js

@@ -0,0 +1,35 @@
+import assert from 'assert'
+import {Vector} from '../Vector.js'
+import {Unknown} from '#core/Type.js'
+import math from '#nanomath'
+
+describe('Vector Type', () => {
+   it('correctly recognizes vectors', () => {
+      assert(Vector.test([3, 4]))
+      assert(!Vector.test({re: 3, im: 4}))
+   })
+   it('can fully type vectors', () => {
+      const {BooleanT, Complex, NumberT} = math.types
+      const typ = math.typeOf
+      const cplx = math.complex
+      assert.strictEqual(typ([3, 4]), Vector(NumberT))
+      assert.strictEqual(typ([true, false]), Vector(BooleanT))
+      assert.strictEqual(typ([3, false]), Vector(Unknown))
+      assert.strictEqual(typ([cplx(3, 4), cplx(5)]), Vector(Complex(NumberT)))
+      assert.strictEqual(typ([[3, 4], [5]]), Vector(Vector(NumberT)))
+   })
+   it('can refine when nested', () => {
+      const {Complex, NumberT} = math.types
+      const cplx = math.complex
+      const VecComplex = Vector(Complex)
+      const vcEx = [cplx(3, 4), cplx(5)]
+      assert(VecComplex.test(vcEx))
+      const vcType = VecComplex.refine(vcEx, math.typeOf)
+      assert.strictEqual(vcType, Vector(Complex(NumberT)))
+      const CplxVec = Complex(Vector)
+      const cvEx = cplx([3,4], [5, 0])
+      assert(CplxVec.test(cvEx))
+      const cvType = CplxVec.refine(cvEx, math.typeOf)
+      assert.strictEqual(cvType, Complex(Vector(NumberT)))
+   })
+})

src/vector/__test__/arithmetic.spec.js

@@ -0,0 +1,82 @@
+import assert from 'assert'
+import math from '#nanomath'
+
+describe('Vector arithmetic functions', () => {
+   const cplx = math.complex
+   const cV = [cplx(3, 4), cplx(4, -3), cplx(-3, -4)]
+   it('distributes abs elementwise', () => {
+      const abs = math.abs
+      assert.deepStrictEqual(abs([-3, 4, -5]), [3, 4, 5])
+      assert.deepStrictEqual(abs(cV), [5, 5, 5])
+      assert.deepStrictEqual(abs([true, -4, cplx(4, 3)]), [1, 4, 5])
+   })
+   it('computes the norm of a vector or matrix', () => {
+      const norm = math.norm
+      assert.strictEqual(norm([-3, 4, -5]), Math.sqrt(50))
+      assert.strictEqual(norm([[1, 2], [4, 2]]), 5)
+      assert.strictEqual(norm(cV), Math.sqrt(75))
+   })
+   it('adds vectors and matrices', () => {
+      const add = math.add
+      assert.deepStrictEqual(add([-3, 4, -5], [8, 0, 2]), [5, 4, -3])
+      assert.deepStrictEqual(
+         add([[1, 2], [4, 2]], [[0, -1], [3, -4]]), [[1, 1], [7, -2]])
+      assert.deepStrictEqual(
+         add([[1, 2], [4, 2]], [0, -1]), [[1, 1], [4, 1]])
+   })
+   it('(pseudo)inverts matrices', () => {
+      const inv = math.invert
+      // inverses
+      assert.deepStrictEqual(inv([3, 4, 5]), [3/50, 2/25, 1/10])
+      assert.deepStrictEqual(inv([[4]]), [[0.25]])
+      assert.deepStrictEqual(inv([[5, 2], [-7, -3]]), [[3, 2], [-7, -5]])
+      assert(math.equal(
+         inv([[3, 0, 2], [2, 1, 0], [1, 4, 2]]),
+         [[1/10, 2/5, -1/10], [-1/5, 1/5, 1/5], [7/20, -3/5, 3/20]]))
+      // pseudoinverses
+      assert.deepStrictEqual(inv([[1, 0], [1, 0]]), [[1/2, 1/2], [0, 0]])
+      assert.deepStrictEqual(
+         inv([[1, 0], [0, 1], [0, 1]]), [[1, 0, 0], [0, 1/2, 1/2]])
+      assert.deepStrictEqual(
+         inv([[1, 0, 0, 0, 2],
+              [0, 0, 3, 0, 0],
+              [0, 0, 0, 0, 0],
+              [0, 4, 0, 0, 0]]),
+         [[1/5,   0, 0,   0],
+          [  0,   0, 0, 1/4],
+          [  0, 1/3, 0,   0],
+          [  0,   0, 0,   0],
+          [2/5,   0, 0,   0]])
+   })
+   it('multiplies vectors and matrices', () => {
+      const mult = math.multiply
+      const pyth = [3, 4, 5]
+      assert.deepStrictEqual(mult(pyth, 2), [6, 8, 10])
+      assert.deepStrictEqual(mult(-3, pyth), [-9, -12, -15])
+      assert.strictEqual(mult(pyth, pyth), 50)
+      const mat23 = [[1, 2, 3], [-3, -2, -1]]
+      assert.deepStrictEqual(mult(mat23, pyth), [26, -22])
+      const mat32 = math.transpose(mat23)
+      assert.deepStrictEqual(mult(pyth, mat32), [26, -22])
+      assert.deepStrictEqual(mult(mat23, mat32), [[14, -10], [-10, 14]])
+      assert.deepStrictEqual(
+         mult(mat32, [[1, 2], [3, 4]]),
+         [[-8, -10], [-4, -4], [0, 2]])
+      assert(math.equal(
+         mult([[3, 0, 2], [2, 1, 0], [1, 4, 2]],
+              [[1/10, 2/5, -1/10], [-1/5, 1/5, 1/5], [7/20, -3/5, 3/20]]),
+         [[1, 0, 0], [0, 1, 0], [0, 0, 1]]))
+   })
+   it('negates a vector', () => {
+      assert.deepStrictEqual(math.negate([-3, 4, -5]), [3, -4, 5])
+   })
+   it('subtracts vectors', () => {
+      assert.deepStrictEqual(math.subtract([-3, 4, -5], [8, 0, 2]), [-11, 4, -7])
+   })
+   it('computes the sum of a vector', () => {
+      const sum = math.sum
+      assert.strictEqual(sum([-3, 4, -5]), -4)
+      assert.deepStrictEqual(sum(cV), cplx(4, -3))
+      assert.strictEqual(sum([4, true, -2]), 3)
+   })
+})

src/vector/__test__/relational.spec.js

@@ -0,0 +1,47 @@
+import assert from 'assert'
+import math from '#nanomath'
+
+const t = true
+const f = false
+describe('Vector relational functions', () => {
+   it('can test two vectors for deep equality', () => {
+      const dEq = math.deepEqual
+      const pyth = [3, 4, 5]
+      assert.strictEqual(dEq(pyth, [3, 4, 5]), t)
+      assert.strictEqual(dEq(pyth, [3, 4, 6]), f)
+      assert.strictEqual(dEq(pyth, [3, 4, [5]]), f)
+      assert.strictEqual(dEq(3, 3 + 1e-13), t)
+      assert.strictEqual(dEq(pyth, [3+1e-13, 4-1e-13, 5]), t)
+      assert.strictEqual(dEq([pyth, pyth], [pyth, [3, 4, 5]]), t)
+      assert.strictEqual(dEq([3], 3), f)
+   })
+   it('performs equal elementwise', () => {
+      const eq = math.equal
+      const pyth = [3, 4, 5]
+      assert.deepStrictEqual(eq(pyth, 4), [f, t, f])
+      assert.deepStrictEqual(eq(5, pyth), [f, f, t])
+      assert.deepStrictEqual(eq(pyth, pyth), [t, t, t])
+      assert.deepStrictEqual(eq(pyth, [3]), [t, f, f])
+      assert.deepStrictEqual(eq([4 + 1e-14], pyth), [f, t, f])
+      assert.deepStrictEqual(eq(pyth, [3, 4]), [t, t, f])
+      assert.deepStrictEqual(eq([3, 2], pyth), [t, f, f])
+      assert.deepStrictEqual(eq([pyth, pyth], [3, 4]), [[t, f, f], [f, t, f]])
+      assert.deepStrictEqual(
+         eq([pyth, pyth], [[5], pyth]), [[f, f, t], [t, t, t]])
+      assert.deepStrictEqual(
+         eq([[1, 2], [3, 4]], [[1, 3], [2, 4]]), [[t, f], [f, t]])
+   })
+   it('performs order comparisons elementwise', () => {
+      const pyth = [3, 4, 5]
+      const ans = [-1, 0, 1]
+      const powers = [2, 4, 8]
+      assert.deepStrictEqual(math.compare(pyth, [5, 4, 3]), ans)
+      assert.deepStrictEqual(math.sign([-Infinity, 0, 3]), ans)
+      assert.deepStrictEqual(math.unequal(pyth, [5, 4 + 1e-14, 5]), [t, f, f])
+      assert.deepStrictEqual(math.larger(pyth, powers), [t, f, f])
+      assert.deepStrictEqual(math.isPositive(ans), [f, f, t])
+      assert.deepStrictEqual(math.largerEq(pyth, powers), [t, t, f])
+      assert.deepStrictEqual(math.smaller(pyth, powers), [f, f, t])
+      assert.deepStrictEqual(math.smallerEq(pyth, powers), [f, t, t])
+   })
+})

src/vector/__test__/type.spec.js

@@ -0,0 +1,45 @@
+import assert from 'assert'
+import {ReturnType, Unknown} from '#core/Type.js'
+import math from '#nanomath'
+
+describe('Vector type functions', () => {
+   it('can construct a vector', () => {
+      const vec = math.vector
+      const {BooleanT, NumberT, Vector} = math.types
+      assert.deepStrictEqual(vec(3, 4, 5), [3, 4, 5])
+      assert.strictEqual(
+         ReturnType(vec.resolve([NumberT, NumberT, NumberT])),
+         Vector(NumberT))
+      assert.deepStrictEqual(vec(3, true), [3, true])
+      assert.strictEqual(
+         ReturnType(vec.resolve([NumberT, BooleanT])),
+         Vector(Unknown))
+   })
+   it('can transpose vectors and matrices', () => {
+      const tsp = math.transpose
+      assert.deepStrictEqual(tsp([3, 4, 5]), [[3], [4], [5]])
+      assert.deepStrictEqual(tsp([[1, 2], [3, 4]]), [[1, 3], [2, 4]])
+      assert.deepStrictEqual(
+         tsp([[1, 2, 3], [4, 5, 6]]), [[1, 4], [2, 5], [3, 6]])
+   })
+   it('can take adjoint (conjugate transpose) of a matrix', () => {
+      const cx = math.complex
+      assert.deepStrictEqual(
+         math.adjoint([[cx(1, 1), cx(2, 2)], [cx(3, 3), cx(4, 4)]]),
+         [[cx(1, -1), cx(3, -3)], [cx(2, -2), cx(4, -4)]])
+   })
+   it('generates identity from an example matrix or a number of rows', () => {
+      const id = math.identity
+      const cx = math.complex
+      assert.deepStrictEqual(id(2), [[1, 0], [0, 1]])
+      assert.deepStrictEqual(id(cx(2)), [[cx(1), cx(0)], [cx(0), cx(1)]])
+      assert.deepStrictEqual(
+         id([[1, 2, 3]]),
+         [[1, 0 , 0], [0, 1, 0], [0, 0, 1]])
+   })
+   it('takes the determinant of a matrix', () => {
+      assert.strictEqual(
+         math.determinant([[6, 1, 1], [4, -2, 5], [2, 8, 7]]),
+         -306)
+   })
+})

src/vector/__test__/utils.spec.js

@@ -0,0 +1,41 @@
+import assert from 'assert'
+import math from '#nanomath'
+
+describe('Vector utility functions', () => {
+   it('can clone a vector', () => {
+      const subj1 = [3, 4]
+      const clone1 = math.clone(subj1)
+      assert.deepStrictEqual(clone1, subj1)
+      assert(clone1 !== subj1)
+      const subj2 = [3, false]
+      const clone2 = math.clone(subj2)
+      assert.deepStrictEqual(clone2, subj2)
+      assert(clone2 !== subj2)
+      const subj3 = [[3, 4], [2, 5]]
+      const clone3 = math.clone(subj3)
+      assert.deepStrictEqual(clone3, subj3)
+      assert(clone3 !== subj3)
+      assert(clone3[0] !== subj3[0])
+      assert(clone3[1] !== subj3[1])
+   })
+   it('performs utilities elementwise', () => {
+      const cplx = math.complex
+      const sink = math.vector(
+         NaN, -Infinity, 7,
+         cplx(3, 4), cplx(3.5, -2.5), cplx(2.2), cplx(cplx(3, 4))
+      )
+      const t = true, f = false
+      assert.deepStrictEqual(math.isnan(sink),        [t, f, f, f, f, f, f])
+      assert.deepStrictEqual(math.isfinite(sink),     [f, f, t, t, t, t, t])
+      assert.deepStrictEqual(math.isInteger(sink),    [f, f, t, t, f, f, t])
+      assert.deepStrictEqual(math.isReal(sink),       [t, t, t, f, f, t, f])
+      assert.deepStrictEqual(math.nonImaginary(sink), [t, t, t, f, f, t, t])
+      assert.deepStrictEqual(math.re(sink), [NaN, -Infinity, 7, 3, 3.5, 2.2, 3])
+      assert.deepStrictEqual(
+         math.im(sink),
+         [0, 0, 0, cplx(0, 4), cplx(0, -2.5), cplx(0, 0), cplx(cplx(0, 4))])
+   })
+})
+
+      
+      

src/vector/all.js

@@ -0,0 +1,7 @@
+export * as typeDefinition from './Vector.js'
+export * as arithmetic from './arithmetic.js'
+export * as logical from './logical.js'
+export * as relational from './relational.js'
+export * as type from './type.js'
+export * as utilities from './utils.js'
+

src/vector/arithmetic.js

@@ -0,0 +1,266 @@
+import {
+   distributeFirst, distributeSecond, promoteBinary, promoteUnary
+} from './helpers.js'
+import {Vector} from './Vector.js'
+
+import {Returns, ReturnType} from '#core/Type.js'
+import {match} from '#core/TypePatterns.js'
+import {ReturnsAs} from '#generic/helpers.js'
+
+export const normsq = match(Vector, (math, V) => {
+   const compNormsq = math.normsq.resolve(V.Component)
+   const sum = math.sum.resolve(Vector(ReturnType(compNormsq)))
+   return ReturnsAs(sum, v => sum(v.map(compNormsq)))
+})
+// abs and norm differ only on Vector (and perhaps other collections) --
+// norm computes overall by the generic formula, whereas abs distributes
+// elementwise:
+export const abs = promoteUnary('abs')
+export const add = promoteBinary('add')
+
+export const dotMultiply = promoteBinary('multiply')
+export const invert = match(Vector, (math, V, strategy) => {
+   if (V.vectorDepth > 2) {
+      throw new TypeError(
+         'invert not implemented for arrays of dimension > 2')
+   }
+   const normsq = math.normsq.resolve(V, strategy)
+   const NormT = ReturnType(normsq)
+   const zNorm = math.zero(NormT)
+   const isRealZ = math.isZero.resolve(NormT, strategy)
+   if (V.vectorDepth === 1) {
+      const invNorm = math.invert.resolve(NormT, strategy)
+      const scalarMult = math.multiply.resolve(
+         [V, ReturnType(invNorm)], strategy)
+      return ReturnsAs(scalarMult, v => {
+         const nsq = normsq(v)
+         if (isRealZ(nsq)) return Array(v.length).fill(zNorm)
+         return scalarMult(v, invNorm(nsq))
+      })
+   }
+   // usual matrix situation, want to find a matrix whose product with v
+   // is the identity, or as close as we can get to that if the rank is
+   // deficient. We use the Moore-Penrose pseudoinverse.
+   const clone = math.clone.resolve(V, strategy)
+   const VComp = V.Component
+   const Elt = VComp.Component
+   const invElt = math.invert.resolve(Elt, strategy)
+   const det = math.determinant.resolve(V, strategy)
+   const neg = math.negate.resolve(Elt, strategy)
+   const multMM = math.multiply.resolve([V, V], strategy)
+   const multMS = math.multiply.resolve([V, ReturnType(invElt)], strategy)
+   const multVS = math.multiply.resolve([VComp, ReturnType(invElt)], strategy)
+   const multSS = math.multiply.resolve([Elt, Elt], strategy)
+   const sub = math.subtract.resolve([Elt, Elt], strategy)
+   const id = math.identity.resolve(V, strategy)
+   const abs = math.abs.resolve(Elt, strategy)
+   const gt = math.larger.resolve([NormT, NormT], strategy)
+   const isEltZ = math.isZero.resolve(Elt, strategy)
+   const zElt = math.zero(Elt)
+   const adj = math.adjoint.resolve(V, strategy)
+   const methods = {
+      invElt, det, isRealZ, neg, multMM, multMS, multVS, multSS,
+      id, abs, isEltZ, zElt, gt, sub, adj, clone
+   }
+   if (ReturnType(abs) !== NormT) {
+      throw new TypeError('type inconsistency in matrix invert')
+   }
+   return Returns(V, m => {
+      const rows = m.length
+      const cols = m[0].length
+      const nsq = normsq(m)
+      if (isRealZ(nsq)) { // all-zero matrix
+         const retval = []
+         for (let ix = 0; ix < cols; ++ix) {
+            retval.push(Array(rows).fill(zNorm))
+         }
+         return retval
+      }
+      if (rows == cols) {
+         // the inv helper will return falsy if not invertible
+         const retval = inv(m, rows, methods)
+         if (retval) return retval
+      }
+
+      return pinv(m, rows, cols, methods)
+   })
+})
+
+// Returns the inverse of a rows×rows matrix or false if not invertible
+// Note: destroys m in the inversion process.
+function inv(
+   origm,
+   rows,
+   {
+      invElt, det, isRealZ, neg, multMS, multVS, multSS,
+      id, abs, isEltZ, zElt, gt, sub, clone
+   }
+) {
+   switch (rows) {
+   case 1: return [[invElt(origm[0][0])]]
+   case 2: {
+      const dt = det(origm)
+      if (isRealZ(dt)) return false
+      const divisor = invElt(dt)
+      const [[a, b], [c, d]] = origm
+      return multMS([[d, neg(b)], [neg(c), a]], divisor)
+   }
+   default: { // Gauss-Jordan elimination
+      const m = clone(origm)
+      const B = id(m)
+      const cols = rows
+      // Loop over columns, performing row reductions:
+      for (let c = 0; c < cols; ++c) {
+         // Pivot: Find row r that has the largest entry in column c, and
+         // swap row c and r:
+         let colMax = abs(m[c][c])
+         let rMax = c
+         for (let r = c + 1; r < rows; ++r) {
+            const mag = abs(m[r][c])
+            if (gt(mag, colMax)) {
+               colMax = mag
+               rMax = r
+            }
+         }
+         if (isRealZ(colMax)) return false
+         if (rMax !== c) {
+            [m[c], m[rMax], B[c], B[rMax]]
+               = [m[rMax], m[c], B[rMax], B[c]]
+         }
+         // Normalize the cth row:
+         const normalizer = invElt(m[c][c])
+         const mc = multVS(m[c], normalizer)
+         m[c] = mc
+         const Bc = multVS(B[c], normalizer)
+         B[c] = Bc
+
+         // Eliminate nonzero values on other rows at column c
+         for (let r = 0; r < rows; ++r) {
+            if (r === c) continue
+            const mr = m[r]
+            const Br = B[r]
+            const mrc = mr[c]
+            if (!isEltZ(mr[c])) {
+               // Subtract Arc times row c from row r to eliminate A[r][c]
+               mr[c] = zElt
+               for (let s = c + 1; s < cols; ++s) {
+                  mr[s] = sub(mr[s], multSS(mrc, mc[s]))
+               }
+               for (let s = 0; s < cols; ++s) {
+                  Br[s] = sub(Br[s], multSS(mrc, Bc[s]))
+               }
+            }
+         }
+      }
+
+      return B
+   }}
+}
+
+// Calculates Moore-Penrose pseudoinverse
+// uses rank factorization per mathjs; SVD appears to be considered better
+// but not worth the effort to implement for this prototype
+function pinv(m, rows, cols, methods) {
+   const {C, F} = rankFactorization(m, rows, cols, methods)
+   const {multMM, adj} = methods
+   const Cstar = adj(C)
+   const Cpinv = multMM(inv(multMM(Cstar, C), Cstar.length, methods), Cstar)
+   const Fstar = adj(F)
+   const Fpinv = multMM(Fstar, inv(multMM(F, Fstar), F.length, methods))
+   return multMM(Fpinv, Cpinv)
+}
+
+// warning: destroys m in computing the row-reduced echelon form
+// TODO: this code should be merged with inv to the extent possible. It's
+// a very similar process.
+function rref(origm, rows, cols, methods) {
+   const {isEltZ, invElt, multVS, zElt, sub, multSS, clone} = methods
+   const m = clone(origm)
+   let lead = -1
+   for (let r = 0; r < rows && ++lead < cols; ++r) {
+      if (cols <= lead) return m
+      let i = r
+      while (isEltZ(m[i][lead])) {
+         if (++i === rows) {
+            i = r
+            if (++lead === cols) return m
+         }
+      }
+
+      if (i !== r) [m[i], m[r]] = [m[r], m[i]]
+
+      let normalizer = invElt(m[r][lead])
+      const mr = multVS(m[r], normalizer)
+      m[r] = mr
+      for (let i = 0; i < rows; ++i) {
+         if (i === r) continue
+         const mi = m[i]
+         const toRemove = mi[lead]
+         mi[lead] = zElt
+         for (let j = lead + 1; j < cols; ++j) {
+            mi[j] = sub(mi[j], multSS(toRemove, mr[j]))
+         }
+      }
+   }
+   return m
+}
+
+function rankFactorization(m, rows, cols, methods) {
+   const RREF = rref(m, rows, cols, methods)
+   const {isEltZ} = methods
+   const rankRows = RREF.map(row => row.some(elt => !isEltZ(elt)))
+   const C = m.map(row => row.filter((_, j) => j < rows && rankRows[j]))
+   const F = RREF.filter((_, i) => rankRows[i])
+   return {C, F}
+}
+
+export const multiply = [
+   distributeFirst('multiply'),
+   distributeSecond('multiply'),
+   match([Vector, Vector], (math, [V, W], strategy) => {
+      const VComp = V.Component
+      if (W.vectorDepth === 1) {
+         if (V.vectorDepth === 1) {
+            const eltWise = math.dotMultiply.resolve([V, W], strategy)
+            const sum = math.sum.resolve(ReturnType(eltWise))
+            return ReturnsAs(sum, (v, w) => sum(eltWise(v, w)))
+         }
+         const compMult = math.multiply.resolve([VComp, W], strategy)
+         return ReturnsAs(
+            Vector(ReturnType(compMult)),
+            (v, w) => v.map(f => compMult(f, w)))
+      }
+      const transpose = math.transpose.resolve(W, strategy)
+      const wrapV = V.vectorDepth === 1
+      const RowV = wrapV ? V : VComp
+      const rowMult = math.multiply.resolve([RowV, W.Component], strategy)
+      let RetType = Vector(ReturnType(rowMult))
+      if (!wrapV) RetType = Vector(RetType)
+      return ReturnsAs(RetType, (v, w) => {
+         if (wrapV) v = [v]
+         w = transpose(w)
+         let retval = v.map(vrow => w.map(wcol => rowMult(vrow, wcol)))
+         if (wrapV) retval = retval[0]
+         return retval
+      })
+   })
+]
+
+export const negate = promoteUnary('negate')
+export const subtract = promoteBinary('subtract')
+
+export const sum = match(Vector, (math, V) => {
+   const add = math.add.resolve([V.Component, V.Component])
+   const haszero = math.haszero(V.Component)
+   const zero = haszero ? math.zero(V.Component) : undefined
+   return ReturnsAs(add, v => {
+      if (v.length === 0) {
+         if (haszero) return zero
+         throw new TypeError(`Can't sum empty ${V}: no zero element`)
+      }
+      let ix = 0
+      let retval = v[ix]
+      while (++ix < v.length) retval = add(retval, v[ix])
+      return retval
+   })
+})

src/vector/helpers.js

@@ -0,0 +1,70 @@
+import {Vector} from './Vector.js'
+import {Returns, ReturnType, Undefined} from '#core/Type.js'
+import {Any, match} from '#core/TypePatterns.js'
+
+export const promoteUnary = name => match(Vector, (math, V, strategy) => {
+   const compOp = math.resolve(name, V.Component, strategy)
+   return Returns(Vector(ReturnType(compOp)), v => v.map(elt => compOp(elt)))
+})
+
+export const distributeFirst = name => match(
+   [Vector, Any],
+   (math, [V, E], strategy) => {
+      const compOp = math.resolve(name, [V.Component, E], strategy)
+      return Returns(
+         Vector(ReturnType(compOp)), (v, e) => v.map(f => compOp(f, e)))
+   })
+
+export const distributeSecond = name => match(
+   [Any, Vector],
+   (math, [E, V], strategy) => {
+      const compOp = math.resolve(name, [E, V.Component], strategy)
+      return Returns(
+         Vector(ReturnType(compOp)), (e, v) => v.map(f => compOp(e, f)))
+   })
+
+export const promoteBinary = name => [
+   distributeFirst(name),
+   distributeSecond(name),
+   match([Vector, Vector], (math, [V, W], strategy) => {
+      const VComp = V.Component
+      const WComp = W.Component
+      // special case: if the vector nesting depths do not match,
+      // we operate between the elements of the deeper one and the entire
+      // more shallow one:
+      if (V.vectorDepth > W.vectorDepth) {
+         const compOp = math.resolve(name, [VComp, W], strategy)
+         return Returns(
+            Vector(ReturnType(compOp)), (v, w) => v.map(f => compOp(f, w)))
+      }
+      if (V.vectorDepth < W.vectorDepth) {
+         const compOp = math.resolve(name, [V, WComp], strategy)
+         return Returns(
+            Vector(ReturnType(compOp)), (v, w) => w.map(f => compOp(v, f)))
+      }
+      const compOp = math.resolve(name, [VComp, WComp], strategy)
+      const opNoV = math.resolve(name, [Undefined, WComp], strategy)
+      const opNoW = math.resolve(name, [VComp, Undefined], strategy)
+      return Returns(
+         Vector(ReturnType(compOp)),
+         (v, w) => {
+            const vInc = Number(v.length > 1)
+            const wInc = Number(w.length >= v.length || w.length > 1)
+            const retval = []
+            let vIx = 0
+            let wIx = 0
+            while ((vInc && vIx < v.length)
+                   || (wInc && wIx < w.length)
+            ) {
+               if (vIx >= v.length) {
+                  retval.push(opNoV(undefined, w[wIx]))
+               } else if (wIx >= w.length) {
+                  retval.push(opNoW(v[vIx], undefined))
+               } else retval.push(compOp(v[vIx], w[wIx]))
+               vIx += vInc
+               wIx += wInc
+            }
+            return retval
+         })
+   })
+]

src/vector/logical.js

@@ -0,0 +1,7 @@
+import {promoteBinary, promoteUnary} from './helpers.js'
+
+export const not = promoteUnary('not')
+export const and = promoteBinary('and')
+export const or = promoteBinary('or')
+
+

src/vector/relational.js

@@ -0,0 +1,95 @@
+import {Vector} from './Vector.js'
+import {promoteBinary} from './helpers.js'
+
+import {Unknown, Returns, ReturnType, Undefined} from '#core/Type.js'
+import {Any, Optional, match} from '#core/TypePatterns.js'
+import {BooleanT} from '#boolean/BooleanT.js'
+
+export const deepEqual = [
+   match([Vector, Any], Returns(BooleanT, () => false)),
+   match([Any, Vector], Returns(BooleanT, () => false)),
+   match([Vector, Vector], (math, [V, W]) => {
+      const compDeep = math.deepEqual.resolve([V.Component, W.Component])
+      return Returns(BooleanT, (v,w) => v === w
+         || (v.length === w.length
+             && v.every((e, i) => compDeep(e, w[i]))))
+   })
+]
+
+export const indistinguishable = [
+   match([Vector, Any, Optional([Any, Any])], (math, [V, E, T]) => {
+      const VComp = V.Component
+      if (T.length === 0) { // no tolerances
+         const same = math.indistinguishable.resolve([VComp, E])
+         return Returns(
+            Vector(ReturnType(same)), (v, e) => v.map(f => same(f, e)))
+      }
+      const [[RT, AT]] = T
+      const same = math.indistinguishable.resolve([VComp, E, RT, AT])
+      return Returns(
+         Vector(ReturnType(same)),
+         (v, e, [[rT, aT]]) => v.map(f => same(f, e, rT, aT)))
+   }),
+   match([Any, Vector, Optional([Any, Any])], (math, [E, V, T]) => {
+      // reimplement to get other order in same so as not to assume
+      // same is symmetric, even though it probably is
+      const VComp = V.Component
+      if (T.length === 0) { // no tolerances
+         const same = math.indistinguishable.resolve([E, VComp])
+         return Returns(
+            Vector(ReturnType(same)), (e, v) => v.map(f => same(e, f)))
+      }
+      const [[RT, AT]] = T
+      const same = math.indistinguishable.resolve([E, VComp, RT, AT])
+      return Returns(
+         Vector(ReturnType(same)),
+         (e, v, [[rT, aT]]) => v.map(f => same(e, f, rT, aT)))
+   }),
+   match([Vector, Vector, Optional([Any, Any])], (math, [V, W, T]) => {
+      const VComp = V.Component
+      const WComp = W.Component
+      const inhomogeneous = VComp === Unknown || WComp === Unknown
+      let same
+      let sameNoV
+      let sameNoW
+      if (inhomogeneous) {
+         same = math.indistinguishable
+         sameNoV = same
+         sameNoW = same
+      } else if (T.length === 0) { // no tolerances
+         same = math.indistinguishable.resolve([VComp, WComp])
+         sameNoV = math.indistinguishable.resolve([Undefined, WComp])
+         sameNoW = math.indistinguishable.resolve([VComp, Undefined])
+      } else {
+         const [[RT, AT]] = T
+         same = math.indistinguishable.resolve([VComp, WComp, RT, AT])
+         sameNoV = math.indistinguishable.resolve([Undefined, WComp, RT, AT])
+         sameNoW = math.indistinguishable.resolve([VComp, Undefined, RT, AT])
+      }
+      return Returns(
+         inhomogeneous ? Vector(Unknown) : Vector(ReturnType(same)),
+         (v, w, [tol = [0, 0]]) => {
+            const [rT, aT] = tol
+            const vInc = Number(v.length > 1)
+            const wInc = Number(w.length >= v.length || w.length > 1)
+            const retval = []
+            let vIx = 0
+            let wIx = 0
+            while ((vInc && vIx < v.length)
+                   || (wInc && wIx < w.length)
+            ) {
+               if (vIx >= v.length) {
+                  retval.push(sameNoV(undefined, w[wIx], rT, aT))
+               } else if (wIx >= w.length) {
+                  retval.push(sameNoW(v[vIx], undefined, rT, aT))
+               } else retval.push(same(v[vIx], w[wIx], rT, aT))
+               vIx += vInc
+               wIx += wInc
+            }
+            return retval
+         })
+   })
+]
+
+export const compare = promoteBinary('compare')
+export const exceeds = promoteBinary('exceeds')

src/vector/type.js

@@ -0,0 +1,133 @@
+import {Vector} from './Vector.js'
+import {OneOf, Returns, ReturnType, Unknown} from '#core/Type.js'
+import {Any, Multiple, match} from '#core/TypePatterns.js'
+
+export const vector = match(Multiple(Any), (math, [TV]) => {
+   if (!TV.length) return Returns(Vector(Unknown), () => [])
+   let CompType = TV[0]
+   if (TV.some(T => T !== CompType)) CompType = Unknown
+   return Returns(Vector(CompType), v => v)
+})
+
+export const determinant = match(Vector(Vector), (math, M, strategy) => {
+   const Elt = M.Component.Component
+   const cloneElt = math.clone.resolve(Elt, strategy)
+   const mult = math.multiply.resolve([Elt, Elt], strategy)
+   const sub = math.subtract.resolve(
+      [ReturnType(mult), ReturnType(mult)], strategy)
+   const isZ = math.isZero.resolve(Elt, strategy)
+   const zElt = math.zero(Elt)
+   const clone = math.clone.resolve(M, strategy)
+   const div = math.divide.resolve([ReturnType(sub), Elt], strategy)
+   const neg = math.negate.resolve(Elt, strategy)
+   return Returns(OneOf(ReturnType(clone), ReturnType(sub), Elt), origm => {
+      const rows = origm.length
+      switch (rows) {
+      case 1: return cloneElt(origm[0][0])
+      case 2: {
+         const [[a, b], [c, d]] = origm
+         return sub(mult(a, d), mult(b, c))
+      }
+      default: { // Bareiss algorithm
+         const m = clone(origm)
+         let negated = false
+         const rowIndices = [...Array(rows).keys()] // track row indices
+         // because the algorithm may swap rows
+         for (let k = 0; k < rows; ++k) {
+            let k_ = rowIndices[k]
+            if (isZ(m[k_][k])) {
+               let _k = k
+               while (++_k < rows) {
+                  if (!isZ(m[rowIndices[_k]][k])) {
+                     k_ = rowIndices[_k]
+                     rowIndices[_k] = rowIndices[k]
+                     rowIndices[k] = k_
+                     negated = !negated
+                     break
+                  }
+               }
+               if (_k === rows) return zElt
+            }
+            const piv = m[k_][k] // we now know nonzero
+            const piv_ = k === 0 ? 1 : m[rowIndices[k-1]][k-1]
+            for (let i = k + 1; i < rows; ++i) {
+               const i_ = rowIndices[i]
+               for (let j = k + 1; j < rows; ++j) {
+                  m[i_][j] = div(
+                     sub(mult(m[i_][j], piv), mult(m[i_][k], m[k_][j])),
+                     piv_)
+               }
+            }
+         }
+         const det = m[rowIndices[rows - 1]][rows - 1]
+         return negated ? neg(det) : det
+      }}
+   })
+})
+
+function identitizer(cols, zero, one) {
+   const retval = []
+   for (let ix = 0; ix < cols; ++ix) {
+      const row = Array(cols).fill(zero)
+      row[ix] = one
+      retval.push(row)
+   }
+   return retval
+}
+
+export const identity = [
+   match(Any, (math, V) => {
+      const toNum = math.number.resolve(V)
+      const zero = math.zero(V)
+      const one = math.one(V)
+      return Returns(Vector(Vector(V)), n => identitizer(toNum(n), zero, one))
+   }),
+   match(Vector, (math, V) => {
+      switch (V.vectorDepth) {
+      case 1: {
+         const Elt = V.Component
+         const one = math.one(Elt)
+         return Returns(math.typeOf(one), () => one)
+      }
+      case 2: {
+         const Elt = V.Component.Component
+         const zero = math.zero(Elt)
+         const one = math.one(Elt)
+         return Returns(V, m => identitizer(
+            m.length ? m[0].length : 0, zero, one))
+      }
+      default:
+         throw new RangeError(
+            `'identity' not implemented on ${V.vectorDepth} dimensional arrays`)
+      }
+   })
+]
+
+// transposes a 2D matrix
+function transposer(wrap, eltFun) {
+   return v => {
+      if (wrap) v = [v]
+      const cols = v.length ? v[0].length : 0
+      const retval = []
+      for (let ix = 0; ix < cols; ++ix) {
+         retval.push(v.map(row => eltFun(row[ix])))
+      }
+      return retval
+   }
+}
+
+export const transpose = match(Vector, (_math, V) => {
+   const wrapV = V.vectorDepth === 1
+   const Mat = wrapV ? Vector(V) : V
+   return Returns(Mat, transposer(wrapV, elt => elt))
+})
+
+// or with conjugation:
+export const adjoint = match(Vector, (math, V, strategy) => {
+   const wrapV = V.vectorDepth === 1
+   const VComp = V.Component
+   const Elt = wrapV ? VComp : VComp.Component
+   const conj = math.conj.resolve(Elt, strategy)
+   const Mat = Vector(Vector(ReturnType(conj)))
+   return Returns(Mat, transposer(wrapV, conj))
+})

src/vector/utils.js

@@ -0,0 +1,10 @@
+import {promoteUnary} from './helpers.js'
+
+export const clone = promoteUnary('clone')
+export const isnan = promoteUnary('isnan')
+export const isfinite = promoteUnary('isfinite')
+export const isInteger = promoteUnary('isInteger')
+export const isReal = promoteUnary('isReal')
+export const nonImaginary = promoteUnary('nonImaginary')
+export const re = promoteUnary('re')
+export const im = promoteUnary('im')