From 8bf167040856b3ae624dbab3506e4c7f4026a78e Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Thu, 13 Feb 2020 12:37:33 +0100
Subject: [PATCH 01/19] split the eigs function into multiple algorithms

---
 src/function/matrix/eigs.js              | 392 ++++-------------------
 src/function/matrix/eigs/complex.js      |  14 +
 src/function/matrix/eigs/realSymetric.js | 306 ++++++++++++++++++
 3 files changed, 380 insertions(+), 332 deletions(-)
 create mode 100644 src/function/matrix/eigs/complex.js
 create mode 100644 src/function/matrix/eigs/realSymetric.js

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 304da41e18..a09d60bf40 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -1,23 +1,34 @@
-import { clone } from '../../utils/object'
 import { factory } from '../../utils/factory'
 import { format } from '../../utils/string'
+import { createRealSymmetric } from './eigs/realSymetric'
+import { createComplex } from './eigs/complex'
 
 const name = 'eigs'
-const dependencies = ['typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'inv', 'bignumber', 'multiply', 'add']
-
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) => {
+const dependencies = [
+  'typed',
+  'matrix',
+  'bignumber',
+  'add',
+  'addScalar',
+  'subtract',
+  'multiply',
+  'multiplyScalar',
+  'abs',
+  'equal',
+  'larger',
+  'atan',
+  'cos',
+  'sin',
+  'inv'
+]
+
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) => {
   /**
-   * Compute eigenvalue and eigenvector of a real symmetric matrix.
-   * Only applicable to two dimensional symmetric matrices. Uses Jacobi
-   * Algorithm. Matrix containing mixed type ('number', 'bignumber', 'fraction')
-   * of elements are not supported. Input matrix or 2D array should contain all elements
-   * of either 'number', 'bignumber' or 'fraction' type. For 'number' and 'fraction', the
-   * eigenvalues are of 'number' type. For 'bignumber' the eigenvalues are of ''bignumber' type.
-   * Eigenvectors are always of 'number' type.
+   * Compute eigenvalues and eigenvectors of a matrix.
    *
    * Syntax:
    *
-   *     math.eigs(x)
+   *     math.eigs(x, [prec])
    *
    * Examples:
    *
@@ -27,352 +38,69 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
    *     const U = ans.vectors
    *     const UTxHxU = math.multiply(math.transpose(U), H, U) // rotates H to the eigen-representation
    *     E[0] == UTxHxU[0][0]  // returns true
+   *
    * See also:
    *
    *     inv
    *
    * @param {Array | Matrix} x  Matrix to be diagonalized
+   * @param {number | BigNumber} [prec] Precision, default value: 1e-15
    * @return {{values: Array, vectors: Array} | {values: Matrix, vectors: Matrix}} Object containing eigenvalues (Array or Matrix) and eigenvectors (2D Array/Matrix).
    */
   const eigs = typed('eigs', {
 
-    Array: function (x) {
-      // check array size
+    'Array': function (x) {
       const mat = matrix(x)
-      const size = mat.size()
-      if (size.length !== 2 || size[0] !== size[1]) {
-        throw new RangeError('Matrix must be square ' +
-          '(size: ' + format(size) + ')')
-      }
-
-      // use dense 2D matrix implementation
-      const ans = checkAndSubmit(mat, size[0])
-      return { values: ans[0], vectors: ans[1] }
+      return computeValuesAndVectors(mat);
     },
 
-    Matrix: function (x) {
-      // use dense 2D array implementation
-      // dense matrix
-      const size = x.size()
-      if (size.length !== 2 || size[0] !== size[1]) {
-        throw new RangeError('Matrix must be square ' +
-          '(size: ' + format(size) + ')')
-      }
-      const ans = checkAndSubmit(x, size[0])
-      return { values: matrix(ans[0]), vectors: matrix(ans[1]) }
-    }
-  })
-
-  // Is the matrix
-  // symmetric ?
-  function isSymmetric (x, n) {
-    for (let i = 0; i < n; i++) {
-      for (let j = i; j < n; j++) {
-        // not symmtric
-        if (!equal(x[i][j], x[j][i])) {
-          throw new TypeError('Input matrix is not symmetric')
-        }
-      }
-    }
-  }
-
-  // check input for possible problems
-  // and perform diagonalization efficiently for
-  // specific type of number
-  function checkAndSubmit (x, n) {
-    let type = x.datatype()
-    // type check
-    if (type === undefined) {
-      type = x.getDataType()
-    }
-    if (type !== 'number' && type !== 'BigNumber' && type !== 'Fraction') {
-      if (type === 'mixed') {
-        throw new TypeError('Mixed matrix element type is not supported')
-      } else {
-        throw new TypeError('Matrix element type not supported (' + type + ')')
-      }
-    } else {
-      isSymmetric(x.toArray(), n)
-    }
-
-    // perform efficient calculation for 'numbers'
-    if (type === 'number') {
-      return diag(x.toArray())
-    } else if (type === 'Fraction') {
-      const xArr = x.toArray()
-      // convert fraction to numbers
-      for (let i = 0; i < n; i++) {
-        for (let j = i; j < n; j++) {
-          xArr[i][j] = xArr[i][j].valueOf()
-          xArr[j][i] = xArr[i][j]
-        }
-      }
-      return diag(x.toArray())
-    } else if (type === 'BigNumber') {
-      return diagBig(x.toArray())
-    }
-  }
-
-  // diagonalization implementation for number (efficient)
-  function diag (x, precision = 1E-12) {
-    const N = x.length
-    const e0 = Math.abs(precision / N)
-    let psi
-    let Sij = new Array(N)
-    // Sij is Identity Matrix
-    for (let i = 0; i < N; i++) {
-      Sij[i] = createArray(N, 0)
-      Sij[i][i] = 1.0
-    }
-    // initial error
-    let Vab = getAij(x)
-    while (Math.abs(Vab[1]) >= Math.abs(e0)) {
-      const i = Vab[0][0]
-      const j = Vab[0][1]
-      psi = getTheta(x[i][i], x[j][j], x[i][j])
-      x = x1(x, psi, i, j)
-      Sij = Sij1(Sij, psi, i, j)
-      Vab = getAij(x)
-    }
-    const Ei = createArray(N, 0) // eigenvalues
-    for (let i = 0; i < N; i++) {
-      Ei[i] = x[i][i]
-    }
-    return sorting(clone(Ei), clone(Sij))
-  }
+    'Array, number|BigNumber': function (x, prec) {
+      const mat = matrix(x)
+      return computeValuesAndVectors(mat, prec);
+    },
 
-  // diagonalization implementation for bigNumber
-  function diagBig (x, precision = 1E-12) {
-    const N = x.length
-    const e0 = abs(precision / N)
-    let psi
-    let Sij = new Array(N)
-    // Sij is Identity Matrix
-    for (let i = 0; i < N; i++) {
-      Sij[i] = createArray(N, 0)
-      Sij[i][i] = 1.0
-    }
-    // initial error
-    let Vab = getAijBig(x)
-    while (abs(Vab[1]) >= abs(e0)) {
-      const i = Vab[0][0]
-      const j = Vab[0][1]
-      psi = getThetaBig(x[i][i], x[j][j], x[i][j])
-      x = x1Big(x, psi, i, j)
-      Sij = Sij1Big(Sij, psi, i, j)
-      Vab = getAijBig(x)
-    }
-    const Ei = createArray(N, 0) // eigenvalues
-    for (let i = 0; i < N; i++) {
-      Ei[i] = x[i][i]
-    }
-    // return [clone(Ei), clone(Sij)]
-    return sorting(clone(Ei), clone(Sij))
-  }
+    'Matrix': function (mat) {
+      return computeValuesAndVectors(mat);
+    },
 
-  // get angle
-  function getTheta (aii, ajj, aij) {
-    let th = 0
-    const denom = (ajj - aii)
-    if (Math.abs(denom) <= 1E-14) {
-      th = Math.PI / 4.0
-    } else {
-      th = 0.5 * Math.atan(2.0 * aij / (ajj - aii))
+    'Matrix, number|BigNumber': function (mat, prec) {
+      return computeValuesAndVectors(mat, prec);
     }
-    return th
-  }
+  })
 
-  // get angle
-  function getThetaBig (aii, ajj, aij) {
-    let th = 0
-    const denom = subtract(ajj, aii)
-    if (abs(denom) <= 1E-14) {
-      th = Math.PI / 4.0
-    } else {
-      th = multiplyScalar(0.5, atan(multiply(2.0, aij, inv(denom))))
-    }
-    return th
-  }
+  const doRealSymetric = createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add });
+  const doComplex = createComplex();
 
-  // update eigvec
-  function Sij1 (Sij, theta, i, j) {
-    const N = Sij.length
-    const c = Math.cos(theta)
-    const s = Math.sin(theta)
-    const Ski = createArray(N, 0)
-    const Skj = createArray(N, 0)
-    for (let k = 0; k < N; k++) {
-      Ski[k] = c * Sij[k][i] - s * Sij[k][j]
-      Skj[k] = s * Sij[k][i] + c * Sij[k][j]
-    }
-    for (let k = 0; k < N; k++) {
-      Sij[k][i] = Ski[k]
-      Sij[k][j] = Skj[k]
-    }
-    return Sij
-  }
-  // update eigvec for overlap
-  function Sij1Big (Sij, theta, i, j) {
-    const N = Sij.length
-    const c = cos(theta)
-    const s = sin(theta)
-    const Ski = createArray(N, 0)
-    const Skj = createArray(N, 0)
-    for (let k = 0; k < N; k++) {
-      Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]))
-      Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]))
-    }
-    for (let k = 0; k < N; k++) {
-      Sij[k][i] = Ski[k]
-      Sij[k][j] = Skj[k]
-    }
-    return Sij
-  }
+  /** @return {boolean} */
+  function isSymmetric(arr, N, prec) {
 
-  // update matrix
-  function x1Big (Hij, theta, i, j) {
-    const N = Hij.length
-    const c = bignumber(cos(theta))
-    const s = bignumber(sin(theta))
-    const c2 = multiplyScalar(c, c)
-    const s2 = multiplyScalar(s, s)
-    const Aki = createArray(N, 0)
-    const Akj = createArray(N, 0)
-    // 2cs Hij
-    const csHij = multiply(2, c, s, Hij[i][j])
-    //  Aii
-    const Aii = addScalar(subtract(multiplyScalar(c2, Hij[i][i]), csHij), multiplyScalar(s2, Hij[j][j]))
-    const Ajj = add(multiplyScalar(s2, Hij[i][i]), csHij, multiplyScalar(c2, Hij[j][j]))
-    // 0  to i
-    for (let k = 0; k < N; k++) {
-      Aki[k] = subtract(multiplyScalar(c, Hij[i][k]), multiplyScalar(s, Hij[j][k]))
-      Akj[k] = addScalar(multiplyScalar(s, Hij[i][k]), multiplyScalar(c, Hij[j][k]))
-    }
-    // Modify Hij
-    Hij[i][i] = Aii
-    Hij[j][j] = Ajj
-    Hij[i][j] = 0
-    Hij[j][i] = 0
-    // 0  to i
-    for (let k = 0; k < N; k++) {
-      if (k !== i && k !== j) {
-        Hij[i][k] = Aki[k]
-        Hij[k][i] = Aki[k]
-        Hij[j][k] = Akj[k]
-        Hij[k][j] = Akj[k]
-      }
-    }
-    return Hij
-  }
+    for (let i = 0; i < N; i++)
+      for (let j = i; j < N; j++)
+        if ( larger( subtract(arr[i][j], arr[j][i]), prec) )
+          return false;
 
-  // update matrix
-  function x1 (Hij, theta, i, j) {
-    const N = Hij.length
-    const c = Math.cos(theta)
-    const s = Math.sin(theta)
-    const c2 = c * c
-    const s2 = s * s
-    const Aki = createArray(N, 0)
-    const Akj = createArray(N, 0)
-    //  Aii
-    const Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j]
-    const Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j]
-    // 0  to i
-    for (let k = 0; k < N; k++) {
-      Aki[k] = c * Hij[i][k] - s * Hij[j][k]
-      Akj[k] = s * Hij[i][k] + c * Hij[j][k]
-    }
-    // Modify Hij
-    Hij[i][i] = Aii
-    Hij[j][j] = Ajj
-    Hij[i][j] = 0
-    Hij[j][i] = 0
-    // 0  to i
-    for (let k = 0; k < N; k++) {
-      if (k !== i && k !== j) {
-        Hij[i][k] = Aki[k]
-        Hij[k][i] = Aki[k]
-        Hij[j][k] = Akj[k]
-        Hij[k][j] = Akj[k]
-      }
-    }
-    return Hij
+    return true;
   }
 
-  // get max off-diagonal value from Upper Diagonal
-  function getAij (Mij) {
-    const N = Mij.length
-    let maxMij = 0
-    let maxIJ = [0, 1]
-    for (let i = 0; i < N; i++) {
-      for (let j = i + 1; j < N; j++) {
-        if (Math.abs(maxMij) < Math.abs(Mij[i][j])) {
-          maxMij = Math.abs(Mij[i][j])
-          maxIJ = [i, j]
-        }
-      }
-    }
-    return [maxIJ, maxMij]
-  }
+  function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
+  {
+    if (prec === undefined)
+      prec = 1e-14;
 
-  // get max off-diagonal value from Upper Diagonal
-  function getAijBig (Mij) {
-    const N = Mij.length
-    let maxMij = 0
-    let maxIJ = [0, 1]
-    for (let i = 0; i < N; i++) {
-      for (let j = i + 1; j < N; j++) {
-        if (abs(maxMij) < abs(Mij[i][j])) {
-          maxMij = abs(Mij[i][j])
-          maxIJ = [i, j]
-        }
-      }
-    }
-    return [maxIJ, maxMij]
-  }
+    const size = mat.size();
 
-  // sort results
-  function sorting (E, S) {
-    const N = E.length
-    const Ef = Array(N)
-    const Sf = Array(N)
-    for (let k = 0; k < N; k++) {
-      Sf[k] = Array(N)
-    }
-    for (let i = 0; i < N; i++) {
-      let minID = 0
-      let minE = E[0]
-      for (let j = 0; j < E.length; j++) {
-        if (E[j] < minE) {
-          minID = j
-          minE = E[minID]
-        }
-      }
-      Ef[i] = E.splice(minID, 1)[0]
-      for (let k = 0; k < N; k++) {
-        Sf[k][i] = S[k][minID]
-        S[k].splice(minID, 1)
-      }
+    if (size.length !== 2 || size[0] !== size[1]) {
+      throw new RangeError('Matrix must be square (size: ' + format(size) + ')')
     }
-    return [clone(Ef), clone(Sf)]
-  }
 
-  /**
-   * Create an array of a certain size and fill all items with an initial value
-   * @param {number} size
-   * @param {number} value
-   * @return {number[]}
-   */
-  function createArray (size, value) {
-    // TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it)
-    const array = new Array(size)
+    const arr = mat.toArray();
+    const N = size[0];
 
-    for (let i = 0; i < size; i++) {
-      array[i] = value
-    }
+    if (isSymmetric(arr, N, prec))
+      return doRealSymetric(mat, N, prec)
 
-    return array
+    return doComplex(mat, prec);
   }
 
-  return eigs
-})
+  return eigs;
+})
\ No newline at end of file
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
new file mode 100644
index 0000000000..198d5ce30e
--- /dev/null
+++ b/src/function/matrix/eigs/complex.js
@@ -0,0 +1,14 @@
+export function createComplex()
+{
+    /**
+     * @param {Matrix} mat
+     * @param {number|BigNumber} prec
+     */
+    function main(mat, prec)
+    {
+        throw new Error('Not implemented yet.')
+    }
+
+
+    return main;
+}
\ No newline at end of file
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
new file mode 100644
index 0000000000..f7b0c763a3
--- /dev/null
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -0,0 +1,306 @@
+import { clone } from '../../../utils/object'
+
+
+export function createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
+
+
+
+  // check input for possible problems
+  // and perform diagonalization efficiently for
+  // specific type of number
+  function checkAndSubmit (x, n, precision = 1E-12) {
+    let type = x.datatype()
+    // vibe check
+    if (type === undefined) {
+      type = x.getDataType()
+    }
+    if (type !== 'number' && type !== 'BigNumber' && type !== 'Fraction') {
+      if (type === 'mixed') {
+        throw new TypeError('Mixed matrix element type is not supported')
+      } else {
+        throw new TypeError('Matrix element type not supported (' + type + ')')
+      }
+    }
+
+    // perform efficient calculation for 'numbers'
+    if (type === 'number') {
+      return diag(x.toArray())
+    } else if (type === 'Fraction') {
+      const xArr = x.toArray()
+      // convert fraction to numbers
+      for (let i = 0; i < n; i++) {
+        for (let j = i; j < n; j++) {
+          xArr[i][j] = xArr[i][j].valueOf()
+          xArr[j][i] = xArr[i][j]
+        }
+      }
+      return diag(x.toArray(), precision)
+    } else if (type === 'BigNumber') {
+      return diagBig(x.toArray(), precision)
+    }
+  }
+
+  // diagonalization implementation for number (efficient)
+  function diag (x, precision) {
+    const N = x.length
+    const e0 = Math.abs(precision / N)
+    let psi
+    let Sij = new Array(N)
+    // Sij is Identity Matrix
+    for (let i = 0; i < N; i++) {
+      Sij[i] = createArray(N, 0)
+      Sij[i][i] = 1.0
+    }
+    // initial error
+    let Vab = getAij(x)
+    while (Math.abs(Vab[1]) >= Math.abs(e0)) {
+      const i = Vab[0][0]
+      const j = Vab[0][1]
+      psi = getTheta(x[i][i], x[j][j], x[i][j])
+      x = x1(x, psi, i, j)
+      Sij = Sij1(Sij, psi, i, j)
+      Vab = getAij(x)
+    }
+    const Ei = createArray(N, 0) // eigenvalues
+    for (let i = 0; i < N; i++) {
+      Ei[i] = x[i][i]
+    }
+    return sorting(clone(Ei), clone(Sij))
+  }
+
+  // diagonalization implementation for bigNumber
+  function diagBig (x, precision) {
+    const N = x.length
+    const e0 = abs(precision / N)
+    let psi
+    let Sij = new Array(N)
+    // Sij is Identity Matrix
+    for (let i = 0; i < N; i++) {
+      Sij[i] = createArray(N, 0)
+      Sij[i][i] = 1.0
+    }
+    // initial error
+    let Vab = getAijBig(x)
+    while (abs(Vab[1]) >= abs(e0)) {
+      const i = Vab[0][0]
+      const j = Vab[0][1]
+      psi = getThetaBig(x[i][i], x[j][j], x[i][j])
+      x = x1Big(x, psi, i, j)
+      Sij = Sij1Big(Sij, psi, i, j)
+      Vab = getAijBig(x)
+    }
+    const Ei = createArray(N, 0) // eigenvalues
+    for (let i = 0; i < N; i++) {
+      Ei[i] = x[i][i]
+    }
+    // return [clone(Ei), clone(Sij)]
+    return sorting(clone(Ei), clone(Sij))
+  }
+
+  // get angle
+  function getTheta (aii, ajj, aij) {
+    let th = 0
+    const denom = (ajj - aii)
+    if (Math.abs(denom) <= 1E-14) {
+      th = Math.PI / 4.0
+    } else {
+      th = 0.5 * Math.atan(2.0 * aij / (ajj - aii))
+    }
+    return th
+  }
+
+  // get angle
+  function getThetaBig (aii, ajj, aij) {
+    let th = 0
+    const denom = subtract(ajj, aii)
+    if (abs(denom) <= 1E-14) {
+      th = Math.PI / 4.0
+    } else {
+      th = multiplyScalar(0.5, atan(multiply(2.0, aij, inv(denom))))
+    }
+    return th
+  }
+
+  // update eigvec
+  function Sij1 (Sij, theta, i, j) {
+    const N = Sij.length
+    const c = Math.cos(theta)
+    const s = Math.sin(theta)
+    const Ski = createArray(N, 0)
+    const Skj = createArray(N, 0)
+    for (let k = 0; k < N; k++) {
+      Ski[k] = c * Sij[k][i] - s * Sij[k][j]
+      Skj[k] = s * Sij[k][i] + c * Sij[k][j]
+    }
+    for (let k = 0; k < N; k++) {
+      Sij[k][i] = Ski[k]
+      Sij[k][j] = Skj[k]
+    }
+    return Sij
+  }
+  // update eigvec for overlap
+  function Sij1Big (Sij, theta, i, j) {
+    const N = Sij.length
+    const c = cos(theta)
+    const s = sin(theta)
+    const Ski = createArray(N, 0)
+    const Skj = createArray(N, 0)
+    for (let k = 0; k < N; k++) {
+      Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]))
+      Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]))
+    }
+    for (let k = 0; k < N; k++) {
+      Sij[k][i] = Ski[k]
+      Sij[k][j] = Skj[k]
+    }
+    return Sij
+  }
+
+  // update matrix
+  function x1Big (Hij, theta, i, j) {
+    const N = Hij.length
+    const c = bignumber(cos(theta))
+    const s = bignumber(sin(theta))
+    const c2 = multiplyScalar(c, c)
+    const s2 = multiplyScalar(s, s)
+    const Aki = createArray(N, 0)
+    const Akj = createArray(N, 0)
+    // 2cs Hij
+    const csHij = multiply(2, c, s, Hij[i][j])
+    //  Aii
+    const Aii = addScalar(subtract(multiplyScalar(c2, Hij[i][i]), csHij), multiplyScalar(s2, Hij[j][j]))
+    const Ajj = add(multiplyScalar(s2, Hij[i][i]), csHij, multiplyScalar(c2, Hij[j][j]))
+    // 0  to i
+    for (let k = 0; k < N; k++) {
+      Aki[k] = subtract(multiplyScalar(c, Hij[i][k]), multiplyScalar(s, Hij[j][k]))
+      Akj[k] = addScalar(multiplyScalar(s, Hij[i][k]), multiplyScalar(c, Hij[j][k]))
+    }
+    // Modify Hij
+    Hij[i][i] = Aii
+    Hij[j][j] = Ajj
+    Hij[i][j] = 0
+    Hij[j][i] = 0
+    // 0  to i
+    for (let k = 0; k < N; k++) {
+      if (k !== i && k !== j) {
+        Hij[i][k] = Aki[k]
+        Hij[k][i] = Aki[k]
+        Hij[j][k] = Akj[k]
+        Hij[k][j] = Akj[k]
+      }
+    }
+    return Hij
+  }
+
+  // update matrix
+  function x1 (Hij, theta, i, j) {
+    const N = Hij.length
+    const c = Math.cos(theta)
+    const s = Math.sin(theta)
+    const c2 = c * c
+    const s2 = s * s
+    const Aki = createArray(N, 0)
+    const Akj = createArray(N, 0)
+    //  Aii
+    const Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j]
+    const Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j]
+    // 0  to i
+    for (let k = 0; k < N; k++) {
+      Aki[k] = c * Hij[i][k] - s * Hij[j][k]
+      Akj[k] = s * Hij[i][k] + c * Hij[j][k]
+    }
+    // Modify Hij
+    Hij[i][i] = Aii
+    Hij[j][j] = Ajj
+    Hij[i][j] = 0
+    Hij[j][i] = 0
+    // 0  to i
+    for (let k = 0; k < N; k++) {
+      if (k !== i && k !== j) {
+        Hij[i][k] = Aki[k]
+        Hij[k][i] = Aki[k]
+        Hij[j][k] = Akj[k]
+        Hij[k][j] = Akj[k]
+      }
+    }
+    return Hij
+  }
+
+  // get max off-diagonal value from Upper Diagonal
+  function getAij (Mij) {
+    const N = Mij.length
+    let maxMij = 0
+    let maxIJ = [0, 1]
+    for (let i = 0; i < N; i++) {
+      for (let j = i + 1; j < N; j++) {
+        if (Math.abs(maxMij) < Math.abs(Mij[i][j])) {
+          maxMij = Math.abs(Mij[i][j])
+          maxIJ = [i, j]
+        }
+      }
+    }
+    return [maxIJ, maxMij]
+  }
+
+  // get max off-diagonal value from Upper Diagonal
+  function getAijBig (Mij) {
+    const N = Mij.length
+    let maxMij = 0
+    let maxIJ = [0, 1]
+    for (let i = 0; i < N; i++) {
+      for (let j = i + 1; j < N; j++) {
+        if (abs(maxMij) < abs(Mij[i][j])) {
+          maxMij = abs(Mij[i][j])
+          maxIJ = [i, j]
+        }
+      }
+    }
+    return [maxIJ, maxMij]
+  }
+
+  // sort results
+  function sorting (E, S) {
+    const N = E.length
+    const Ef = Array(N)
+    const Sf = Array(N)
+    for (let k = 0; k < N; k++) {
+      Sf[k] = Array(N)
+    }
+    for (let i = 0; i < N; i++) {
+      let minID = 0
+      let minE = E[0]
+      for (let j = 0; j < E.length; j++) {
+        if (E[j] < minE) {
+          minID = j
+          minE = E[minID]
+        }
+      }
+      Ef[i] = E.splice(minID, 1)[0]
+      for (let k = 0; k < N; k++) {
+        Sf[k][i] = S[k][minID]
+        S[k].splice(minID, 1)
+      }
+    }
+    return [clone(Ef), clone(Sf)]
+  }
+
+  /**
+   * Create an array of a certain size and fill all items with an initial value
+   * @param {number} size
+   * @param {number} value
+   * @return {number[]}
+   */
+  function createArray (size, value) {
+    // TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it)
+    const array = new Array(size)
+
+    for (let i = 0; i < size; i++) {
+      array[i] = value
+    }
+
+    return array
+  }
+
+  return checkAndSubmit;
+
+}

From f80749c761a33712809c6da2e305ca156eee6efe Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Fri, 14 Feb 2020 13:46:27 +0100
Subject: [PATCH 02/19] moved checks and coersions to eigs.js, made them more
 robust

---
 src/function/matrix/eigs.js              | 148 +++++++++++++++++++----
 src/function/matrix/eigs/complex.js      |   3 +-
 src/function/matrix/eigs/realSymetric.js |  48 +++-----
 src/utils/array.js                       |   2 +-
 4 files changed, 143 insertions(+), 58 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index a09d60bf40..9ba62401b3 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -1,5 +1,6 @@
 import { factory } from '../../utils/factory'
 import { format } from '../../utils/string'
+import { typeOf, isNumber, isBigNumber, isComplex, isFraction } from '../../utils/is'
 import { createRealSymmetric } from './eigs/realSymetric'
 import { createComplex } from './eigs/complex'
 
@@ -7,7 +8,9 @@ const name = 'eigs'
 const dependencies = [
   'typed',
   'matrix',
+  'number',
   'bignumber',
+  'complex',
   'add',
   'addScalar',
   'subtract',
@@ -16,13 +19,15 @@ const dependencies = [
   'abs',
   'equal',
   'larger',
+  're',
+  'im',
   'atan',
   'cos',
   'sin',
   'inv'
 ]
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, matrix, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, number, bignumber, multiply, add }) => {
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
    *
@@ -51,56 +56,153 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
 
     'Array': function (x) {
       const mat = matrix(x)
-      return computeValuesAndVectors(mat);
+      return computeValuesAndVectors(mat)
     },
 
     'Array, number|BigNumber': function (x, prec) {
       const mat = matrix(x)
-      return computeValuesAndVectors(mat, prec);
+      return computeValuesAndVectors(mat, prec)
     },
 
     'Matrix': function (mat) {
-      return computeValuesAndVectors(mat);
+      return computeValuesAndVectors(mat)
     },
 
     'Matrix, number|BigNumber': function (mat, prec) {
-      return computeValuesAndVectors(mat, prec);
+      return computeValuesAndVectors(mat, prec)
     }
   })
 
-  const doRealSymetric = createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add });
+  const doRealSymetric = createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add });
   const doComplex = createComplex();
 
-  /** @return {boolean} */
-  function isSymmetric(arr, N, prec) {
-
-    for (let i = 0; i < N; i++)
-      for (let j = i; j < N; j++)
-        if ( larger( subtract(arr[i][j], arr[j][i]), prec) )
-          return false;
-
-    return true;
-  }
 
   function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
   {
     if (prec === undefined)
       prec = 1e-14;
 
-    const size = mat.size();
+    const size = mat.size()
 
     if (size.length !== 2 || size[0] !== size[1]) {
       throw new RangeError('Matrix must be square (size: ' + format(size) + ')')
     }
 
-    const arr = mat.toArray();
-    const N = size[0];
+    const arr = mat.toArray()
+    const N = size[0]
+
+
+    if (isReal(arr, N, prec))
+    {
+      coerceReal(arr, N)
+
+      if (isSymmetric(arr, N, prec))
+      {
+        let type = coerceTypes(mat, arr, N)
+        return doRealSymetric(mat, N, prec, type)
+      }
+    }
+
+
+    let type = coerceTypes(mat, arr, N)
+    return doComplex(mat, prec, type)
+  }
+
+
+
+  /** @return {boolean} */
+  function isSymmetric(arr, N, prec) {
+
+    for (let i = 0; i < N; i++)
+    for (let j = i; j < N; j++)
+      if ( larger( subtract(arr[i][j], arr[j][i]), prec) )
+        return false
+
+    return true
+  }
+
+  /** @return {boolean} */
+  function isReal(arr, N, prec) {
+    for (let i = 0; i < N; i++)
+    for (let j = 0; j < N; j++)
+      if ( larger( im(arr[i][j]), prec) )
+        return false
+
+    return true
+  }
 
-    if (isSymmetric(arr, N, prec))
-      return doRealSymetric(mat, N, prec)
 
-    return doComplex(mat, prec);
+  function coerceReal(arr, N) {
+    for (let i = 0; i < N; i++)
+    for (let j = 0; j < N; j++)
+      arr[i][j] = re(arr[i][j])
+  }
+
+
+  /** @return {'number' | 'BigNumber' | 'Complex'} */
+  function coerceTypes(mat, arr, N) {
+
+    /** @type {string} */
+    let type = mat.datatype()
+
+    if (type == 'number' || type == 'BigNumber' || type == 'Complex') {
+      return type
+    }
+
+    let hasNumber  = false
+    let hasBig     = false
+    let hasComplex = false
+
+    for (let i = 0; i < N; i++)
+    for (let j = 0; j < N; j++)
+    {
+      el = arr[i][j]
+
+      if (isNumber(el) || isFraction(el))
+        hasNumber = true
+
+      else if (isBigNumber(el))
+        hasBig = true
+
+      else if (isComplex(el))
+        hasComplex = true
+
+      else
+        throw TypeError("Unsupported type in Matrix: " + typeOf(el))
+    }
+
+    if (hasBig && hasComplex)
+      console.warn("Complex BigNumbers not supported, this operation will lose precission.");
+
+    if (hasComplex)
+    {
+      for (let i = 0; i < N; i++)
+      for (let j = 0; j < N; j++)
+        arr[i][j] = complex(arr[i][j])
+
+      return 'Complex'
+    }
+
+    if (hasBig)
+    {
+      for (let i = 0; i < N; i++)
+      for (let j = 0; j < N; j++)
+        arr[i][j] = bignumber(arr[i][j])
+
+      return 'BigNumber'
+    }
+
+    if (hasNumber)
+    {
+      for (let i = 0; i < N; i++)
+      for (let j = 0; j < N; j++)
+        arr[i][j] = number(arr[i][j])
+
+      return 'number'
+    }
+
+    else throw TypeError('Matrix contains unsupported types only.')
   }
 
-  return eigs;
+  return eigs
 })
\ No newline at end of file
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 198d5ce30e..fdfbbeeea7 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -3,8 +3,9 @@ export function createComplex()
     /**
      * @param {Matrix} mat
      * @param {number|BigNumber} prec
+     * @param {'number'|'BigNumber'|'Complex'} type
      */
-    function main(mat, prec)
+    function main(mat, prec, type)
     {
         throw new Error('Not implemented yet.')
     }
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index f7b0c763a3..b9ceefdd69 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -5,39 +5,21 @@ export function createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos
 
 
 
-  // check input for possible problems
-  // and perform diagonalization efficiently for
-  // specific type of number
-  function checkAndSubmit (x, n, precision = 1E-12) {
-    let type = x.datatype()
-    // vibe check
-    if (type === undefined) {
-      type = x.getDataType()
-    }
-    if (type !== 'number' && type !== 'BigNumber' && type !== 'Fraction') {
-      if (type === 'mixed') {
-        throw new TypeError('Mixed matrix element type is not supported')
-      } else {
-        throw new TypeError('Matrix element type not supported (' + type + ')')
-      }
-    }
+  /**
+   * @param {number[] | BigNumber[]} arr
+   * @param {number} N
+   * @param {number} prec
+   * @param {'number' | 'BigNumber'} type
+   */
+  function main (arr, N, prec = 1E-12, type) {
 
-    // perform efficient calculation for 'numbers'
-    if (type === 'number') {
-      return diag(x.toArray())
-    } else if (type === 'Fraction') {
-      const xArr = x.toArray()
-      // convert fraction to numbers
-      for (let i = 0; i < n; i++) {
-        for (let j = i; j < n; j++) {
-          xArr[i][j] = xArr[i][j].valueOf()
-          xArr[j][i] = xArr[i][j]
-        }
-      }
-      return diag(x.toArray(), precision)
-    } else if (type === 'BigNumber') {
-      return diagBig(x.toArray(), precision)
-    }
+    if (type === 'number')
+      return diag(arr, prec)
+
+    if (type === 'BigNumber')
+      return diagBig(arr, prec)
+
+    throw TypeError('Unsupported data type: ' + type)
   }
 
   // diagonalization implementation for number (efficient)
@@ -301,6 +283,6 @@ export function createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos
     return array
   }
 
-  return checkAndSubmit;
+  return main;
 
 }
diff --git a/src/utils/array.js b/src/utils/array.js
index 2f3330a4e0..3b1ee48724 100644
--- a/src/utils/array.js
+++ b/src/utils/array.js
@@ -520,7 +520,7 @@ export function generalize (a) {
  * This method does not validate Array Matrix shape
  * @param {Array} array
  * @param {function} typeOf   Callback function to use to determine the type of a value
- * @return string
+ * @return {string}
  */
 export function getArrayDataType (array, typeOf) {
   let type // to hold type info

From ddc643e7e28acac0c8a69f7a8c25c0b749068071 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Fri, 14 Feb 2020 18:47:01 +0100
Subject: [PATCH 03/19] fix little bugs, make im and re more robust

---
 src/function/complex/im.js                   |  4 ++++
 src/function/complex/re.js                   |  4 ++++
 src/function/matrix/eigs.js                  | 21 ++++++++++++--------
 src/function/matrix/eigs/realSymetric.js     | 12 +++++++++--
 test/unit-tests/function/matrix/eigs.test.js |  2 +-
 5 files changed, 32 insertions(+), 11 deletions(-)

diff --git a/src/function/complex/im.js b/src/function/complex/im.js
index 04d3845a73..cec897266c 100644
--- a/src/function/complex/im.js
+++ b/src/function/complex/im.js
@@ -41,6 +41,10 @@ export const createIm = /* #__PURE__ */ factory(name, dependencies, ({ typed })
       return x.mul(0)
     },
 
+    Fraction: function (x) {
+      return x.mul(0)
+    },
+
     Complex: function (x) {
       return x.im
     },
diff --git a/src/function/complex/re.js b/src/function/complex/re.js
index 8b63fe7d75..6ad059d85e 100644
--- a/src/function/complex/re.js
+++ b/src/function/complex/re.js
@@ -41,6 +41,10 @@ export const createRe = /* #__PURE__ */ factory(name, dependencies, ({ typed })
       return x
     },
 
+    Fraction: function (x) {
+      return x
+    },
+
     Complex: function (x) {
       return x.re
     },
diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 9ba62401b3..529a55a301 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -8,6 +8,7 @@ const name = 'eigs'
 const dependencies = [
   'typed',
   'matrix',
+  'column',
   'number',
   'bignumber',
   'complex',
@@ -27,7 +28,7 @@ const dependencies = [
   'inv'
 ]
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, matrix, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, number, bignumber, multiply, add }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, matrix, column, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
    *
@@ -50,7 +51,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
    *
    * @param {Array | Matrix} x  Matrix to be diagonalized
    * @param {number | BigNumber} [prec] Precision, default value: 1e-15
-   * @return {{values: Array, vectors: Array} | {values: Matrix, vectors: Matrix}} Object containing eigenvalues (Array or Matrix) and eigenvectors (2D Array/Matrix).
+   * @return {{values: Array, vectors: Matrix[]}} Object containing an array of eigenvalues and an array of eigenvectors.
    */
   const eigs = typed('eigs', {
 
@@ -73,8 +74,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     }
   })
 
-  const doRealSymetric = createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add });
-  const doComplex = createComplex();
+  const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
+  const doComplex = createComplex()
 
 
   function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
@@ -82,6 +83,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     if (prec === undefined)
       prec = 1e-14;
 
+    prec = bignumber(prec)
+
     const size = mat.size()
 
     if (size.length !== 2 || size[0] !== size[1]) {
@@ -99,7 +102,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
       if (isSymmetric(arr, N, prec))
       {
         let type = coerceTypes(mat, arr, N)
-        return doRealSymetric(mat, N, prec, type)
+        return doRealSymetric(arr, N, prec, type)
       }
     }
 
@@ -115,7 +118,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
 
     for (let i = 0; i < N; i++)
     for (let j = i; j < N; j++)
-      if ( larger( subtract(arr[i][j], arr[j][i]), prec) )
+      // FIXME do proper comparison of bignum and frac
+      if ( larger( bignumber(abs(subtract(arr[i][j], arr[j][i]))), prec) )
         return false
 
     return true
@@ -125,7 +129,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   function isReal(arr, N, prec) {
     for (let i = 0; i < N; i++)
     for (let j = 0; j < N; j++)
-      if ( larger( im(arr[i][j]), prec) )
+      // FIXME do proper comparison of bignum and frac
+      if ( larger( bignumber(abs(im(arr[i][j]))), prec) )
         return false
 
     return true
@@ -156,7 +161,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     for (let i = 0; i < N; i++)
     for (let j = 0; j < N; j++)
     {
-      el = arr[i][j]
+      const el = arr[i][j]
 
       if (isNumber(el) || isFraction(el))
         hasNumber = true
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index b9ceefdd69..0834529679 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -1,7 +1,7 @@
 import { clone } from '../../../utils/object'
 
 
-export function createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
+export function createRealSymmetric({ addScalar, subtract, column, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
 
 
 
@@ -263,7 +263,15 @@ export function createRealSymmetric({ addScalar, subtract, equal, abs, atan, cos
         S[k].splice(minID, 1)
       }
     }
-    return [clone(Ef), clone(Sf)]
+
+    const values = clone(Ef)
+    const vectors = []
+
+    for (let i = 0; i < Sf.length; i++)
+      vectors[i] = column(Sf, i)
+
+    // FIXME vectors are always Array[], never Matrix[]
+    return { values, vectors }
   }
 
   /**
diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js
index 8a0dddc20e..a1dcf6c8c0 100644
--- a/test/unit-tests/function/matrix/eigs.test.js
+++ b/test/unit-tests/function/matrix/eigs.test.js
@@ -33,7 +33,7 @@ describe('eigs', function () {
       [[1, 0.1], [0.1, 1]]).values, [0.9, 1.1]
     )
     approx.deepEqual(eigs(
-      math.matrix([[1, 0.1], [0.1, 1]])).values, math.matrix([0.9, 1.1])
+      math.matrix([[1, 0.1], [0.1, 1]])).values, [0.9, 1.1]
     )
     approx.deepEqual(eigs(
       [[5, 2.3], [2.3, 1]]).values, [-0.04795013082563382, 6.047950130825635]

From 3d11fbda04078dd73087e2d3cc67a26285220d06 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Fri, 14 Feb 2020 21:40:19 +0100
Subject: [PATCH 04/19] Implemented matrix balancing algorithm

---
 src/function/matrix/eigs.js              |  11 ++-
 src/function/matrix/eigs/complex.js      | 114 ++++++++++++++++++++++-
 src/function/matrix/eigs/realSymetric.js |   2 +-
 3 files changed, 117 insertions(+), 10 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 529a55a301..bd928b9cb7 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -9,6 +9,7 @@ const dependencies = [
   'typed',
   'matrix',
   'column',
+  'diag',
   'number',
   'bignumber',
   'complex',
@@ -28,7 +29,7 @@ const dependencies = [
   'inv'
 ]
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, matrix, column, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, matrix, column, diag, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
    *
@@ -75,7 +76,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   })
 
   const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplex = createComplex()
+  const doComplex = createComplex({ add, multiply, abs, bignumber, diag })
 
 
   function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
@@ -108,7 +109,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
 
 
     let type = coerceTypes(mat, arr, N)
-    return doComplex(mat, prec, type)
+    return doComplex(arr, N, prec, type)
   }
 
 
@@ -118,7 +119,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
 
     for (let i = 0; i < N; i++)
     for (let j = i; j < N; j++)
-      // FIXME do proper comparison of bignum and frac
+      // TODO proper comparison of bignum and frac
       if ( larger( bignumber(abs(subtract(arr[i][j], arr[j][i]))), prec) )
         return false
 
@@ -129,7 +130,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   function isReal(arr, N, prec) {
     for (let i = 0; i < N; i++)
     for (let j = 0; j < N; j++)
-      // FIXME do proper comparison of bignum and frac
+      // TODO proper comparison of bignum and frac
       if ( larger( bignumber(abs(im(arr[i][j]))), prec) )
         return false
 
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index fdfbbeeea7..72a4060b79 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,15 +1,121 @@
-export function createComplex()
-{
+export function createComplex({add, multiply, abs, bignumber, diag}) {
     /**
-     * @param {Matrix} mat
+     * @param {number[][]} arr
+     * @param {number} N
      * @param {number|BigNumber} prec
      * @param {'number'|'BigNumber'|'Complex'} type
      */
-    function main(mat, prec, type)
+    function main(arr, N, prec, type)
     {
+        // TODO check if any row/col are zero except the diagonal
+
+        let R = balance(arr, N, prec, type)
+
+        console.log(arr, R)
+
         throw new Error('Not implemented yet.')
     }
 
+    function balance(arr, N, prec, type)
+    {
+        const big  = type === 'BigNumber'
+        const cplx = type === 'Complex'
+
+        const bigZero = bignumber(0)
+        const bigOne  = bignumber(1)
+
+        // base of the floating-point arithmetic
+        const radix = big ? bignumber(10) : 2
+        const radixSq = multiply(radix, radix)
+
+        // the diagonal transformation matrix R
+        const Rdiag = Array(+N).fill(big ? bigOne : 1)
+
+        // this isn't the only time we loop thru the matrix...
+        let last = false
+
+        while (!last)
+        {
+            // ...haha I'm joking! unless...
+            last = true
+
+            for (let i = 0; i < N; i++)
+            {
+
+                // compute the taxicab norm of i-th column and row
+                let colNorm = big ? bigZero : 0
+                let rowNorm = big ? bigZero : 0
+
+                for (let j = 0; j < N; j++)
+                {
+                    if (i === j) continue;
+                    let c = big||cplx ? arr[i][j].abs() : Math.abs(arr[i][j])
+                    let r = big||cplx ? arr[j][i].abs() : Math.abs(arr[j][i])
+                    colNorm = big ? colNorm.add(c) : colNorm + c
+                    rowNorm = big ? rowNorm.add(c) : rowNorm + r
+                }
+
+
+                if (big ? !colNorm.equals(0) && !rowNorm.equals(0) : colNorm !== 0 && rowNorm !== 0)
+                {
+                    // find integer power closest to balancing the matrix
+                    let f = big ? bignumber(1) : 1
+                    let c = colNorm
+
+                    const rowDivRadix = big ? rowNorm.div(radix) : rowNorm/radix
+                    const rowMulRadix = big ? rowNorm.mul(radix) : rowNorm*radix
+
+                    if (!big) {
+                        while (c < rowDivRadix) {
+                            c *= radixSq
+                            f *= radix
+                        }
+                        while (c > rowMulRadix) {
+                            c /= radixSq
+                            f /= radix
+                        }
+                    } else {
+                        while (c.lessThan(rowDivRadix)) {
+                            c = c.mul(radixSq)
+                            f = f.mul(radix)
+                        }
+                        while (c.greaterThan(rowMulRadix)) {
+                            c = c.div(radixSq)
+                            f = f.div(radix)
+                        }
+                    }
+
+
+                    // check whether balancing is needed
+                    const condition = !big
+                        ? (c + rowNorm)/f < 0.95*(colNorm + rowNorm)
+                        : c.add(rowNorm).div(f).lessThan(colNorm.add(rowNorm).mul(0.95))
+
+
+                    // apply balancing similarity transformation
+                    if (condition)
+                    {
+                        // okay bro we should definitely loop once again
+                        last = false
+
+                        let g = big ? bigOne.div(f) : 1/f
+                        R[i] = big ? R[i].mul(f) : R[i]*f
+
+                        for (let j = 0; j < N; j++)
+                        {
+                            arr[i][j] = multiply(arr[i][j], g)
+                            arr[j][i] = multiply(arr[j][i], f)
+                        }
+                    }
+
+                }
+            }
+        }
+
+        // return the diagonal transformation matrix
+        return diag(Rdiag);
+    }
+
 
     return main;
 }
\ No newline at end of file
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index 0834529679..dd160303ca 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -270,7 +270,7 @@ export function createRealSymmetric({ addScalar, subtract, column, abs, atan, co
     for (let i = 0; i < Sf.length; i++)
       vectors[i] = column(Sf, i)
 
-    // FIXME vectors are always Array[], never Matrix[]
+    // !FIXME vectors are always Array[], never Matrix[]
     return { values, vectors }
   }
 

From 97209abbfc0ed2cbbf8aafbd5c104edecd3e7a76 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Sat, 15 Feb 2020 11:19:00 +0100
Subject: [PATCH 05/19] fix typos

---
 src/function/matrix/eigs/complex.js | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 72a4060b79..5a994c64c3 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -99,12 +99,13 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
                         last = false
 
                         let g = big ? bigOne.div(f) : 1/f
-                        R[i] = big ? R[i].mul(f) : R[i]*f
+                        Rdiag[i] = big ? R[i].mul(f) : Rdiag[i]*f
 
                         for (let j = 0; j < N; j++)
                         {
-                            arr[i][j] = multiply(arr[i][j], g)
-                            arr[j][i] = multiply(arr[j][i], f)
+                            if (i === j) continue;
+                            arr[i][j] = multiply(arr[i][j], f)
+                            arr[j][i] = multiply(arr[j][i], g)
                         }
                     }
 

From 72f8e9d858583b8add89c2ce4b907e46034a2c82 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Sat, 15 Feb 2020 15:37:50 +0100
Subject: [PATCH 06/19] a draft of reduction to Hessenberg matrix

---
 src/function/matrix/eigs/complex.js | 72 +++++++++++++++++++++++++++++
 1 file changed, 72 insertions(+)

diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 5a994c64c3..a89aba5390 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -11,6 +11,12 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
 
         let R = balance(arr, N, prec, type)
 
+        // TODO if magnitudes of elements vary over many orders,
+        // move greatest elements to the top left corner
+
+        R = reduceToHessenberg(arg, N, prec, type).mul(R)
+
+
         console.log(arr, R)
 
         throw new Error('Not implemented yet.')
@@ -118,5 +124,71 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
     }
 
 
+    function reduceToHessenberg()
+    {
+        const big  = type === 'BigNumber'
+        const cplx = type === 'Complex'
+
+        const bigZero = bignumber(0)
+        const bigOne  = bignumber(1)
+
+        for (let i = 0; i < N-2; i++)
+        {
+            // Find the largest subdiag element in the i-th col
+
+            let maxIndex = 0
+            let max = big ? bigZero : 0
+
+            for (let j = i + 1; j < N; j++)
+            {
+                let el = abs(arr[j][i])
+                if (big ? max.abs().lessThan(el) : Math.abs(max) < el) {
+                    max = el
+                    maxIndex = i
+                }
+            }
+
+            // This col is pivoted, no need to do anything
+            if (big ? max.equals(bigZero) : max === 0)
+                continue
+
+
+            // Interchange maxIndex-th and (i+1)-th row
+            const tmp1 = arr[maxIndex]
+            arr[maxIndex] = arr[i+1]
+            arr[i+1] = tmp1
+
+            // Interchange maxIndex-th and (i+1)-th column
+            for (let j = 0; j < N; j++)
+            {
+                const tmp2 = arr[j][maxIndex]
+                arr[j][maxIndex] = arr[j][i+1]
+                arr[j][i+1] = tmp2
+            }
+
+            // TODO keep track of transformations
+
+            // Reduce following rows and columns
+            for (let j = i + 2; j < N; j++)
+            {
+                let n = !big ? arr[j][i] / div[i+1][i] : arr[j][i].div(arr[i+1][i])
+
+                // row
+                for (let k = 0; k < N; k++)
+                    arr[j][k] = !big ? arr[j][k] - n*arr[j][i+1] : arr[j][k].sub(n.mul(arr[j][i+1]))
+
+                // column
+                for (let k = 0; k < N; k++)
+                    arr[k][j] = !big ? arr[k][j] + n*arr[i+1][j] : arr[k][j].add(n.mul(arr[i+1][j]))
+
+                // TODO keep track of transformations
+            }
+        }
+
+        // !FIXME
+        return diag(Array(N).fill(1))
+    }
+
+
     return main;
 }
\ No newline at end of file

From aa4ee48bcf91f82556e65c2e9f9f688e9c4f170a Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Sun, 16 Feb 2020 17:44:11 +0100
Subject: [PATCH 07/19] finished implementation of reduction to Hessenberg

---
 src/function/matrix/eigs/complex.js | 73 +++++++++++++++++++----------
 1 file changed, 47 insertions(+), 26 deletions(-)

diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index a89aba5390..8542f2de84 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -14,14 +14,19 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
         // TODO if magnitudes of elements vary over many orders,
         // move greatest elements to the top left corner
 
-        R = reduceToHessenberg(arg, N, prec, type).mul(R)
-
-
+        reduceToHessenberg(arr, N, prec, type, R)
         console.log(arr, R)
 
         throw new Error('Not implemented yet.')
     }
 
+    /**
+     * @param {number[][]} arr
+     * @param {number} N
+     * @param {number} prec
+     * @param {'number'|'BigNumber'|'Complex'} type
+     * @returns {number[][]}
+     */
     function balance(arr, N, prec, type)
     {
         const big  = type === 'BigNumber'
@@ -49,6 +54,7 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
             {
 
                 // compute the taxicab norm of i-th column and row
+                // TODO optimize for complex numbers
                 let colNorm = big ? bigZero : 0
                 let rowNorm = big ? bigZero : 0
 
@@ -120,11 +126,18 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
         }
 
         // return the diagonal transformation matrix
-        return diag(Rdiag);
+        return diag(Rdiag)
     }
 
 
-    function reduceToHessenberg()
+    /**
+     * @param {number[][]} arr
+     * @param {number} N
+     * @param {number} prec
+     * @param {'number'|'BigNumber'|'Complex'} type
+     * @param {number[][]} R
+     */
+    function reduceToHessenberg(arr, N, prec, type, R)
     {
         const big  = type === 'BigNumber'
         const cplx = type === 'Complex'
@@ -144,7 +157,7 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
                 let el = abs(arr[j][i])
                 if (big ? max.abs().lessThan(el) : Math.abs(max) < el) {
                     max = el
-                    maxIndex = i
+                    maxIndex = j
                 }
             }
 
@@ -152,41 +165,49 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
             if (big ? max.equals(bigZero) : max === 0)
                 continue
 
+            if (maxIndex !== i+1)
+            {
+                // Interchange maxIndex-th and (i+1)-th row
+                const tmp1 = arr[maxIndex]
+                arr[maxIndex] = arr[i+1]
+                arr[i+1] = tmp1
 
-            // Interchange maxIndex-th and (i+1)-th row
-            const tmp1 = arr[maxIndex]
-            arr[maxIndex] = arr[i+1]
-            arr[i+1] = tmp1
+                // Interchange maxIndex-th and (i+1)-th column
+                for (let j = 0; j < N; j++)
+                {
+                    const tmp2 = arr[j][maxIndex]
+                    arr[j][maxIndex] = arr[j][i+1]
+                    arr[j][i+1] = tmp2
+                }
 
-            // Interchange maxIndex-th and (i+1)-th column
-            for (let j = 0; j < N; j++)
-            {
-                const tmp2 = arr[j][maxIndex]
-                arr[j][maxIndex] = arr[j][i+1]
-                arr[j][i+1] = tmp2
+                // keep track of transformations
+                const tmp3 = R[maxIndex]
+                R[maxIndex] = R[i+1]
+                R[i+1] = tmp3
             }
 
-            // TODO keep track of transformations
-
             // Reduce following rows and columns
             for (let j = i + 2; j < N; j++)
             {
-                let n = !big ? arr[j][i] / div[i+1][i] : arr[j][i].div(arr[i+1][i])
+                let n = !big ? arr[j][i] / max : arr[j][i].div(max)
 
-                // row
+                if (n == 0) continue;
+
+                // from j-th row subtract n-times (i+1)th row
                 for (let k = 0; k < N; k++)
-                    arr[j][k] = !big ? arr[j][k] - n*arr[j][i+1] : arr[j][k].sub(n.mul(arr[j][i+1]))
+                    arr[j][k] = !big ? arr[j][k] - n*arr[i+1][k] : arr[j][k].sub(n.mul(arr[i+1][k]))
 
-                // column
+                // to (i+1)th column add n-times j-th column
                 for (let k = 0; k < N; k++)
-                    arr[k][j] = !big ? arr[k][j] + n*arr[i+1][j] : arr[k][j].add(n.mul(arr[i+1][j]))
+                    arr[k][i+1] = !big ? arr[k][i+1] + n*arr[k][j] : arr[k][i+1].add(n.mul(arr[k][j]))
 
-                // TODO keep track of transformations
+                // keep track of transformations
+                for (let k = 0; k < N; k++)
+                    R[j][k] = !big ? R[j][k] - n*R[i+1][k] : R[j][k].sub(n.mul(R[i+1][k]))
             }
         }
 
-        // !FIXME
-        return diag(Array(N).fill(1))
+        return R
     }
 
 

From a32d6e8efdabed6cfe542fa23281867cd3150970 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Thu, 20 Feb 2020 16:51:53 +0100
Subject: [PATCH 08/19] fix Hessenberg elimination for complex numbers

---
 src/function/matrix/eigs.js         |  2 +-
 src/function/matrix/eigs/complex.js | 18 +++++++++---------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index bd928b9cb7..532b8ea737 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -76,7 +76,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   })
 
   const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplex = createComplex({ add, multiply, abs, bignumber, diag })
+  const doComplex = createComplex({ add, subtract, multiply, abs, bignumber, diag })
 
 
   function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 8542f2de84..7d5b2c6bb1 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,4 +1,4 @@
-export function createComplex({add, multiply, abs, bignumber, diag}) {
+export function createComplex({add, subtract, multiply, abs, bignumber, diag}) {
     /**
      * @param {number[][]} arr
      * @param {number} N
@@ -141,9 +141,9 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
     {
         const big  = type === 'BigNumber'
         const cplx = type === 'Complex'
+        const num = !big && !cplx
 
         const bigZero = bignumber(0)
-        const bigOne  = bignumber(1)
 
         for (let i = 0; i < N-2; i++)
         {
@@ -154,15 +154,15 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
 
             for (let j = i + 1; j < N; j++)
             {
-                let el = abs(arr[j][i])
-                if (big ? max.abs().lessThan(el) : Math.abs(max) < el) {
+                let el = arr[j][i]
+                if (big ? abs(max).lessThan(abs(el)) : abs(max) < abs(el)) {
                     max = el
                     maxIndex = j
                 }
             }
 
             // This col is pivoted, no need to do anything
-            if (big ? max.equals(bigZero) : max === 0)
+            if (big ? abs(max).lessThan(prec) : abs(max) < prec)
                 continue
 
             if (maxIndex !== i+1)
@@ -189,21 +189,21 @@ export function createComplex({add, multiply, abs, bignumber, diag}) {
             // Reduce following rows and columns
             for (let j = i + 2; j < N; j++)
             {
-                let n = !big ? arr[j][i] / max : arr[j][i].div(max)
+                let n = num ? arr[j][i] / max : arr[j][i].div(max)
 
                 if (n == 0) continue;
 
                 // from j-th row subtract n-times (i+1)th row
                 for (let k = 0; k < N; k++)
-                    arr[j][k] = !big ? arr[j][k] - n*arr[i+1][k] : arr[j][k].sub(n.mul(arr[i+1][k]))
+                    arr[j][k] = num ? arr[j][k] - n*arr[i+1][k] : arr[j][k].sub(n.mul(arr[i+1][k]))
 
                 // to (i+1)th column add n-times j-th column
                 for (let k = 0; k < N; k++)
-                    arr[k][i+1] = !big ? arr[k][i+1] + n*arr[k][j] : arr[k][i+1].add(n.mul(arr[k][j]))
+                    arr[k][i+1] = num ? arr[k][i+1] + n*arr[k][j] : arr[k][i+1].add(n.mul(arr[k][j]))
 
                 // keep track of transformations
                 for (let k = 0; k < N; k++)
-                    R[j][k] = !big ? R[j][k] - n*R[i+1][k] : R[j][k].sub(n.mul(R[i+1][k]))
+                    R[j][k] = num ? R[j][k] - n*R[i+1][k] : subtract(R[j][k], n.mul(R[i+1][k]))
             }
         }
 

From 1e4a667315a37d2243660ea0513a9777660a89bd Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Thu, 20 Feb 2020 22:54:46 +0100
Subject: [PATCH 09/19] implemented non-shifted explicit QR algorithm for real
 matrices

---
 src/function/matrix/eigs.js         |   8 ++-
 src/function/matrix/eigs/complex.js | 107 +++++++++++++++++++++++++++-
 2 files changed, 109 insertions(+), 6 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 532b8ea737..78247d90d9 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -23,13 +23,15 @@ const dependencies = [
   'larger',
   're',
   'im',
+  'sqrt',
   'atan',
   'cos',
   'sin',
-  'inv'
+  'inv',
+  'qr'
 ]
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, matrix, column, diag, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, qr, sqrt, matrix, column, diag, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
    *
@@ -76,7 +78,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   })
 
   const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplex = createComplex({ add, subtract, multiply, abs, bignumber, diag })
+  const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr })
 
 
   function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 7d5b2c6bb1..3081c92c46 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,4 +1,4 @@
-export function createComplex({add, subtract, multiply, abs, bignumber, diag}) {
+export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr}) {
     /**
      * @param {number[][]} arr
      * @param {number} N
@@ -7,6 +7,11 @@ export function createComplex({add, subtract, multiply, abs, bignumber, diag}) {
      */
     function main(arr, N, prec, type)
     {
+        // TODO implement QR for complex matrices
+        if (type === 'Complex')
+            throw new TypeError('Complex matrices not yet supported')
+
+
         // TODO check if any row/col are zero except the diagonal
 
         let R = balance(arr, N, prec, type)
@@ -15,9 +20,12 @@ export function createComplex({add, subtract, multiply, abs, bignumber, diag}) {
         // move greatest elements to the top left corner
 
         reduceToHessenberg(arr, N, prec, type, R)
-        console.log(arr, R)
 
-        throw new Error('Not implemented yet.')
+        let values = iterateUntilDiagonal(arr, N, prec, type, R)
+
+        // TODO find vectors
+
+        return { values }
     }
 
     /**
@@ -210,6 +218,99 @@ export function createComplex({add, subtract, multiply, abs, bignumber, diag}) {
         return R
     }
 
+    function iterateUntilDiagonal(arr, N, prec, type, R)
+    {
+        const big = type === 'BigNumber'
+        const lambdas = []
+        let lastConvergenceBefore = 0;
+
+        while (lastConvergenceBefore <= 30)
+        {
+            lastConvergenceBefore += 1
+
+            // check the dimensions of matrix before factoring
+            if (N == 1) {
+                lambdas.push(arr[0][0])
+                break;
+            }
+            if (N == 2) {
+                let ll = eigenvalues2x2(
+                    arr[0][0], arr[0][1],
+                    arr[1][0], arr[1][1]
+                )
+                lambdas.push(...ll)
+                break;
+            }
+
+
+            // perform the factorization
+
+            let k = 0; // TODO set close to an eigenvalue
+
+            for (let i = 0; i < N; i++)
+                arr[i][i] -= k
+
+            // TODO do an implicit QR transformation
+            let {Q, R} = qr(arr)
+            arr = multiply(R, Q)
+
+            for (let i = 0; i < N; i++)
+                arr[i][i] += k
+
+
+            // The last element converged to an eigenvalue
+            if (big ? abs(arr[N-1][N-2]).lessThan(prec) : abs(arr[N-1][N-2]) < prec)
+            {
+                lastConvergenceBefore = 0
+                lambdas.push(arr[N-1][N-1])
+
+                // reduce the matrix size
+                N -= 1
+                arr.pop()
+                for (let i = 0; i < N; i++)
+                    arr[i].pop()
+            }
+
+            // The last 2x2 block matrix converged
+            else if (big ? abs(arr[N-2][N-3]).lessThan(prec) : abs(arr[N-2][N-3]) < prec)
+            {
+                lastConvergenceBefore = 0
+                let ll = eigenvalues2x2(
+                    arr[N-2][N-2], arr[N-2][N-1],
+                    arr[N-1][N-2], arr[N-1][N-1]
+                )
+                lambdas.push(...ll)
+
+                // reduce the matrix size
+                N -= 2
+                arr.pop()
+                arr.pop()
+                for (let i = 0; i < N; i++){
+                    arr[i].pop()
+                    arr[i].pop()
+                }
+            }
+        }
+
+        return lambdas
+    }
+
+
+    /**
+     * Compute the eigenvalues of an 2x2 matrix
+     * @return {[number,number]}
+     */
+    function eigenvalues2x2(a,b,c,d)
+    {
+        // λ± = ½ trA ± √( tr²A - 4 detA )
+        const trA = addScalar(a, d)
+        const detA = subtract(multiply(a, d), multiply(b, c))
+        const x = multiply(trA, 0.5)
+        const y = multiply(sqrt( subtract(multiply(trA, trA), multiply(4, detA)) ), 0.5)
+
+        return [addScalar(x, y), subtract(x, y)]
+    }
+
 
     return main;
 }
\ No newline at end of file

From 58ad520fb8754dc074684f22578ebffd9f390b87 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Sat, 22 Feb 2020 20:59:21 +0100
Subject: [PATCH 10/19] implemented vector computation, won't work untill
 usolve is fixed

---
 src/function/matrix/eigs.js         |   7 +-
 src/function/matrix/eigs/complex.js | 328 +++++++++++++++++++++++-----
 2 files changed, 277 insertions(+), 58 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 78247d90d9..b5a65ac285 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -28,10 +28,11 @@ const dependencies = [
   'cos',
   'sin',
   'inv',
-  'qr'
+  'qr',
+  'usolve'
 ]
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, qr, sqrt, matrix, column, diag, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, qr, usolve, sqrt, matrix, column, diag, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
    *
@@ -78,7 +79,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   })
 
   const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr })
+  const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr, inv, usolve })
 
 
   function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 3081c92c46..30f6e72119 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,12 +1,17 @@
-export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr}) {
+import { clone } from '../../../utils/object'
+
+export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumber, diag, inv, qr, usolve}) {
     /**
-     * @param {number[][]} arr
-     * @param {number} N
-     * @param {number|BigNumber} prec
+     * @param {number[][]} arr the matrix to find eigenvalues of
+     * @param {number} N size of the matrix
+     * @param {number|BigNumber} prec precision, anything lower will be considered zero
      * @param {'number'|'BigNumber'|'Complex'} type
+     * @param {boolean} findVectors should we find eigenvectors?
      */
-    function main(arr, N, prec, type)
+    function main(arr, N, prec, type, findVectors)
     {
+        if (findVectors === undefined) findVectors = true
+
         // TODO implement QR for complex matrices
         if (type === 'Complex')
             throw new TypeError('Complex matrices not yet supported')
@@ -14,18 +19,39 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
 
         // TODO check if any row/col are zero except the diagonal
 
-        let R = balance(arr, N, prec, type)
+        // make sure corresponding rows and columns have similar magnitude
+        // important because of numerical stability
+        let R = balance(arr, N, prec, type, findVectors)
+
+        // R is the row transformation matrix
+        // A' = R A R⁻¹, A is the original matrix
+        // (if findVectors is false, R is undefined)
 
         // TODO if magnitudes of elements vary over many orders,
         // move greatest elements to the top left corner
 
-        reduceToHessenberg(arr, N, prec, type, R)
+        // using similarity transformations, reduce the matrix
+        // to Hessenberg form (upper triangular plus one subdiagonal row)
+        // updates the transformation matrix R with new row operationsq
+        reduceToHessenberg(arr, N, prec, type, findVectors, R)
+
+        // find eigenvalues
+        let { values, C } = iterateUntilTriangular(arr, N, prec, type, findVectors)
 
-        let values = iterateUntilDiagonal(arr, N, prec, type, R)
+        // values is the list of eigenvalues, C is the column
+        // transformation matrix that transforms the hessenberg
+        // matrix to upper triangular
 
-        // TODO find vectors
+        // compose transformations A → hess. and hess. → triang.
+        C = multiply( inv(R), C )
 
-        return { values }
+
+        let vectors
+
+        if (findVectors)
+            vectors = findEigenvectors(arr, N, C, values)
+
+        return { values, vectors }
     }
 
     /**
@@ -35,7 +61,7 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
      * @param {'number'|'BigNumber'|'Complex'} type
      * @returns {number[][]}
      */
-    function balance(arr, N, prec, type)
+    function balance(arr, N, prec, type, findVectors)
     {
         const big  = type === 'BigNumber'
         const cplx = type === 'Complex'
@@ -48,7 +74,9 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
         const radixSq = multiply(radix, radix)
 
         // the diagonal transformation matrix R
-        const Rdiag = Array(+N).fill(big ? bigOne : 1)
+        let Rdiag
+        if (findVectors)
+            Rdiag = Array(N).fill(big ? bigOne : 1)
 
         // this isn't the only time we loop thru the matrix...
         let last = false
@@ -76,10 +104,14 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
                 }
 
 
-                if (big ? !colNorm.equals(0) && !rowNorm.equals(0) : colNorm !== 0 && rowNorm !== 0)
+                if (!big ? (colNorm !== 0 && rowNorm !== 0) : (!colNorm.equals(0) && !rowNorm.equals(0)) )
                 {
                     // find integer power closest to balancing the matrix
-                    let f = big ? bignumber(1) : 1
+                    // (we want to scale only by integer powers of radix,
+                    // so that we don't lose any precision due to round-off)
+
+
+                    let f = big ? bigOne : 1
                     let c = colNorm
 
                     const rowDivRadix = big ? rowNorm.div(radix) : rowNorm/radix
@@ -115,11 +147,11 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
                     // apply balancing similarity transformation
                     if (condition)
                     {
-                        // okay bro we should definitely loop once again
+                        // we should loop once again to check whether
+                        // another rebalancing is needed
                         last = false
 
                         let g = big ? bigOne.div(f) : 1/f
-                        Rdiag[i] = big ? R[i].mul(f) : Rdiag[i]*f
 
                         for (let j = 0; j < N; j++)
                         {
@@ -127,13 +159,17 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
                             arr[i][j] = multiply(arr[i][j], f)
                             arr[j][i] = multiply(arr[j][i], g)
                         }
+
+                        // keep track of transformations
+                        if (findVectors)
+                            Rdiag[i] = big ? R[i].mul(f) : Rdiag[i]*f
                     }
 
                 }
             }
         }
 
-        // return the diagonal transformation matrix
+        // return the diagonal row transformation matrix
         return diag(Rdiag)
     }
 
@@ -143,9 +179,10 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
      * @param {number} N
      * @param {number} prec
      * @param {'number'|'BigNumber'|'Complex'} type
-     * @param {number[][]} R
+     * @param {boolean} findVectors
+     * @param {number[][]} R the row transformation matrix that will be modified
      */
-    function reduceToHessenberg(arr, N, prec, type, R)
+    function reduceToHessenberg(arr, N, prec, type, findVectors, R)
     {
         const big  = type === 'BigNumber'
         const cplx = type === 'Complex'
@@ -189,9 +226,11 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
                 }
 
                 // keep track of transformations
-                const tmp3 = R[maxIndex]
-                R[maxIndex] = R[i+1]
-                R[i+1] = tmp3
+                if (findVectors) {
+                    const tmp3 = R[maxIndex]
+                    R[maxIndex] = R[i+1]
+                    R[i+1] = tmp3
+                }
             }
 
             // Reduce following rows and columns
@@ -210,89 +249,178 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
                     arr[k][i+1] = num ? arr[k][i+1] + n*arr[k][j] : arr[k][i+1].add(n.mul(arr[k][j]))
 
                 // keep track of transformations
-                for (let k = 0; k < N; k++)
-                    R[j][k] = num ? R[j][k] - n*R[i+1][k] : subtract(R[j][k], n.mul(R[i+1][k]))
+                if (findVectors)
+                    for (let k = 0; k < N; k++)
+                        R[j][k] = num ? R[j][k] - n*R[i+1][k] : subtract(R[j][k], n.mul(R[i+1][k]))
             }
         }
 
         return R
     }
 
-    function iterateUntilDiagonal(arr, N, prec, type, R)
+    /**
+     * @returns {{values: values, C: Matrix}}
+     * @see Press, Wiliams: Numerical recipes in Fortran 77
+     * @see https://en.wikipedia.org/wiki/QR_algorithm
+     */
+    function iterateUntilTriangular(A, N, prec, type, findVectors)
     {
         const big = type === 'BigNumber'
+
+        // The Francis Algorithm
+        // The core idea of this algorithm is that doing successive
+        // A' = Q⁺AQ transformations will eventually converge to block-
+        // upper-triangular with diagonal blocks either 1x1 or 2x2.
+        // The Q here is the one from the QR decomposition, A = QR.
+        // Since the eigenvalues of a block-upper-triangular matrix are
+        // the eigenvalues of its diagonal blocks and we know how to find
+        // eigenvalues of a 2x2 matrix, we know the eigenvalues of A.
+
+
+        let arr = clone(A)
+
+        // the list of converged eigenvalues
         const lambdas = []
+
+        // size of arr, which will get smaller as eigenvalues converge
+        let n = N;
+
+        // the diagonal of the block-diagonal matrix that turns
+        // converged 2x2 matrices into upper triangular matrices
+        const Sdiag = []
+
+        // N×N matrix describing the overall transformation done during the QR algorithm
+        let Qtotal = findVectors ? diag(Array(N).fill(1)) : undefined
+
+        // n×n matrix describing the QR transformations done since last convergence
+        let Qpartial = findVectors ? diag(Array(n).fill(1)) : undefined
+
+        // last eigenvalue converged before this many steps
         let lastConvergenceBefore = 0;
 
+
         while (lastConvergenceBefore <= 30)
         {
             lastConvergenceBefore += 1
 
-            // check the dimensions of matrix before factoring
-            if (N == 1) {
-                lambdas.push(arr[0][0])
-                break;
-            }
-            if (N == 2) {
-                let ll = eigenvalues2x2(
-                    arr[0][0], arr[0][1],
-                    arr[1][0], arr[1][1]
-                )
-                lambdas.push(...ll)
-                break;
-            }
-
+            // TODO if the convergence is slow, do something clever
 
-            // perform the factorization
+            // Perform the factorization
 
             let k = 0; // TODO set close to an eigenvalue
 
-            for (let i = 0; i < N; i++)
+            for (let i = 0; i < n; i++)
                 arr[i][i] -= k
 
             // TODO do an implicit QR transformation
             let {Q, R} = qr(arr)
             arr = multiply(R, Q)
 
-            for (let i = 0; i < N; i++)
+            for (let i = 0; i < n; i++)
                 arr[i][i] += k
 
 
-            // The last element converged to an eigenvalue
-            if (big ? abs(arr[N-1][N-2]).lessThan(prec) : abs(arr[N-1][N-2]) < prec)
+            // keep track of transformations
+            if (findVectors)
+                Qpartial = multiply(Qpartial, Q)
+
+
+            // The rightmost diagonal element converged to an eigenvalue
+            if (n == 1 || (big ? abs(arr[n-1][n-2]).lessThan(prec) : abs(arr[n-1][n-2]) < prec))
             {
                 lastConvergenceBefore = 0
-                lambdas.push(arr[N-1][N-1])
+                lambdas.push(arr[n-1][n-1])
+
+                // keep track of transformations
+                if (findVectors) {
+                    Sdiag.unshift( [[1]] )
+                    inflateMatrix(Qpartial, N)
+                    Qtotal = multiply(Qtotal, Qpartial)
+
+                    if (n > 1)
+                        Qpartial = diag( Array(n-1).fill(1) )
+                }
 
                 // reduce the matrix size
-                N -= 1
+                n -= 1
                 arr.pop()
-                for (let i = 0; i < N; i++)
+                for (let i = 0; i < n; i++)
                     arr[i].pop()
             }
 
-            // The last 2x2 block matrix converged
-            else if (big ? abs(arr[N-2][N-3]).lessThan(prec) : abs(arr[N-2][N-3]) < prec)
+            // The rightmost diagonal 2x2 block converged
+            else if (n == 2 || (big ? abs(arr[n-2][n-3]).lessThan(prec) : abs(arr[n-2][n-3]) < prec))
             {
                 lastConvergenceBefore = 0
                 let ll = eigenvalues2x2(
-                    arr[N-2][N-2], arr[N-2][N-1],
-                    arr[N-1][N-2], arr[N-1][N-1]
+                    arr[n-2][n-2], arr[n-2][n-1],
+                    arr[n-1][n-2], arr[n-1][n-1]
                 )
                 lambdas.push(...ll)
 
+                // keep track of transformations
+                if (findVectors) {
+                    Sdiag.unshift(jordanBase2x2(
+                        arr[n-2][n-2], arr[n-2][n-1],
+                        arr[n-1][n-2], arr[n-1][n-1],
+                        ll[0], ll[1], prec, type
+                    ))
+                    inflateMatrix(Qpartial, N)
+                    Qtotal = multiply(Qtotal, Qpartial)
+
+                    if (n > 2)
+                        Qpartial = diag( Array(n-2).fill(1) )
+                }
+
                 // reduce the matrix size
-                N -= 2
+                n -= 2
                 arr.pop()
                 arr.pop()
-                for (let i = 0; i < N; i++){
+                for (let i = 0; i < n; i++) {
                     arr[i].pop()
                     arr[i].pop()
                 }
             }
+
+            if (n == 0) break;
+        }
+
+        // the algorithm didn't converge
+        if (lastConvergenceBefore > 30)
+            throw Error('The eigenvalues failed to converge. Only found these eigenvalues: '+values.join(', '))
+
+        // combine the overall QR transformation Qtotal with the subsequent
+        // transformation S that turns the diagonal 2x2 blocks to upper triangular
+        let C = findVectors ? multiply( Qtotal, blockDiag(Sdiag, N) ) : undefined
+
+        return { values: lambdas, C }
+    }
+
+
+    /**
+     * @param {Matrix} A original matrix
+     * @param {number} N size of A
+     * @param {Matrix} C column transformation matrix that turns A into upper triangular
+     * @param {number[]} values array of eigenvalues of A
+     * @returns {Matrix[]} eigenvalues
+     */
+    function findEigenvectors(A, N, C, values)
+    {
+        const Cinv = inv(C)
+        const U = multiply( Cinv, A, C )
+
+        const b = Array(N).fill(0)
+        const E = diag(Array(N).fill(1))
+
+        let vectors = []
+
+        for (const l of values) {
+            // !FIXME this is a mock implementation
+            const v = usolve( subtract(U, multiply(l, E)), b )
+            vectors.push( multiply(C, v) )
         }
 
-        return lambdas
+        return vectors
     }
 
 
@@ -302,7 +430,7 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
      */
     function eigenvalues2x2(a,b,c,d)
     {
-        // λ± = ½ trA ± √( tr²A - 4 detA )
+        // λ± = ½ trA ± ½ √( tr²A - 4 detA )
         const trA = addScalar(a, d)
         const detA = subtract(multiply(a, d), multiply(b, c))
         const x = multiply(trA, 0.5)
@@ -312,5 +440,95 @@ export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumb
     }
 
 
+    /**
+     * For an 2x2 matrix compute the transformation matrix S,
+     * so that SAS⁻¹ is an upper triangular matrix
+     * @return {[[number,number],[number,number]]}
+     * @see https://math.berkeley.edu/~ogus/old/Math_54-05/webfoils/jordan.pdf
+     * @see http://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html
+     */
+    function jordanBase2x2(a,b,c,d,l1,l2,prec,type)
+    {
+        const big = type === 'BigNumber'
+
+        const b0 = big ? 0 : bignumber(0)
+        const b1 = big ? 1 : bignumber(1)
+
+
+        // matrix is already upper triangular
+        // return an identity matrix
+        if (big ? abs(c).lessThan(prec) : abs(c) < prec)
+            return [ [b1, b0], [b0, b1] ]
+
+        // matrix is diagonalizable
+        // return its eigenvectors as columns
+        if (big ? abs(l1.sub(l2)).greaterThan(prec) : abs(subtract(l1, l2)) > prec)
+            return [ [ subtract(l1, d), subtract(l2, d) ], [ c, c ] ]
+
+
+        // matrix is not diagonalizable
+        // compute off-diagonal elements of N = A - λI
+        // N₁₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ )
+        // N₁₂ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ )
+
+        let na = subtract(a, l1)
+        let nb = subtract(b, l1)
+        let nc = subtract(c, l1)
+        let nd = subtract(d, l1)
+
+        if (big ? abs(nb).lessThan(prec) : abs(nb) < prec)
+            return [ [na, b1], [nc, b0] ]
+        else
+            return [ [nb, b0], [nd, b1] ]
+    }
+
+
+    /**
+     * Enlarge the matrix from n×n to N×N, setting the new
+     * elements to 1 on diagonal and 0 elsewhere
+     */
+    function inflateMatrix(arr, N)
+    {
+        // add columns
+        for (let i = 0; i < arr.length; i++)
+            arr[i].push( ...Array(N-arr[i].length).fill(0) )
+
+        // add rows
+        for (let i = arr.length; i < N; i++) {
+            arr.push( Array(N).fill(0) )
+            arr[i][i] = 1
+        }
+
+        return arr
+    }
+
+
+    /**
+     * Create a block-diagonal matrix with the given square matrices on the diagonal
+     * @param {Matrix[] | number[][][]} arr array of matrices to be placed on the diagonal
+     * @param {number} N the size of the resulting matrix
+     */
+    function blockDiag(arr, N)
+    {
+        const M = []
+        for (let i = 0; i < N; i++)
+            M[i] = Array(N).fill(0)
+
+        let I = 0;
+        for (const sub of arr)
+        {
+            const n = sub.length
+
+            for (let i = 0; i < n; i++)
+            for (let j = 0; j < n; j++)
+                M[I+i][I+j] = sub[i][j]
+
+            I += n
+        }
+
+        return M
+    }
+
+
     return main;
 }
\ No newline at end of file

From e7bf6daf2d20000cf32ca131333f1782ccc2bf33 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Sat, 22 Feb 2020 21:38:29 +0100
Subject: [PATCH 11/19] refactored to match yarn lint

---
 src/function/matrix/eigs.js              | 155 ++--
 src/function/matrix/eigs/complex.js      | 889 +++++++++++------------
 src/function/matrix/eigs/realSymetric.js |  19 +-
 3 files changed, 522 insertions(+), 541 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index b5a65ac285..d22faad013 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -59,7 +59,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
    */
   const eigs = typed('eigs', {
 
-    'Array': function (x) {
+    Array: function (x) {
       const mat = matrix(x)
       return computeValuesAndVectors(mat)
     },
@@ -69,7 +69,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
       return computeValuesAndVectors(mat, prec)
     },
 
-    'Matrix': function (mat) {
+    Matrix: function (mat) {
       return computeValuesAndVectors(mat)
     },
 
@@ -81,11 +81,10 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
   const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
   const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr, inv, usolve })
 
-
-  function computeValuesAndVectors(/**@type {Matrix}*/ mat, prec)
-  {
-    if (prec === undefined)
-      prec = 1e-14;
+  function computeValuesAndVectors (mat, prec) {
+    if (prec === undefined) {
+      prec = 1e-14
+    }
 
     prec = bignumber(prec)
 
@@ -98,120 +97,120 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     const arr = mat.toArray()
     const N = size[0]
 
-
-    if (isReal(arr, N, prec))
-    {
+    if (isReal(arr, N, prec)) {
       coerceReal(arr, N)
 
-      if (isSymmetric(arr, N, prec))
-      {
-        let type = coerceTypes(mat, arr, N)
+      if (isSymmetric(arr, N, prec)) {
+        const type = coerceTypes(mat, arr, N)
         return doRealSymetric(arr, N, prec, type)
       }
     }
 
-
-    let type = coerceTypes(mat, arr, N)
+    const type = coerceTypes(mat, arr, N)
     return doComplex(arr, N, prec, type)
   }
 
-
-
   /** @return {boolean} */
-  function isSymmetric(arr, N, prec) {
-
-    for (let i = 0; i < N; i++)
-    for (let j = i; j < N; j++)
-      // TODO proper comparison of bignum and frac
-      if ( larger( bignumber(abs(subtract(arr[i][j], arr[j][i]))), prec) )
-        return false
+  function isSymmetric (arr, N, prec) {
+    for (let i = 0; i < N; i++) {
+      for (let j = i; j < N; j++) {
+        // TODO proper comparison of bignum and frac
+        if (larger(bignumber(abs(subtract(arr[i][j], arr[j][i]))), prec)) {
+          return false
+        }
+      }
+    }
 
     return true
   }
 
   /** @return {boolean} */
-  function isReal(arr, N, prec) {
-    for (let i = 0; i < N; i++)
-    for (let j = 0; j < N; j++)
-      // TODO proper comparison of bignum and frac
-      if ( larger( bignumber(abs(im(arr[i][j]))), prec) )
-        return false
+  function isReal (arr, N, prec) {
+    for (let i = 0; i < N; i++) {
+      for (let j = 0; j < N; j++) {
+        // TODO proper comparison of bignum and frac
+        if (larger(bignumber(abs(im(arr[i][j]))), prec)) {
+          return false
+        }
+      }
+    }
 
     return true
   }
 
-
-  function coerceReal(arr, N) {
-    for (let i = 0; i < N; i++)
-    for (let j = 0; j < N; j++)
-      arr[i][j] = re(arr[i][j])
+  function coerceReal (arr, N) {
+    for (let i = 0; i < N; i++) {
+      for (let j = 0; j < N; j++) {
+        arr[i][j] = re(arr[i][j])
+      }
+    }
   }
 
-
   /** @return {'number' | 'BigNumber' | 'Complex'} */
-  function coerceTypes(mat, arr, N) {
-
+  function coerceTypes (mat, arr, N) {
     /** @type {string} */
-    let type = mat.datatype()
+    const type = mat.datatype()
 
-    if (type == 'number' || type == 'BigNumber' || type == 'Complex') {
+    if (type === 'number' || type === 'BigNumber' || type === 'Complex') {
       return type
     }
 
-    let hasNumber  = false
-    let hasBig     = false
+    let hasNumber = false
+    let hasBig = false
     let hasComplex = false
 
-    for (let i = 0; i < N; i++)
-    for (let j = 0; j < N; j++)
-    {
-      const el = arr[i][j]
-
-      if (isNumber(el) || isFraction(el))
-        hasNumber = true
-
-      else if (isBigNumber(el))
-        hasBig = true
-
-      else if (isComplex(el))
-        hasComplex = true
-
-      else
-        throw TypeError("Unsupported type in Matrix: " + typeOf(el))
+    for (let i = 0; i < N; i++) {
+      for (let j = 0; j < N; j++) {
+        const el = arr[i][j]
+
+        if (isNumber(el) || isFraction(el)) {
+          hasNumber = true
+        } else if (isBigNumber(el)) {
+          hasBig = true
+        } else if (isComplex(el)) {
+          hasComplex = true
+        } else {
+          throw TypeError('Unsupported type in Matrix: ' + typeOf(el))
+        }
+      }
     }
 
-    if (hasBig && hasComplex)
-      console.warn("Complex BigNumbers not supported, this operation will lose precission.");
+    if (hasBig && hasComplex) {
+      console.warn('Complex BigNumbers not supported, this operation will lose precission.')
+    }
 
-    if (hasComplex)
-    {
-      for (let i = 0; i < N; i++)
-      for (let j = 0; j < N; j++)
-        arr[i][j] = complex(arr[i][j])
+    if (hasComplex) {
+      for (let i = 0; i < N; i++) {
+        for (let j = 0; j < N; j++) {
+          arr[i][j] = complex(arr[i][j])
+        }
+      }
 
       return 'Complex'
     }
 
-    if (hasBig)
-    {
-      for (let i = 0; i < N; i++)
-      for (let j = 0; j < N; j++)
-        arr[i][j] = bignumber(arr[i][j])
+    if (hasBig) {
+      for (let i = 0; i < N; i++) {
+        for (let j = 0; j < N; j++) {
+          arr[i][j] = bignumber(arr[i][j])
+        }
+      }
 
       return 'BigNumber'
     }
 
-    if (hasNumber)
-    {
-      for (let i = 0; i < N; i++)
-      for (let j = 0; j < N; j++)
-        arr[i][j] = number(arr[i][j])
+    if (hasNumber) {
+      for (let i = 0; i < N; i++) {
+        for (let j = 0; j < N; j++) {
+          arr[i][j] = number(arr[i][j])
+        }
+      }
 
       return 'number'
+    } else {
+      throw TypeError('Matrix contains unsupported types only.')
     }
-
-    else throw TypeError('Matrix contains unsupported types only.')
   }
 
   return eigs
-})
\ No newline at end of file
+})
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 30f6e72119..20728a1ca7 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,534 +1,519 @@
 import { clone } from '../../../utils/object'
 
-export function createComplex({addScalar, subtract, multiply, sqrt, abs, bignumber, diag, inv, qr, usolve}) {
-    /**
-     * @param {number[][]} arr the matrix to find eigenvalues of
-     * @param {number} N size of the matrix
-     * @param {number|BigNumber} prec precision, anything lower will be considered zero
-     * @param {'number'|'BigNumber'|'Complex'} type
-     * @param {boolean} findVectors should we find eigenvectors?
-     */
-    function main(arr, N, prec, type, findVectors)
-    {
-        if (findVectors === undefined) findVectors = true
-
-        // TODO implement QR for complex matrices
-        if (type === 'Complex')
-            throw new TypeError('Complex matrices not yet supported')
-
+export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, inv, qr, usolve }) {
+  /**
+   * @param {number[][]} arr the matrix to find eigenvalues of
+   * @param {number} N size of the matrix
+   * @param {number|BigNumber} prec precision, anything lower will be considered zero
+   * @param {'number'|'BigNumber'|'Complex'} type
+   * @param {boolean} findVectors should we find eigenvectors?
+   */
+  function main (arr, N, prec, type, findVectors) {
+    if (findVectors === undefined) {
+      findVectors = true
+    }
 
-        // TODO check if any row/col are zero except the diagonal
+    // TODO implement QR for complex matrices
+    if (type === 'Complex') {
+      throw new TypeError('Complex matrices not yet supported')
+    }
 
-        // make sure corresponding rows and columns have similar magnitude
-        // important because of numerical stability
-        let R = balance(arr, N, prec, type, findVectors)
+    // TODO check if any row/col are zero except the diagonal
 
-        // R is the row transformation matrix
-        // A' = R A R⁻¹, A is the original matrix
-        // (if findVectors is false, R is undefined)
+    // make sure corresponding rows and columns have similar magnitude
+    // important because of numerical stability
+    const R = balance(arr, N, prec, type, findVectors)
 
-        // TODO if magnitudes of elements vary over many orders,
-        // move greatest elements to the top left corner
+    // R is the row transformation matrix
+    // A' = R A R⁻¹, A is the original matrix
+    // (if findVectors is false, R is undefined)
 
-        // using similarity transformations, reduce the matrix
-        // to Hessenberg form (upper triangular plus one subdiagonal row)
-        // updates the transformation matrix R with new row operationsq
-        reduceToHessenberg(arr, N, prec, type, findVectors, R)
+    // TODO if magnitudes of elements vary over many orders,
+    // move greatest elements to the top left corner
 
-        // find eigenvalues
-        let { values, C } = iterateUntilTriangular(arr, N, prec, type, findVectors)
+    // using similarity transformations, reduce the matrix
+    // to Hessenberg form (upper triangular plus one subdiagonal row)
+    // updates the transformation matrix R with new row operationsq
+    reduceToHessenberg(arr, N, prec, type, findVectors, R)
 
-        // values is the list of eigenvalues, C is the column
-        // transformation matrix that transforms the hessenberg
-        // matrix to upper triangular
+    // find eigenvalues
+    let { values, C } = iterateUntilTriangular(arr, N, prec, type, findVectors)
 
-        // compose transformations A → hess. and hess. → triang.
-        C = multiply( inv(R), C )
+    // values is the list of eigenvalues, C is the column
+    // transformation matrix that transforms the hessenberg
+    // matrix to upper triangular
 
+    // compose transformations A → hess. and hess. → triang.
+    C = multiply(inv(R), C)
 
-        let vectors
+    let vectors
 
-        if (findVectors)
-            vectors = findEigenvectors(arr, N, C, values)
+    if (findVectors) {
+      vectors = findEigenvectors(arr, N, C, values)
+    }
 
-        return { values, vectors }
+    return { values, vectors }
+  }
+
+  /**
+   * @param {number[][]} arr
+   * @param {number} N
+   * @param {number} prec
+   * @param {'number'|'BigNumber'|'Complex'} type
+   * @returns {number[][]}
+   */
+  function balance (arr, N, prec, type, findVectors) {
+    const big = type === 'BigNumber'
+    const cplx = type === 'Complex'
+
+    const bigZero = bignumber(0)
+    const bigOne = bignumber(1)
+
+    // base of the floating-point arithmetic
+    const radix = big ? bignumber(10) : 2
+    const radixSq = multiply(radix, radix)
+
+    // the diagonal transformation matrix R
+    let Rdiag
+    if (findVectors) {
+      Rdiag = Array(N).fill(big ? bigOne : 1)
     }
 
-    /**
-     * @param {number[][]} arr
-     * @param {number} N
-     * @param {number} prec
-     * @param {'number'|'BigNumber'|'Complex'} type
-     * @returns {number[][]}
-     */
-    function balance(arr, N, prec, type, findVectors)
-    {
-        const big  = type === 'BigNumber'
-        const cplx = type === 'Complex'
-
-        const bigZero = bignumber(0)
-        const bigOne  = bignumber(1)
-
-        // base of the floating-point arithmetic
-        const radix = big ? bignumber(10) : 2
-        const radixSq = multiply(radix, radix)
-
-        // the diagonal transformation matrix R
-        let Rdiag
-        if (findVectors)
-            Rdiag = Array(N).fill(big ? bigOne : 1)
-
-        // this isn't the only time we loop thru the matrix...
-        let last = false
-
-        while (!last)
-        {
-            // ...haha I'm joking! unless...
-            last = true
-
-            for (let i = 0; i < N; i++)
-            {
-
-                // compute the taxicab norm of i-th column and row
-                // TODO optimize for complex numbers
-                let colNorm = big ? bigZero : 0
-                let rowNorm = big ? bigZero : 0
-
-                for (let j = 0; j < N; j++)
-                {
-                    if (i === j) continue;
-                    let c = big||cplx ? arr[i][j].abs() : Math.abs(arr[i][j])
-                    let r = big||cplx ? arr[j][i].abs() : Math.abs(arr[j][i])
-                    colNorm = big ? colNorm.add(c) : colNorm + c
-                    rowNorm = big ? rowNorm.add(c) : rowNorm + r
-                }
-
-
-                if (!big ? (colNorm !== 0 && rowNorm !== 0) : (!colNorm.equals(0) && !rowNorm.equals(0)) )
-                {
-                    // find integer power closest to balancing the matrix
-                    // (we want to scale only by integer powers of radix,
-                    // so that we don't lose any precision due to round-off)
-
-
-                    let f = big ? bigOne : 1
-                    let c = colNorm
-
-                    const rowDivRadix = big ? rowNorm.div(radix) : rowNorm/radix
-                    const rowMulRadix = big ? rowNorm.mul(radix) : rowNorm*radix
-
-                    if (!big) {
-                        while (c < rowDivRadix) {
-                            c *= radixSq
-                            f *= radix
-                        }
-                        while (c > rowMulRadix) {
-                            c /= radixSq
-                            f /= radix
-                        }
-                    } else {
-                        while (c.lessThan(rowDivRadix)) {
-                            c = c.mul(radixSq)
-                            f = f.mul(radix)
-                        }
-                        while (c.greaterThan(rowMulRadix)) {
-                            c = c.div(radixSq)
-                            f = f.div(radix)
-                        }
-                    }
-
-
-                    // check whether balancing is needed
-                    const condition = !big
-                        ? (c + rowNorm)/f < 0.95*(colNorm + rowNorm)
-                        : c.add(rowNorm).div(f).lessThan(colNorm.add(rowNorm).mul(0.95))
-
-
-                    // apply balancing similarity transformation
-                    if (condition)
-                    {
-                        // we should loop once again to check whether
-                        // another rebalancing is needed
-                        last = false
-
-                        let g = big ? bigOne.div(f) : 1/f
-
-                        for (let j = 0; j < N; j++)
-                        {
-                            if (i === j) continue;
-                            arr[i][j] = multiply(arr[i][j], f)
-                            arr[j][i] = multiply(arr[j][i], g)
-                        }
-
-                        // keep track of transformations
-                        if (findVectors)
-                            Rdiag[i] = big ? R[i].mul(f) : Rdiag[i]*f
-                    }
-
-                }
-            }
+    // this isn't the only time we loop thru the matrix...
+    let last = false
+
+    while (!last) {
+      // ...haha I'm joking! unless...
+      last = true
+
+      for (let i = 0; i < N; i++) {
+        // compute the taxicab norm of i-th column and row
+        // TODO optimize for complex numbers
+        let colNorm = big ? bigZero : 0
+        let rowNorm = big ? bigZero : 0
+
+        for (let j = 0; j < N; j++) {
+          if (i === j) continue
+          const c = (big || cplx) ? arr[i][j].abs() : Math.abs(arr[i][j])
+          const r = (big || cplx) ? arr[j][i].abs() : Math.abs(arr[j][i])
+          colNorm = big ? colNorm.add(c) : colNorm + c
+          rowNorm = big ? rowNorm.add(c) : rowNorm + r
         }
 
-        // return the diagonal row transformation matrix
-        return diag(Rdiag)
-    }
+        if (!big ? (colNorm !== 0 && rowNorm !== 0) : (!colNorm.equals(0) && !rowNorm.equals(0))) {
+          // find integer power closest to balancing the matrix
+          // (we want to scale only by integer powers of radix,
+          // so that we don't lose any precision due to round-off)
 
+          let f = big ? bigOne : 1
+          let c = colNorm
 
-    /**
-     * @param {number[][]} arr
-     * @param {number} N
-     * @param {number} prec
-     * @param {'number'|'BigNumber'|'Complex'} type
-     * @param {boolean} findVectors
-     * @param {number[][]} R the row transformation matrix that will be modified
-     */
-    function reduceToHessenberg(arr, N, prec, type, findVectors, R)
-    {
-        const big  = type === 'BigNumber'
-        const cplx = type === 'Complex'
-        const num = !big && !cplx
-
-        const bigZero = bignumber(0)
-
-        for (let i = 0; i < N-2; i++)
-        {
-            // Find the largest subdiag element in the i-th col
-
-            let maxIndex = 0
-            let max = big ? bigZero : 0
-
-            for (let j = i + 1; j < N; j++)
-            {
-                let el = arr[j][i]
-                if (big ? abs(max).lessThan(abs(el)) : abs(max) < abs(el)) {
-                    max = el
-                    maxIndex = j
-                }
-            }
-
-            // This col is pivoted, no need to do anything
-            if (big ? abs(max).lessThan(prec) : abs(max) < prec)
-                continue
+          const rowDivRadix = big ? rowNorm.div(radix) : rowNorm / radix
+          const rowMulRadix = big ? rowNorm.mul(radix) : rowNorm * radix
 
-            if (maxIndex !== i+1)
-            {
-                // Interchange maxIndex-th and (i+1)-th row
-                const tmp1 = arr[maxIndex]
-                arr[maxIndex] = arr[i+1]
-                arr[i+1] = tmp1
-
-                // Interchange maxIndex-th and (i+1)-th column
-                for (let j = 0; j < N; j++)
-                {
-                    const tmp2 = arr[j][maxIndex]
-                    arr[j][maxIndex] = arr[j][i+1]
-                    arr[j][i+1] = tmp2
-                }
-
-                // keep track of transformations
-                if (findVectors) {
-                    const tmp3 = R[maxIndex]
-                    R[maxIndex] = R[i+1]
-                    R[i+1] = tmp3
-                }
+          if (!big) {
+            while (c < rowDivRadix) {
+              c *= radixSq
+              f *= radix
+            }
+            while (c > rowMulRadix) {
+              c /= radixSq
+              f /= radix
+            }
+          } else {
+            while (c.lessThan(rowDivRadix)) {
+              c = c.mul(radixSq)
+              f = f.mul(radix)
+            }
+            while (c.greaterThan(rowMulRadix)) {
+              c = c.div(radixSq)
+              f = f.div(radix)
             }
+          }
 
-            // Reduce following rows and columns
-            for (let j = i + 2; j < N; j++)
-            {
-                let n = num ? arr[j][i] / max : arr[j][i].div(max)
+          // check whether balancing is needed
+          const condition = !big
+            ? (c + rowNorm) / f < 0.95 * (colNorm + rowNorm)
+            : c.add(rowNorm).div(f).lessThan(colNorm.add(rowNorm).mul(0.95))
 
-                if (n == 0) continue;
+          // apply balancing similarity transformation
+          if (condition) {
+            // we should loop once again to check whether
+            // another rebalancing is needed
+            last = false
 
-                // from j-th row subtract n-times (i+1)th row
-                for (let k = 0; k < N; k++)
-                    arr[j][k] = num ? arr[j][k] - n*arr[i+1][k] : arr[j][k].sub(n.mul(arr[i+1][k]))
+            const g = big ? bigOne.div(f) : 1 / f
 
-                // to (i+1)th column add n-times j-th column
-                for (let k = 0; k < N; k++)
-                    arr[k][i+1] = num ? arr[k][i+1] + n*arr[k][j] : arr[k][i+1].add(n.mul(arr[k][j]))
+            for (let j = 0; j < N; j++) {
+              if (i === j) {
+                continue
+              }
+              arr[i][j] = multiply(arr[i][j], f)
+              arr[j][i] = multiply(arr[j][i], g)
+            }
 
-                // keep track of transformations
-                if (findVectors)
-                    for (let k = 0; k < N; k++)
-                        R[j][k] = num ? R[j][k] - n*R[i+1][k] : subtract(R[j][k], n.mul(R[i+1][k]))
+            // keep track of transformations
+            if (findVectors) {
+              Rdiag[i] = big ? Rdiag[i].mul(f) : Rdiag[i] * f
             }
+          }
         }
-
-        return R
+      }
     }
 
-    /**
-     * @returns {{values: values, C: Matrix}}
-     * @see Press, Wiliams: Numerical recipes in Fortran 77
-     * @see https://en.wikipedia.org/wiki/QR_algorithm
-     */
-    function iterateUntilTriangular(A, N, prec, type, findVectors)
-    {
-        const big = type === 'BigNumber'
-
-        // The Francis Algorithm
-        // The core idea of this algorithm is that doing successive
-        // A' = Q⁺AQ transformations will eventually converge to block-
-        // upper-triangular with diagonal blocks either 1x1 or 2x2.
-        // The Q here is the one from the QR decomposition, A = QR.
-        // Since the eigenvalues of a block-upper-triangular matrix are
-        // the eigenvalues of its diagonal blocks and we know how to find
-        // eigenvalues of a 2x2 matrix, we know the eigenvalues of A.
-
+    // return the diagonal row transformation matrix
+    return diag(Rdiag)
+  }
+
+  /**
+   * @param {number[][]} arr
+   * @param {number} N
+   * @param {number} prec
+   * @param {'number'|'BigNumber'|'Complex'} type
+   * @param {boolean} findVectors
+   * @param {number[][]} R the row transformation matrix that will be modified
+   */
+  function reduceToHessenberg (arr, N, prec, type, findVectors, R) {
+    const big = type === 'BigNumber'
+    const cplx = type === 'Complex'
+    const num = !big && !cplx
+
+    const bigZero = bignumber(0)
+
+    for (let i = 0; i < N - 2; i++) {
+      // Find the largest subdiag element in the i-th col
+
+      let maxIndex = 0
+      let max = big ? bigZero : 0
+
+      for (let j = i + 1; j < N; j++) {
+        const el = arr[j][i]
+        if (big ? abs(max).lessThan(abs(el)) : abs(max) < abs(el)) {
+          max = el
+          maxIndex = j
+        }
+      }
+
+      // This col is pivoted, no need to do anything
+      if (big ? abs(max).lessThan(prec) : abs(max) < prec) {
+        continue
+      }
+
+      if (maxIndex !== i + 1) {
+        // Interchange maxIndex-th and (i+1)-th row
+        const tmp1 = arr[maxIndex]
+        arr[maxIndex] = arr[i + 1]
+        arr[i + 1] = tmp1
+
+        // Interchange maxIndex-th and (i+1)-th column
+        for (let j = 0; j < N; j++) {
+          const tmp2 = arr[j][maxIndex]
+          arr[j][maxIndex] = arr[j][i + 1]
+          arr[j][i + 1] = tmp2
+        }
 
-        let arr = clone(A)
+        // keep track of transformations
+        if (findVectors) {
+          const tmp3 = R[maxIndex]
+          R[maxIndex] = R[i + 1]
+          R[i + 1] = tmp3
+        }
+      }
 
-        // the list of converged eigenvalues
-        const lambdas = []
+      // Reduce following rows and columns
+      for (let j = i + 2; j < N; j++) {
+        const n = num ? arr[j][i] / max : arr[j][i].div(max)
 
-        // size of arr, which will get smaller as eigenvalues converge
-        let n = N;
+        if (n === 0) {
+          continue
+        }
 
-        // the diagonal of the block-diagonal matrix that turns
-        // converged 2x2 matrices into upper triangular matrices
-        const Sdiag = []
+        // from j-th row subtract n-times (i+1)th row
+        for (let k = 0; k < N; k++) {
+          arr[j][k] = num ? arr[j][k] - n * arr[i + 1][k] : arr[j][k].sub(n.mul(arr[i + 1][k]))
+        }
 
-        // N×N matrix describing the overall transformation done during the QR algorithm
-        let Qtotal = findVectors ? diag(Array(N).fill(1)) : undefined
+        // to (i+1)th column add n-times j-th column
+        for (let k = 0; k < N; k++) {
+          arr[k][i + 1] = num ? arr[k][i + 1] + n * arr[k][j] : arr[k][i + 1].add(n.mul(arr[k][j]))
+        }
 
-        // n×n matrix describing the QR transformations done since last convergence
-        let Qpartial = findVectors ? diag(Array(n).fill(1)) : undefined
+        // keep track of transformations
+        if (findVectors) {
+          for (let k = 0; k < N; k++) {
+            R[j][k] = num ? R[j][k] - n * R[i + 1][k] : subtract(R[j][k], n.mul(R[i + 1][k]))
+          }
+        }
+      }
+    }
 
-        // last eigenvalue converged before this many steps
-        let lastConvergenceBefore = 0;
+    return R
+  }
 
+  /**
+   * @returns {{values: values, C: Matrix}}
+   * @see Press, Wiliams: Numerical recipes in Fortran 77
+   * @see https://en.wikipedia.org/wiki/QR_algorithm
+   */
+  function iterateUntilTriangular (A, N, prec, type, findVectors) {
+    const big = type === 'BigNumber'
 
-        while (lastConvergenceBefore <= 30)
-        {
-            lastConvergenceBefore += 1
+    // The Francis Algorithm
+    // The core idea of this algorithm is that doing successive
+    // A' = Q⁺AQ transformations will eventually converge to block-
+    // upper-triangular with diagonal blocks either 1x1 or 2x2.
+    // The Q here is the one from the QR decomposition, A = QR.
+    // Since the eigenvalues of a block-upper-triangular matrix are
+    // the eigenvalues of its diagonal blocks and we know how to find
+    // eigenvalues of a 2x2 matrix, we know the eigenvalues of A.
 
-            // TODO if the convergence is slow, do something clever
+    let arr = clone(A)
 
-            // Perform the factorization
+    // the list of converged eigenvalues
+    const lambdas = []
 
-            let k = 0; // TODO set close to an eigenvalue
+    // size of arr, which will get smaller as eigenvalues converge
+    let n = N
 
-            for (let i = 0; i < n; i++)
-                arr[i][i] -= k
+    // the diagonal of the block-diagonal matrix that turns
+    // converged 2x2 matrices into upper triangular matrices
+    const Sdiag = []
 
-            // TODO do an implicit QR transformation
-            let {Q, R} = qr(arr)
-            arr = multiply(R, Q)
+    // N×N matrix describing the overall transformation done during the QR algorithm
+    let Qtotal = findVectors ? diag(Array(N).fill(1)) : undefined
 
-            for (let i = 0; i < n; i++)
-                arr[i][i] += k
+    // n×n matrix describing the QR transformations done since last convergence
+    let Qpartial = findVectors ? diag(Array(n).fill(1)) : undefined
 
+    // last eigenvalue converged before this many steps
+    let lastConvergenceBefore = 0
 
-            // keep track of transformations
-            if (findVectors)
-                Qpartial = multiply(Qpartial, Q)
-
-
-            // The rightmost diagonal element converged to an eigenvalue
-            if (n == 1 || (big ? abs(arr[n-1][n-2]).lessThan(prec) : abs(arr[n-1][n-2]) < prec))
-            {
-                lastConvergenceBefore = 0
-                lambdas.push(arr[n-1][n-1])
-
-                // keep track of transformations
-                if (findVectors) {
-                    Sdiag.unshift( [[1]] )
-                    inflateMatrix(Qpartial, N)
-                    Qtotal = multiply(Qtotal, Qpartial)
-
-                    if (n > 1)
-                        Qpartial = diag( Array(n-1).fill(1) )
-                }
-
-                // reduce the matrix size
-                n -= 1
-                arr.pop()
-                for (let i = 0; i < n; i++)
-                    arr[i].pop()
-            }
+    while (lastConvergenceBefore <= 30) {
+      lastConvergenceBefore += 1
 
-            // The rightmost diagonal 2x2 block converged
-            else if (n == 2 || (big ? abs(arr[n-2][n-3]).lessThan(prec) : abs(arr[n-2][n-3]) < prec))
-            {
-                lastConvergenceBefore = 0
-                let ll = eigenvalues2x2(
-                    arr[n-2][n-2], arr[n-2][n-1],
-                    arr[n-1][n-2], arr[n-1][n-1]
-                )
-                lambdas.push(...ll)
-
-                // keep track of transformations
-                if (findVectors) {
-                    Sdiag.unshift(jordanBase2x2(
-                        arr[n-2][n-2], arr[n-2][n-1],
-                        arr[n-1][n-2], arr[n-1][n-1],
-                        ll[0], ll[1], prec, type
-                    ))
-                    inflateMatrix(Qpartial, N)
-                    Qtotal = multiply(Qtotal, Qpartial)
-
-                    if (n > 2)
-                        Qpartial = diag( Array(n-2).fill(1) )
-                }
-
-                // reduce the matrix size
-                n -= 2
-                arr.pop()
-                arr.pop()
-                for (let i = 0; i < n; i++) {
-                    arr[i].pop()
-                    arr[i].pop()
-                }
-            }
+      // TODO if the convergence is slow, do something clever
 
-            if (n == 0) break;
-        }
+      // Perform the factorization
 
-        // the algorithm didn't converge
-        if (lastConvergenceBefore > 30)
-            throw Error('The eigenvalues failed to converge. Only found these eigenvalues: '+values.join(', '))
+      const k = 0 // TODO set close to an eigenvalue
 
-        // combine the overall QR transformation Qtotal with the subsequent
-        // transformation S that turns the diagonal 2x2 blocks to upper triangular
-        let C = findVectors ? multiply( Qtotal, blockDiag(Sdiag, N) ) : undefined
+      for (let i = 0; i < n; i++) {
+        arr[i][i] -= k
+      }
 
-        return { values: lambdas, C }
-    }
+      // TODO do an implicit QR transformation
+      const { Q, R } = qr(arr)
+      arr = multiply(R, Q)
 
+      for (let i = 0; i < n; i++) {
+        arr[i][i] += k
+      }
 
-    /**
-     * @param {Matrix} A original matrix
-     * @param {number} N size of A
-     * @param {Matrix} C column transformation matrix that turns A into upper triangular
-     * @param {number[]} values array of eigenvalues of A
-     * @returns {Matrix[]} eigenvalues
-     */
-    function findEigenvectors(A, N, C, values)
-    {
-        const Cinv = inv(C)
-        const U = multiply( Cinv, A, C )
+      // keep track of transformations
+      if (findVectors) {
+        Qpartial = multiply(Qpartial, Q)
+      }
 
-        const b = Array(N).fill(0)
-        const E = diag(Array(N).fill(1))
+      // The rightmost diagonal element converged to an eigenvalue
+      if (n === 1 || (big ? abs(arr[n - 1][n - 2]).lessThan(prec) : abs(arr[n - 1][n - 2]) < prec)) {
+        lastConvergenceBefore = 0
+        lambdas.push(arr[n - 1][n - 1])
 
-        let vectors = []
+        // keep track of transformations
+        if (findVectors) {
+          Sdiag.unshift([[1]])
+          inflateMatrix(Qpartial, N)
+          Qtotal = multiply(Qtotal, Qpartial)
 
-        for (const l of values) {
-            // !FIXME this is a mock implementation
-            const v = usolve( subtract(U, multiply(l, E)), b )
-            vectors.push( multiply(C, v) )
+          if (n > 1) {
+            Qpartial = diag(Array(n - 1).fill(1))
+          }
         }
 
-        return vectors
-    }
+        // reduce the matrix size
+        n -= 1
+        arr.pop()
+        for (let i = 0; i < n; i++) {
+          arr[i].pop()
+        }
 
+      // The rightmost diagonal 2x2 block converged
+      } else if (n === 2 || (big ? abs(arr[n - 2][n - 3]).lessThan(prec) : abs(arr[n - 2][n - 3]) < prec)) {
+        lastConvergenceBefore = 0
+        const ll = eigenvalues2x2(
+          arr[n - 2][n - 2], arr[n - 2][n - 1],
+          arr[n - 1][n - 2], arr[n - 1][n - 1]
+        )
+        lambdas.push(...ll)
+
+        // keep track of transformations
+        if (findVectors) {
+          Sdiag.unshift(jordanBase2x2(
+            arr[n - 2][n - 2], arr[n - 2][n - 1],
+            arr[n - 1][n - 2], arr[n - 1][n - 1],
+            ll[0], ll[1], prec, type
+          ))
+          inflateMatrix(Qpartial, N)
+          Qtotal = multiply(Qtotal, Qpartial)
+
+          if (n > 2) {
+            Qpartial = diag(Array(n - 2).fill(1))
+          }
+        }
 
-    /**
-     * Compute the eigenvalues of an 2x2 matrix
-     * @return {[number,number]}
-     */
-    function eigenvalues2x2(a,b,c,d)
-    {
-        // λ± = ½ trA ± ½ √( tr²A - 4 detA )
-        const trA = addScalar(a, d)
-        const detA = subtract(multiply(a, d), multiply(b, c))
-        const x = multiply(trA, 0.5)
-        const y = multiply(sqrt( subtract(multiply(trA, trA), multiply(4, detA)) ), 0.5)
+        // reduce the matrix size
+        n -= 2
+        arr.pop()
+        arr.pop()
+        for (let i = 0; i < n; i++) {
+          arr[i].pop()
+          arr[i].pop()
+        }
+      }
 
-        return [addScalar(x, y), subtract(x, y)]
+      if (n === 0) {
+        break
+      }
     }
 
+    // the algorithm didn't converge
+    if (lastConvergenceBefore > 30) {
+      throw Error('The eigenvalues failed to converge. Only found these eigenvalues: ' + lambdas.join(', '))
+    }
 
-    /**
-     * For an 2x2 matrix compute the transformation matrix S,
-     * so that SAS⁻¹ is an upper triangular matrix
-     * @return {[[number,number],[number,number]]}
-     * @see https://math.berkeley.edu/~ogus/old/Math_54-05/webfoils/jordan.pdf
-     * @see http://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html
-     */
-    function jordanBase2x2(a,b,c,d,l1,l2,prec,type)
-    {
-        const big = type === 'BigNumber'
-
-        const b0 = big ? 0 : bignumber(0)
-        const b1 = big ? 1 : bignumber(1)
-
-
-        // matrix is already upper triangular
-        // return an identity matrix
-        if (big ? abs(c).lessThan(prec) : abs(c) < prec)
-            return [ [b1, b0], [b0, b1] ]
-
-        // matrix is diagonalizable
-        // return its eigenvectors as columns
-        if (big ? abs(l1.sub(l2)).greaterThan(prec) : abs(subtract(l1, l2)) > prec)
-            return [ [ subtract(l1, d), subtract(l2, d) ], [ c, c ] ]
-
-
-        // matrix is not diagonalizable
-        // compute off-diagonal elements of N = A - λI
-        // N₁₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ )
-        // N₁₂ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ )
+    // combine the overall QR transformation Qtotal with the subsequent
+    // transformation S that turns the diagonal 2x2 blocks to upper triangular
+    const C = findVectors ? multiply(Qtotal, blockDiag(Sdiag, N)) : undefined
+
+    return { values: lambdas, C }
+  }
+
+  /**
+   * @param {Matrix} A original matrix
+   * @param {number} N size of A
+   * @param {Matrix} C column transformation matrix that turns A into upper triangular
+   * @param {number[]} values array of eigenvalues of A
+   * @returns {Matrix[]} eigenvalues
+   */
+  function findEigenvectors (A, N, C, values) {
+    const Cinv = inv(C)
+    const U = multiply(Cinv, A, C)
+
+    const b = Array(N).fill(0)
+    const E = diag(Array(N).fill(1))
+
+    const vectors = []
+
+    for (const l of values) {
+      // !FIXME this is a mock implementation
+      const v = usolve(subtract(U, multiply(l, E)), b)
+      vectors.push(multiply(C, v))
+    }
 
-        let na = subtract(a, l1)
-        let nb = subtract(b, l1)
-        let nc = subtract(c, l1)
-        let nd = subtract(d, l1)
+    return vectors
+  }
+
+  /**
+   * Compute the eigenvalues of an 2x2 matrix
+   * @return {[number,number]}
+   */
+  function eigenvalues2x2 (a, b, c, d) {
+    // λ± = ½ trA ± ½ √( tr²A - 4 detA )
+    const trA = addScalar(a, d)
+    const detA = subtract(multiply(a, d), multiply(b, c))
+    const x = multiply(trA, 0.5)
+    const y = multiply(sqrt(subtract(multiply(trA, trA), multiply(4, detA))), 0.5)
+
+    return [addScalar(x, y), subtract(x, y)]
+  }
+
+  /**
+   * For an 2x2 matrix compute the transformation matrix S,
+   * so that SAS⁻¹ is an upper triangular matrix
+   * @return {[[number,number],[number,number]]}
+   * @see https://math.berkeley.edu/~ogus/old/Math_54-05/webfoils/jordan.pdf
+   * @see http://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html
+   */
+  function jordanBase2x2 (a, b, c, d, l1, l2, prec, type) {
+    const big = type === 'BigNumber'
+
+    const b0 = big ? 0 : bignumber(0)
+    const b1 = big ? 1 : bignumber(1)
+
+    // matrix is already upper triangular
+    // return an identity matrix
+    if (big ? abs(c).lessThan(prec) : abs(c) < prec) {
+      return [[b1, b0], [b0, b1]]
+    }
 
-        if (big ? abs(nb).lessThan(prec) : abs(nb) < prec)
-            return [ [na, b1], [nc, b0] ]
-        else
-            return [ [nb, b0], [nd, b1] ]
+    // matrix is diagonalizable
+    // return its eigenvectors as columns
+    if (big ? abs(l1.sub(l2)).greaterThan(prec) : abs(subtract(l1, l2)) > prec) {
+      return [[subtract(l1, d), subtract(l2, d)], [c, c]]
     }
 
+    // matrix is not diagonalizable
+    // compute off-diagonal elements of N = A - λI
+    // N₁₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ )
+    // N₁₂ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ )
 
-    /**
-     * Enlarge the matrix from n×n to N×N, setting the new
-     * elements to 1 on diagonal and 0 elsewhere
-     */
-    function inflateMatrix(arr, N)
-    {
-        // add columns
-        for (let i = 0; i < arr.length; i++)
-            arr[i].push( ...Array(N-arr[i].length).fill(0) )
-
-        // add rows
-        for (let i = arr.length; i < N; i++) {
-            arr.push( Array(N).fill(0) )
-            arr[i][i] = 1
-        }
+    const na = subtract(a, l1)
+    const nb = subtract(b, l1)
+    const nc = subtract(c, l1)
+    const nd = subtract(d, l1)
 
-        return arr
+    if (big ? abs(nb).lessThan(prec) : abs(nb) < prec) {
+      return [[na, b1], [nc, b0]]
+    } else {
+      return [[nb, b0], [nd, b1]]
+    }
+  }
+
+  /**
+   * Enlarge the matrix from n×n to N×N, setting the new
+   * elements to 1 on diagonal and 0 elsewhere
+   */
+  function inflateMatrix (arr, N) {
+    // add columns
+    for (let i = 0; i < arr.length; i++) {
+      arr[i].push(...Array(N - arr[i].length).fill(0))
     }
 
+    // add rows
+    for (let i = arr.length; i < N; i++) {
+      arr.push(Array(N).fill(0))
+      arr[i][i] = 1
+    }
 
-    /**
-     * Create a block-diagonal matrix with the given square matrices on the diagonal
-     * @param {Matrix[] | number[][][]} arr array of matrices to be placed on the diagonal
-     * @param {number} N the size of the resulting matrix
-     */
-    function blockDiag(arr, N)
-    {
-        const M = []
-        for (let i = 0; i < N; i++)
-            M[i] = Array(N).fill(0)
-
-        let I = 0;
-        for (const sub of arr)
-        {
-            const n = sub.length
+    return arr
+  }
+
+  /**
+   * Create a block-diagonal matrix with the given square matrices on the diagonal
+   * @param {Matrix[] | number[][][]} arr array of matrices to be placed on the diagonal
+   * @param {number} N the size of the resulting matrix
+   */
+  function blockDiag (arr, N) {
+    const M = []
+    for (let i = 0; i < N; i++) {
+      M[i] = Array(N).fill(0)
+    }
 
-            for (let i = 0; i < n; i++)
-            for (let j = 0; j < n; j++)
-                M[I+i][I+j] = sub[i][j]
+    let I = 0
+    for (const sub of arr) {
+      const n = sub.length
 
-            I += n
+      for (let i = 0; i < n; i++) {
+        for (let j = 0; j < n; j++) {
+          M[I + i][I + j] = sub[i][j]
         }
+      }
 
-        return M
+      I += n
     }
 
+    return M
+  }
 
-    return main;
-}
\ No newline at end of file
+  return main
+}
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index dd160303ca..1277912d15 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -1,10 +1,6 @@
 import { clone } from '../../../utils/object'
 
-
-export function createRealSymmetric({ addScalar, subtract, column, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
-
-
-
+export function createRealSymmetric ({ addScalar, subtract, column, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
   /**
    * @param {number[] | BigNumber[]} arr
    * @param {number} N
@@ -12,12 +8,13 @@ export function createRealSymmetric({ addScalar, subtract, column, abs, atan, co
    * @param {'number' | 'BigNumber'} type
    */
   function main (arr, N, prec = 1E-12, type) {
-
-    if (type === 'number')
+    if (type === 'number') {
       return diag(arr, prec)
+    }
 
-    if (type === 'BigNumber')
+    if (type === 'BigNumber') {
       return diagBig(arr, prec)
+    }
 
     throw TypeError('Unsupported data type: ' + type)
   }
@@ -267,8 +264,9 @@ export function createRealSymmetric({ addScalar, subtract, column, abs, atan, co
     const values = clone(Ef)
     const vectors = []
 
-    for (let i = 0; i < Sf.length; i++)
+    for (let i = 0; i < Sf.length; i++) {
       vectors[i] = column(Sf, i)
+    }
 
     // !FIXME vectors are always Array[], never Matrix[]
     return { values, vectors }
@@ -291,6 +289,5 @@ export function createRealSymmetric({ addScalar, subtract, column, abs, atan, co
     return array
   }
 
-  return main;
-
+  return main
 }

From 5adea51b728427c7580669dbc9de337a26fbac9b Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Wed, 11 Mar 2020 16:50:07 +0100
Subject: [PATCH 12/19] some minor changes

---
 src/function/matrix/eigs.js                  | 7 ++++---
 src/function/matrix/eigs/realSymetric.js     | 4 ++--
 test/unit-tests/function/matrix/eigs.test.js | 5 +----
 3 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index d22faad013..6462f9b6cd 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -29,10 +29,11 @@ const dependencies = [
   'sin',
   'inv',
   'qr',
-  'usolve'
+  'usolve',
+  'flatten'
 ]
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, qr, usolve, sqrt, matrix, column, diag, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed, re, im, qr, usolve, flatten, sqrt, matrix, column, diag, addScalar, subtract, equal, abs, larger, atan, cos, sin, multiplyScalar, inv, complex, number, bignumber, multiply, add }) => {
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
    *
@@ -78,7 +79,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     }
   })
 
-  const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
+  const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
   const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr, inv, usolve })
 
   function computeValuesAndVectors (mat, prec) {
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index 1277912d15..e3373904fa 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -1,6 +1,6 @@
 import { clone } from '../../../utils/object'
 
-export function createRealSymmetric ({ addScalar, subtract, column, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
+export function createRealSymmetric ({ addScalar, subtract, column, flatten, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
   /**
    * @param {number[] | BigNumber[]} arr
    * @param {number} N
@@ -265,7 +265,7 @@ export function createRealSymmetric ({ addScalar, subtract, column, abs, atan, c
     const vectors = []
 
     for (let i = 0; i < Sf.length; i++) {
-      vectors[i] = column(Sf, i)
+      vectors[i] = flatten(column(Sf, i))
     }
 
     // !FIXME vectors are always Array[], never Matrix[]
diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js
index a1dcf6c8c0..f99b0b386b 100644
--- a/test/unit-tests/function/matrix/eigs.test.js
+++ b/test/unit-tests/function/matrix/eigs.test.js
@@ -11,12 +11,9 @@ describe('eigs', function () {
     assert.throws(function () { eigs([4, 5, 6]) }, /Matrix must be square/)
     assert.throws(function () { eigs(1.0) }, /TypeError: Unexpected type of argument/)
     assert.throws(function () { eigs('random') }, /TypeError: Unexpected type of argument/)
-    assert.throws(function () { eigs(math.matrix([[1, 2], [2.1, 3]])) }, /Input matrix is not symmetric/)
   })
   it('should only accept a matrix with valid element type', function () {
-    assert.throws(function () { eigs([['x', 2], [4, 5]]) }, /Mixed matrix element type is not supported/)
-    assert.throws(function () { eigs([[1, 2], [2, '5']]) }, /Mixed matrix element type is not supported/)
-    assert.throws(function () { eigs([['1', '2'], ['2', '5']]) }, /Matrix element type not supported/)
+    assert.throws(function () { eigs([['x', 2], [4, 5]]) }, /Cannot convert "x" to a number/)
   })
   it('eigenvalue check for diagonal matrix', function () {
     // trivial test

From 58641c900958e95d47cb5fcc3ee9d77daa34da50 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Sat, 26 Sep 2020 18:25:51 +0200
Subject: [PATCH 13/19] solve merge conflicts

---
 src/function/matrix/eigs.js         | 10 +++++-----
 src/function/matrix/eigs/complex.js |  8 ++++----
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 8b31cffffc..40fa4e834b 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -6,9 +6,12 @@ import { createComplex } from './eigs/complex'
 
 const name = 'eigs'
 
-const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'inv', 'bignumber', 'multiply', 'add']
+const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'complex', 'sqrt', 'diag', 'qr', 'usolveAll', 'im', 're']
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add, larger, column, flatten, complex, sqrt, diag, qr, usolveAll, im, re }) => {
+
+  const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
+  const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr, inv, usolveAll })
 
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
@@ -58,9 +61,6 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
     }
   })
 
-  const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr, inv, usolve })
-
   function computeValuesAndVectors (mat, prec) {
     if (prec === undefined) {
       prec = 1e-14
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 20728a1ca7..cc4e4f2f34 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,6 +1,6 @@
 import { clone } from '../../../utils/object'
 
-export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, inv, qr, usolve }) {
+export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, inv, qr, usolveAll }) {
   /**
    * @param {number[][]} arr the matrix to find eigenvalues of
    * @param {number} N size of the matrix
@@ -405,9 +405,9 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
     const vectors = []
 
     for (const l of values) {
-      // !FIXME this is a mock implementation
-      const v = usolve(subtract(U, multiply(l, E)), b)
-      vectors.push(multiply(C, v))
+      // ! FIXME might fail for imprecise eigenvalues, replace with an iterative eigenvector algorithm
+      const V = usolveAll(subtract(U, multiply(l, E)), b)
+      for (const v of V) { vectors.push(multiply(C, v)) }
     }
 
     return vectors

From a40e8470d30a1a018fa19ef976d351e466b1f676 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Sat, 26 Sep 2020 23:30:16 +0200
Subject: [PATCH 14/19] refactored and re-fixed #1789

---
 src/function/matrix/eigs.js              |   8 +-
 src/function/matrix/eigs/complex.js      | 130 +++++++++++------------
 src/function/matrix/eigs/realSymetric.js |  34 +++---
 3 files changed, 80 insertions(+), 92 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 40fa4e834b..c5d04a63dd 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -6,12 +6,12 @@ import { createComplex } from './eigs/complex'
 
 const name = 'eigs'
 
-const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'complex', 'sqrt', 'diag', 'qr', 'usolveAll', 'im', 're']
+const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolveAll', 'im', 're', 'smaller']
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add, larger, column, flatten, complex, sqrt, diag, qr, usolveAll, im, re }) => {
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, qr, usolveAll, im, re, smaller }) => {
 
-  const doRealSymetric = createRealSymmetric({ addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplex = createComplex({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, qr, inv, usolveAll })
+  const doRealSymetric = createRealSymmetric({ config, addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
+  const doComplex = createComplex({ config, addScalar, subtract, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, qr, inv, usolveAll, equal, complex, larger, smaller })
 
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index cc4e4f2f34..b1df625a8b 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,6 +1,6 @@
 import { clone } from '../../../utils/object'
 
-export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignumber, diag, inv, qr, usolveAll }) {
+export function createComplex ({ addScalar, subtract, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, inv, qr, usolveAll, equal, complex, larger, smaller }) {
   /**
    * @param {number[][]} arr the matrix to find eigenvalues of
    * @param {number} N size of the matrix
@@ -15,7 +15,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
 
     // TODO implement QR for complex matrices
     if (type === 'Complex') {
-      throw new TypeError('Complex matrices not yet supported')
+      //throw new TypeError('Complex matrices not yet supported')
     }
 
     // TODO check if any row/col are zero except the diagonal
@@ -66,17 +66,17 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
     const big = type === 'BigNumber'
     const cplx = type === 'Complex'
 
-    const bigZero = bignumber(0)
-    const bigOne = bignumber(1)
+    const zero = big ? bignumber(0) : cplx ? complex(0) : 0
+    const one  = big ? bignumber(1) : cplx ? complex(1) : 1
 
     // base of the floating-point arithmetic
     const radix = big ? bignumber(10) : 2
-    const radixSq = multiply(radix, radix)
+    const radixSq = multiplyScalar(radix, radix)
 
     // the diagonal transformation matrix R
     let Rdiag
     if (findVectors) {
-      Rdiag = Array(N).fill(big ? bigOne : 1)
+      Rdiag = Array(N).fill(one)
     }
 
     // this isn't the only time we loop thru the matrix...
@@ -89,18 +89,18 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
       for (let i = 0; i < N; i++) {
         // compute the taxicab norm of i-th column and row
         // TODO optimize for complex numbers
-        let colNorm = big ? bigZero : 0
-        let rowNorm = big ? bigZero : 0
+        let colNorm = zero
+        let rowNorm = zero
 
         for (let j = 0; j < N; j++) {
           if (i === j) continue
-          const c = (big || cplx) ? arr[i][j].abs() : Math.abs(arr[i][j])
-          const r = (big || cplx) ? arr[j][i].abs() : Math.abs(arr[j][i])
-          colNorm = big ? colNorm.add(c) : colNorm + c
-          rowNorm = big ? rowNorm.add(c) : rowNorm + r
+          const c = abs(arr[i][j])
+          const r = abs(arr[j][i])
+          colNorm = addScalar(colNorm, c)
+          rowNorm = addScalar(rowNorm, c)
         }
 
-        if (!big ? (colNorm !== 0 && rowNorm !== 0) : (!colNorm.equals(0) && !rowNorm.equals(0))) {
+        if (!equal(colNorm, 0) && !equal(rowNorm, 0)) {
           // find integer power closest to balancing the matrix
           // (we want to scale only by integer powers of radix,
           // so that we don't lose any precision due to round-off)
@@ -108,33 +108,21 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
           let f = big ? bigOne : 1
           let c = colNorm
 
-          const rowDivRadix = big ? rowNorm.div(radix) : rowNorm / radix
-          const rowMulRadix = big ? rowNorm.mul(radix) : rowNorm * radix
+          const rowDivRadix = divideScalar(rowNorm, radix)
+          const rowMulRadix = multiplyScalar(rowNorm, radix)
 
-          if (!big) {
-            while (c < rowDivRadix) {
-              c *= radixSq
-              f *= radix
-            }
-            while (c > rowMulRadix) {
-              c /= radixSq
-              f /= radix
-            }
-          } else {
-            while (c.lessThan(rowDivRadix)) {
-              c = c.mul(radixSq)
-              f = f.mul(radix)
-            }
-            while (c.greaterThan(rowMulRadix)) {
-              c = c.div(radixSq)
-              f = f.div(radix)
-            }
+          while (smaller(c, rowDivRadix)) {
+            c = multiplyScalar(c, radixSq)
+            f = multiplyScalar(f, radix)
+          }
+          while (larger(c, rowMulRadix)) {
+            c = divideScalar(c, radixSq)
+            f = divideScalar(f, radix)
           }
 
           // check whether balancing is needed
-          const condition = !big
-            ? (c + rowNorm) / f < 0.95 * (colNorm + rowNorm)
-            : c.add(rowNorm).div(f).lessThan(colNorm.add(rowNorm).mul(0.95))
+          // condition = (c + rowNorm) / f < 0.95 * (colNorm + rowNorm)
+          const condition = smaller( divideScalar(addScalar(c, rowNorm), f), multiplyScalar(addScalar(colNorm, rowNorm), 0.95))
 
           // apply balancing similarity transformation
           if (condition) {
@@ -142,19 +130,19 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
             // another rebalancing is needed
             last = false
 
-            const g = big ? bigOne.div(f) : 1 / f
+            const g = divideScalar(1, f)
 
             for (let j = 0; j < N; j++) {
               if (i === j) {
                 continue
               }
-              arr[i][j] = multiply(arr[i][j], f)
-              arr[j][i] = multiply(arr[j][i], g)
+              arr[i][j] = multiplyScalar(arr[i][j], f)
+              arr[j][i] = multiplyScalar(arr[j][i], g)
             }
 
             // keep track of transformations
             if (findVectors) {
-              Rdiag[i] = big ? Rdiag[i].mul(f) : Rdiag[i] * f
+              Rdiag[i] = multiplyScalar(Rdiag[i], f)
             }
           }
         }
@@ -176,26 +164,25 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
   function reduceToHessenberg (arr, N, prec, type, findVectors, R) {
     const big = type === 'BigNumber'
     const cplx = type === 'Complex'
-    const num = !big && !cplx
 
-    const bigZero = bignumber(0)
+    const zero = big ? bignumber(0) : cplx ? complex(0) : 0
 
     for (let i = 0; i < N - 2; i++) {
       // Find the largest subdiag element in the i-th col
 
       let maxIndex = 0
-      let max = big ? bigZero : 0
+      let max = zero
 
       for (let j = i + 1; j < N; j++) {
         const el = arr[j][i]
-        if (big ? abs(max).lessThan(abs(el)) : abs(max) < abs(el)) {
+        if (smaller(abs(max), abs(el))) {
           max = el
           maxIndex = j
         }
       }
 
       // This col is pivoted, no need to do anything
-      if (big ? abs(max).lessThan(prec) : abs(max) < prec) {
+      if (smaller(abs(max), prec)) {
         continue
       }
 
@@ -222,7 +209,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
 
       // Reduce following rows and columns
       for (let j = i + 2; j < N; j++) {
-        const n = num ? arr[j][i] / max : arr[j][i].div(max)
+        const n = divideScalar(arr[j][i], max)
 
         if (n === 0) {
           continue
@@ -230,18 +217,18 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
 
         // from j-th row subtract n-times (i+1)th row
         for (let k = 0; k < N; k++) {
-          arr[j][k] = num ? arr[j][k] - n * arr[i + 1][k] : arr[j][k].sub(n.mul(arr[i + 1][k]))
+          arr[j][k] = subtract(arr[j][k], multiplyScalar(n, arr[i + 1][k]))
         }
 
         // to (i+1)th column add n-times j-th column
         for (let k = 0; k < N; k++) {
-          arr[k][i + 1] = num ? arr[k][i + 1] + n * arr[k][j] : arr[k][i + 1].add(n.mul(arr[k][j]))
+          arr[k][i + 1] = addScalar(arr[k][i + 1], multiplyScalar(n, arr[k][j]))
         }
 
         // keep track of transformations
         if (findVectors) {
           for (let k = 0; k < N; k++) {
-            R[j][k] = num ? R[j][k] - n * R[i + 1][k] : subtract(R[j][k], n.mul(R[i + 1][k]))
+            R[j][k] = subtract(R[j][k], multiplyScalar(n, R[i + 1][k]))
           }
         }
       }
@@ -257,6 +244,9 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
    */
   function iterateUntilTriangular (A, N, prec, type, findVectors) {
     const big = type === 'BigNumber'
+    const cplx = type === 'Complex'
+
+    const one = big ? bignumber(1) : cplx ? complex(1) : 1
 
     // The Francis Algorithm
     // The core idea of this algorithm is that doing successive
@@ -280,10 +270,10 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
     const Sdiag = []
 
     // N×N matrix describing the overall transformation done during the QR algorithm
-    let Qtotal = findVectors ? diag(Array(N).fill(1)) : undefined
+    let Qtotal = findVectors ? diag(Array(N).fill(one)) : undefined
 
     // n×n matrix describing the QR transformations done since last convergence
-    let Qpartial = findVectors ? diag(Array(n).fill(1)) : undefined
+    let Qpartial = findVectors ? diag(Array(n).fill(one)) : undefined
 
     // last eigenvalue converged before this many steps
     let lastConvergenceBefore = 0
@@ -298,7 +288,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
       const k = 0 // TODO set close to an eigenvalue
 
       for (let i = 0; i < n; i++) {
-        arr[i][i] -= k
+        arr[i][i] = subtract(arr[i][i], k)
       }
 
       // TODO do an implicit QR transformation
@@ -306,7 +296,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
       arr = multiply(R, Q)
 
       for (let i = 0; i < n; i++) {
-        arr[i][i] += k
+        arr[i][i] = addScalar(arr[i][i], k)
       }
 
       // keep track of transformations
@@ -315,7 +305,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
       }
 
       // The rightmost diagonal element converged to an eigenvalue
-      if (n === 1 || (big ? abs(arr[n - 1][n - 2]).lessThan(prec) : abs(arr[n - 1][n - 2]) < prec)) {
+      if (n === 1 || smaller(abs(arr[n - 1][n - 2]), prec)) {
         lastConvergenceBefore = 0
         lambdas.push(arr[n - 1][n - 1])
 
@@ -326,7 +316,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
           Qtotal = multiply(Qtotal, Qpartial)
 
           if (n > 1) {
-            Qpartial = diag(Array(n - 1).fill(1))
+            Qpartial = diag(Array(n - 1).fill(one))
           }
         }
 
@@ -338,7 +328,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
         }
 
       // The rightmost diagonal 2x2 block converged
-      } else if (n === 2 || (big ? abs(arr[n - 2][n - 3]).lessThan(prec) : abs(arr[n - 2][n - 3]) < prec)) {
+      } else if (n === 2 || smaller(abs(arr[n - 2][n - 3]), prec)) {
         lastConvergenceBefore = 0
         const ll = eigenvalues2x2(
           arr[n - 2][n - 2], arr[n - 2][n - 1],
@@ -357,7 +347,7 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
           Qtotal = multiply(Qtotal, Qpartial)
 
           if (n > 2) {
-            Qpartial = diag(Array(n - 2).fill(1))
+            Qpartial = diag(Array(n - 2).fill(one))
           }
         }
 
@@ -405,7 +395,8 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
     const vectors = []
 
     for (const l of values) {
-      // ! FIXME might fail for imprecise eigenvalues, replace with an iterative eigenvector algorithm
+      // TODO replace with an iterative eigenvector algorithm
+      // (this one might fail for imprecise eigenvalues)
       const V = usolveAll(subtract(U, multiply(l, E)), b)
       for (const v of V) { vectors.push(multiply(C, v)) }
     }
@@ -420,9 +411,9 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
   function eigenvalues2x2 (a, b, c, d) {
     // λ± = ½ trA ± ½ √( tr²A - 4 detA )
     const trA = addScalar(a, d)
-    const detA = subtract(multiply(a, d), multiply(b, c))
-    const x = multiply(trA, 0.5)
-    const y = multiply(sqrt(subtract(multiply(trA, trA), multiply(4, detA))), 0.5)
+    const detA = subtract(multiplyScalar(a, d), multiplyScalar(b, c))
+    const x = multiplyScalar(trA, 0.5)
+    const y = multiplyScalar(sqrt(subtract(multiplyScalar(trA, trA), multiplyScalar(4, detA))), 0.5)
 
     return [addScalar(x, y), subtract(x, y)]
   }
@@ -436,19 +427,20 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
    */
   function jordanBase2x2 (a, b, c, d, l1, l2, prec, type) {
     const big = type === 'BigNumber'
+    const cplx = type === 'Complex'
 
-    const b0 = big ? 0 : bignumber(0)
-    const b1 = big ? 1 : bignumber(1)
+    const zero = big ?  bignumber(0) : cplx ? complex(0) : 0
+    const one = big ?  bignumber(1) : cplx ? complex(1) : 1
 
     // matrix is already upper triangular
     // return an identity matrix
-    if (big ? abs(c).lessThan(prec) : abs(c) < prec) {
-      return [[b1, b0], [b0, b1]]
+    if (smaller(abs(c), prec)) {
+      return [[one, zero], [zero, one]]
     }
 
     // matrix is diagonalizable
     // return its eigenvectors as columns
-    if (big ? abs(l1.sub(l2)).greaterThan(prec) : abs(subtract(l1, l2)) > prec) {
+    if (larger(abs(subtract(l1, l2)), prec)) {
       return [[subtract(l1, d), subtract(l2, d)], [c, c]]
     }
 
@@ -462,10 +454,10 @@ export function createComplex ({ addScalar, subtract, multiply, sqrt, abs, bignu
     const nc = subtract(c, l1)
     const nd = subtract(d, l1)
 
-    if (big ? abs(nb).lessThan(prec) : abs(nb) < prec) {
-      return [[na, b1], [nc, b0]]
+    if (smaller(abs(nb), prec)) {
+      return [[na, one], [nc, zero]]
     } else {
-      return [[nb, b0], [nd, b1]]
+      return [[nb, zero], [nd, one]]
     }
   }
 
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index e3373904fa..27e55dc70d 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -1,13 +1,13 @@
 import { clone } from '../../../utils/object'
 
-export function createRealSymmetric ({ addScalar, subtract, column, flatten, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
+export function createRealSymmetric ({ config, addScalar, subtract, column, flatten, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
   /**
    * @param {number[] | BigNumber[]} arr
    * @param {number} N
    * @param {number} prec
    * @param {'number' | 'BigNumber'} type
    */
-  function main (arr, N, prec = 1E-12, type) {
+  function main (arr, N, prec = config.epsilon, type) {
     if (type === 'number') {
       return diag(arr, prec)
     }
@@ -78,26 +78,22 @@ export function createRealSymmetric ({ addScalar, subtract, column, flatten, abs
 
   // get angle
   function getTheta (aii, ajj, aij) {
-    let th = 0
     const denom = (ajj - aii)
-    if (Math.abs(denom) <= 1E-14) {
-      th = Math.PI / 4.0
+    if (Math.abs(denom) <= config.epsilon) {
+      return Math.PI / 4.0
     } else {
-      th = 0.5 * Math.atan(2.0 * aij / (ajj - aii))
+      return 0.5 * Math.atan(2.0 * aij / (ajj - aii))
     }
-    return th
   }
 
   // get angle
   function getThetaBig (aii, ajj, aij) {
-    let th = 0
     const denom = subtract(ajj, aii)
-    if (abs(denom) <= 1E-14) {
-      th = Math.PI / 4.0
+    if (abs(denom) <= config.epsilon) {
+      return bignumber(-1).acos().div(4)
     } else {
-      th = multiplyScalar(0.5, atan(multiply(2.0, aij, inv(denom))))
+      return multiplyScalar(0.5, atan(multiply(2.0, aij, inv(denom))))
     }
-    return th
   }
 
   // update eigvec
@@ -122,8 +118,8 @@ export function createRealSymmetric ({ addScalar, subtract, column, flatten, abs
     const N = Sij.length
     const c = cos(theta)
     const s = sin(theta)
-    const Ski = createArray(N, 0)
-    const Skj = createArray(N, 0)
+    const Ski = createArray(N, bignumber(0))
+    const Skj = createArray(N, bignumber(0))
     for (let k = 0; k < N; k++) {
       Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]))
       Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]))
@@ -142,10 +138,10 @@ export function createRealSymmetric ({ addScalar, subtract, column, flatten, abs
     const s = bignumber(sin(theta))
     const c2 = multiplyScalar(c, c)
     const s2 = multiplyScalar(s, s)
-    const Aki = createArray(N, 0)
-    const Akj = createArray(N, 0)
+    const Aki = createArray(N, bignumber(0))
+    const Akj = createArray(N, bignumber(0))
     // 2cs Hij
-    const csHij = multiply(2, c, s, Hij[i][j])
+    const csHij = multiply(bignumber(2), c, s, Hij[i][j])
     //  Aii
     const Aii = addScalar(subtract(multiplyScalar(c2, Hij[i][i]), csHij), multiplyScalar(s2, Hij[j][j]))
     const Ajj = add(multiplyScalar(s2, Hij[i][i]), csHij, multiplyScalar(c2, Hij[j][j]))
@@ -157,8 +153,8 @@ export function createRealSymmetric ({ addScalar, subtract, column, flatten, abs
     // Modify Hij
     Hij[i][i] = Aii
     Hij[j][j] = Ajj
-    Hij[i][j] = 0
-    Hij[j][i] = 0
+    Hij[i][j] = bignumber(0)
+    Hij[j][i] = bignumber(0)
     // 0  to i
     for (let k = 0; k < N; k++) {
       if (k !== i && k !== j) {

From a00b7cb0741271a51385f140ece3e32f335f1892 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Mon, 8 Mar 2021 10:43:30 +0100
Subject: [PATCH 15/19] some old uncommited changes

---
 src/function/matrix/eigs.js         |  7 +--
 src/function/matrix/eigs/complex.js | 66 ++++++++++++++++++++++++-----
 2 files changed, 59 insertions(+), 14 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index c5d04a63dd..391e410df8 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -63,7 +63,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
 
   function computeValuesAndVectors (mat, prec) {
     if (prec === undefined) {
-      prec = 1e-14
+      prec = config.epsilon
     }
 
     prec = bignumber(prec)
@@ -81,8 +81,9 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
       coerceReal(arr, N)
 
       if (isSymmetric(arr, N, prec)) {
-        const type = coerceTypes(mat, arr, N)
-        return doRealSymetric(arr, N, prec, type)
+        // ! FIXME uncomment the following lines
+        //const type = coerceTypes(mat, arr, N)
+        //return doRealSymetric(arr, N, prec, type)
       }
     }
 
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index b1df625a8b..e74d5bb8ff 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -13,11 +13,6 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
       findVectors = true
     }
 
-    // TODO implement QR for complex matrices
-    if (type === 'Complex') {
-      //throw new TypeError('Complex matrices not yet supported')
-    }
-
     // TODO check if any row/col are zero except the diagonal
 
     // make sure corresponding rows and columns have similar magnitude
@@ -366,6 +361,9 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
       }
     }
 
+    // standard sorting
+    lambdas.sort((a, b) => +subtract(abs(a), abs(b)))
+
     // the algorithm didn't converge
     if (lastConvergenceBefore > 30) {
       throw Error('The eigenvalues failed to converge. Only found these eigenvalues: ' + lambdas.join(', '))
@@ -389,16 +387,44 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     const Cinv = inv(C)
     const U = multiply(Cinv, A, C)
 
+    // turn values into a kind of "multiset"
+    // this way it is easier to find eigenvectors
+    const uniqueValues = []
+    const multiplicities = []
+
+    for (const λ of values) {
+      const i = indexOf(uniqueValues, λ, equal)
+
+      if (i === -1) {
+        uniqueValues.push(λ)
+        multiplicities.push(1)
+      } else {
+        multiplicities[i] += 1
+      }
+    }
+
+    // find eigenvectors by solving U − λE = 0
+    // TODO replace with an iterative eigenvector algorithm
+    // (this one might fail for imprecise eigenvalues)
+
+    const vectors = []
+    const len = uniqueValues.length
     const b = Array(N).fill(0)
     const E = diag(Array(N).fill(1))
 
-    const vectors = []
+    for (let i = 0; i < len; i++) {
+      const λ = uniqueValues[i]
+
+      let solutions = usolveAll(subtract(U, multiply(λ, E)), b)
+      solutions = solutions.map( v => multiply(C, v) )
 
-    for (const l of values) {
-      // TODO replace with an iterative eigenvector algorithm
-      // (this one might fail for imprecise eigenvalues)
-      const V = usolveAll(subtract(U, multiply(l, E)), b)
-      for (const v of V) { vectors.push(multiply(C, v)) }
+      solutions.shift() // ignore the null vector
+
+      if (solutions.length < multiplicities[i]) {
+        //
+      }
+
+      vectors.push(...solutions)
     }
 
     return vectors
@@ -507,5 +533,23 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     return M
   }
 
+
+  /**
+   * Finds the index of an element in an array using a custom equality function
+   * @template T
+   * @param {Array<T>} arr array in which to search
+   * @param {T} el the element to find
+   * @param {function(T, T): boolean} fn the equality function, first argument is an element of `arr`, the second is always `el`
+   * @returns {number} the index of `el`, or -1 when it's not in `arr`
+   */
+  function indexOf(arr, el, fn) {
+    for (let i = 0; i < arr.length; i++) {
+      if (fn(arr[i], el)) {
+        return i
+      }
+    }
+    return -1
+  }
+
   return main
 }

From 1f2a21d9105e67eae47356ef7310ff3dbee45a2b Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Mon, 8 Mar 2021 13:40:30 +0100
Subject: [PATCH 16/19] fix small problems introduced by merging

---
 src/function/matrix/eigs.js                  |  6 +++---
 src/function/matrix/eigs/complex.js          | 14 +++++++++-----
 src/function/matrix/eigs/realSymetric.js     |  4 ++--
 test/unit-tests/function/matrix/eigs.test.js |  2 +-
 4 files changed, 15 insertions(+), 11 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index d30c619a06..af5f203da4 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -1,6 +1,8 @@
-import { clone } from '../../utils/object.js'
 import { factory } from '../../utils/factory.js'
 import { format } from '../../utils/string.js'
+import { createComplex } from './eigs/complex.js'
+import { createRealSymmetric } from './eigs/realSymetric.js'
+import { typeOf, isNumber, isBigNumber, isComplex, isFraction } from '../../utils/is.js'
 
 const name = 'eigs'
 
@@ -64,8 +66,6 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
       prec = config.epsilon
     }
 
-    prec = bignumber(prec)
-
     const size = mat.size()
 
     if (size.length !== 2 || size[0] !== size[1]) {
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index e74d5bb8ff..78b30d49a5 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,4 +1,4 @@
-import { clone } from '../../../utils/object'
+import { clone } from '../../../utils/object.js'
 
 export function createComplex ({ addScalar, subtract, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, inv, qr, usolveAll, equal, complex, larger, smaller }) {
   /**
@@ -100,7 +100,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
           // (we want to scale only by integer powers of radix,
           // so that we don't lose any precision due to round-off)
 
-          let f = big ? bigOne : 1
+          let f = one
           let c = colNorm
 
           const rowDivRadix = divideScalar(rowNorm, radix)
@@ -162,6 +162,8 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
 
     const zero = big ? bignumber(0) : cplx ? complex(0) : 0
 
+    if (big) { prec = bignumber(prec) }
+
     for (let i = 0; i < N - 2; i++) {
       // Find the largest subdiag element in the i-th col
 
@@ -243,6 +245,8 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
 
     const one = big ? bignumber(1) : cplx ? complex(1) : 1
 
+    if (big) { prec = bignumber(prec) }
+
     // The Francis Algorithm
     // The core idea of this algorithm is that doing successive
     // A' = Q⁺AQ transformations will eventually converge to block-
@@ -273,7 +277,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     // last eigenvalue converged before this many steps
     let lastConvergenceBefore = 0
 
-    while (lastConvergenceBefore <= 30) {
+    while (lastConvergenceBefore <= 100) {
       lastConvergenceBefore += 1
 
       // TODO if the convergence is slow, do something clever
@@ -365,7 +369,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     lambdas.sort((a, b) => +subtract(abs(a), abs(b)))
 
     // the algorithm didn't converge
-    if (lastConvergenceBefore > 30) {
+    if (lastConvergenceBefore > 100) {
       throw Error('The eigenvalues failed to converge. Only found these eigenvalues: ' + lambdas.join(', '))
     }
 
@@ -421,7 +425,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
       solutions.shift() // ignore the null vector
 
       if (solutions.length < multiplicities[i]) {
-        //
+        // !FIXME
       }
 
       vectors.push(...solutions)
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index 27e55dc70d..04477e10a9 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -1,4 +1,4 @@
-import { clone } from '../../../utils/object'
+import { clone } from '../../../utils/object.js'
 
 export function createRealSymmetric ({ config, addScalar, subtract, column, flatten, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
   /**
@@ -245,7 +245,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, column, flat
       let minID = 0
       let minE = E[0]
       for (let j = 0; j < E.length; j++) {
-        if (E[j] < minE) {
+        if (abs(E[j]) < abs(minE)) {
           minID = j
           minE = E[minID]
         }
diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js
index 78aae29156..9bb4439580 100644
--- a/test/unit-tests/function/matrix/eigs.test.js
+++ b/test/unit-tests/function/matrix/eigs.test.js
@@ -53,7 +53,7 @@ describe('eigs', function () {
         [-3.8571699139231796, 0.7047577966772156, 0.9122549659760404, 0.9232933211541949],
         [2.852995822026198, 0.9122549659760404, 1.6598316026960402, -1.2931270747054358],
         [4.1957619745869845, 0.9232933211541949, -1.2931270747054358, -4.665994662426116]]).values,
-    [-8.687249803623432, -0.9135495807127523, 2.26552473288741, 5.6502090685149735]
+    [-0.9135495807127523, 2.26552473288741, 5.6502090685149735, -8.687249803623432]
     )
   })
 

From ac1222ed4ca97d6b1c6cca71e15fef39772f3cfc Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Wed, 31 Mar 2021 20:41:08 +0200
Subject: [PATCH 17/19] done some polishing

---
 src/function/arithmetic/norm.js              |  2 +-
 src/function/matrix/eigs.js                  | 22 +++++----
 src/function/matrix/eigs/complex.js          | 47 +++++++++++++-------
 src/function/matrix/eigs/realSymetric.js     | 21 +++------
 test/unit-tests/function/matrix/eigs.test.js | 18 +++++---
 5 files changed, 61 insertions(+), 49 deletions(-)

diff --git a/src/function/arithmetic/norm.js b/src/function/arithmetic/norm.js
index 9ef2632b1f..7a86097d58 100644
--- a/src/function/arithmetic/norm.js
+++ b/src/function/arithmetic/norm.js
@@ -240,7 +240,7 @@ export const createNorm = /* #__PURE__ */ factory(
       const tx = ctranspose(x)
       const squaredX = multiply(tx, x)
       const eigenVals = eigs(squaredX).values
-      const rho = eigenVals.get([eigenVals.size()[0] - 1])
+      const rho = eigenVals[eigenVals.length - 1]
       return abs(sqrt(rho))
     }
 
diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index af5f203da4..e3a5c7a624 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -6,12 +6,11 @@ import { typeOf, isNumber, isBigNumber, isComplex, isFraction } from '../../util
 
 const name = 'eigs'
 
-const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolveAll', 'im', 're', 'smaller']
-
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, qr, usolveAll, im, re, smaller }) => {
-
+// The absolute state of math.js's dependency system:
+const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolveAll', 'im', 're', 'smaller', 'round', 'log10', 'transpose']
+export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, qr, usolveAll, im, re, smaller, round, log10, transpose }) => {
   const doRealSymetric = createRealSymmetric({ config, addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
-  const doComplex = createComplex({ config, addScalar, subtract, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, qr, inv, usolveAll, equal, complex, larger, smaller })
+  const doComplex = createComplex({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, qr, inv, usolveAll, equal, complex, larger, smaller, round, log10, transpose })
 
   /**
    * Compute eigenvalues and eigenvectors of a matrix.
@@ -22,12 +21,13 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
    *
    * Examples:
    *
+   *     const { eigs, multiply, column, transpose } = math
    *     const H = [[5, 2.3], [2.3, 1]]
-   *     const ans = math.eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}
+   *     const ans = eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}
    *     const E = ans.values
    *     const U = ans.vectors
-   *     math.multiply(H, math.column(U, 0)) // returns math.multiply(E[0], math.column(U, 0))
-   *     const UTxHxU = math.multiply(math.transpose(U), H, U) // rotates H to the eigen-representation
+   *     multiply(H, column(U, 0)) // returns multiply(E[0], column(U, 0))
+   *     const UTxHxU = multiply(transpose(U), H, U) // rotates H to the eigen-representation
    *     E[0] == UTxHxU[0][0]  // returns true
    *
    * See also:
@@ -79,9 +79,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
       coerceReal(arr, N)
 
       if (isSymmetric(arr, N, prec)) {
-        // ! FIXME uncomment the following lines
-        //const type = coerceTypes(mat, arr, N)
-        //return doRealSymetric(arr, N, prec, type)
+        const type = coerceTypes(mat, arr, N)
+        return doRealSymetric(arr, N, prec, type)
       }
     }
 
@@ -103,7 +102,6 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
     return true
   }
 
-
   /** @return {boolean} */
   function isReal (arr, N, prec) {
     for (let i = 0; i < N; i++) {
diff --git a/src/function/matrix/eigs/complex.js b/src/function/matrix/eigs/complex.js
index 78b30d49a5..c1eb9c5c02 100644
--- a/src/function/matrix/eigs/complex.js
+++ b/src/function/matrix/eigs/complex.js
@@ -1,6 +1,6 @@
 import { clone } from '../../../utils/object.js'
 
-export function createComplex ({ addScalar, subtract, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, inv, qr, usolveAll, equal, complex, larger, smaller }) {
+export function createComplex ({ addScalar, subtract, flatten, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, inv, qr, usolveAll, equal, complex, larger, smaller, round, log10, transpose }) {
   /**
    * @param {number[][]} arr the matrix to find eigenvalues of
    * @param {number} N size of the matrix
@@ -44,7 +44,8 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     let vectors
 
     if (findVectors) {
-      vectors = findEigenvectors(arr, N, C, values)
+      vectors = findEigenvectors(arr, N, C, values, prec)
+      vectors = transpose((vectors)) // vectors are columns of a matrix
     }
 
     return { values, vectors }
@@ -62,7 +63,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     const cplx = type === 'Complex'
 
     const zero = big ? bignumber(0) : cplx ? complex(0) : 0
-    const one  = big ? bignumber(1) : cplx ? complex(1) : 1
+    const one = big ? bignumber(1) : cplx ? complex(1) : 1
 
     // base of the floating-point arithmetic
     const radix = big ? bignumber(10) : 2
@@ -90,7 +91,6 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
         for (let j = 0; j < N; j++) {
           if (i === j) continue
           const c = abs(arr[i][j])
-          const r = abs(arr[j][i])
           colNorm = addScalar(colNorm, c)
           rowNorm = addScalar(rowNorm, c)
         }
@@ -117,7 +117,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
 
           // check whether balancing is needed
           // condition = (c + rowNorm) / f < 0.95 * (colNorm + rowNorm)
-          const condition = smaller( divideScalar(addScalar(c, rowNorm), f), multiplyScalar(addScalar(colNorm, rowNorm), 0.95))
+          const condition = smaller(divideScalar(addScalar(c, rowNorm), f), multiplyScalar(addScalar(colNorm, rowNorm), 0.95))
 
           // apply balancing similarity transformation
           if (condition) {
@@ -370,7 +370,10 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
 
     // the algorithm didn't converge
     if (lastConvergenceBefore > 100) {
-      throw Error('The eigenvalues failed to converge. Only found these eigenvalues: ' + lambdas.join(', '))
+      const err = Error('The eigenvalues failed to converge. Only found these eigenvalues: ' + lambdas.join(', '))
+      err.values = lambdas
+      err.vectors = []
+      throw err
     }
 
     // combine the overall QR transformation Qtotal with the subsequent
@@ -387,7 +390,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
    * @param {number[]} values array of eigenvalues of A
    * @returns {Matrix[]} eigenvalues
    */
-  function findEigenvectors (A, N, C, values) {
+  function findEigenvectors (A, N, C, values, prec) {
     const Cinv = inv(C)
     const U = multiply(Cinv, A, C)
 
@@ -399,8 +402,12 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     for (const λ of values) {
       const i = indexOf(uniqueValues, λ, equal)
 
+      // a dirty trick that helps us find more vectors
+      // TODO with iterative algorithm this can be removed
+      const roundedλ = round(λ, subtract(-1, log10(prec)))
+
       if (i === -1) {
-        uniqueValues.push(λ)
+        uniqueValues.push(roundedλ)
         multiplicities.push(1)
       } else {
         multiplicities[i] += 1
@@ -416,19 +423,30 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     const b = Array(N).fill(0)
     const E = diag(Array(N).fill(1))
 
+    // eigenvalues for which usolve failed (due to numerical error)
+    const failedLambdas = []
+
     for (let i = 0; i < len; i++) {
       const λ = uniqueValues[i]
 
       let solutions = usolveAll(subtract(U, multiply(λ, E)), b)
-      solutions = solutions.map( v => multiply(C, v) )
+      solutions = solutions.map(v => multiply(C, v))
 
       solutions.shift() // ignore the null vector
 
+      // looks like we missed something
       if (solutions.length < multiplicities[i]) {
-        // !FIXME
+        failedLambdas.push(λ)
       }
 
-      vectors.push(...solutions)
+      vectors.push(...solutions.map(v => flatten(v)))
+    }
+
+    if (failedLambdas.length !== 0) {
+      const err = new Error('Failed to find eigenvectors for the following eigenvalues: ' + failedLambdas.join(', '))
+      err.values = values
+      err.vectors = vectors
+      throw err
     }
 
     return vectors
@@ -459,8 +477,8 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     const big = type === 'BigNumber'
     const cplx = type === 'Complex'
 
-    const zero = big ?  bignumber(0) : cplx ? complex(0) : 0
-    const one = big ?  bignumber(1) : cplx ? complex(1) : 1
+    const zero = big ? bignumber(0) : cplx ? complex(0) : 0
+    const one = big ? bignumber(1) : cplx ? complex(1) : 1
 
     // matrix is already upper triangular
     // return an identity matrix
@@ -537,7 +555,6 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
     return M
   }
 
-
   /**
    * Finds the index of an element in an array using a custom equality function
    * @template T
@@ -546,7 +563,7 @@ export function createComplex ({ addScalar, subtract, multiply, multiplyScalar,
    * @param {function(T, T): boolean} fn the equality function, first argument is an element of `arr`, the second is always `el`
    * @returns {number} the index of `el`, or -1 when it's not in `arr`
    */
-  function indexOf(arr, el, fn) {
+  function indexOf (arr, el, fn) {
     for (let i = 0; i < arr.length; i++) {
       if (fn(arr[i], el)) {
         return i
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymetric.js
index 04477e10a9..57c215f2c6 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymetric.js
@@ -1,6 +1,6 @@
 import { clone } from '../../../utils/object.js'
 
-export function createRealSymmetric ({ config, addScalar, subtract, column, flatten, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
+export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
   /**
    * @param {number[] | BigNumber[]} arr
    * @param {number} N
@@ -236,10 +236,11 @@ export function createRealSymmetric ({ config, addScalar, subtract, column, flat
   // sort results
   function sorting (E, S) {
     const N = E.length
-    const Ef = Array(N)
-    const Sf = Array(N)
+    const values = Array(N)
+    const vectors = Array(N)
+
     for (let k = 0; k < N; k++) {
-      Sf[k] = Array(N)
+      vectors[k] = Array(N)
     }
     for (let i = 0; i < N; i++) {
       let minID = 0
@@ -250,21 +251,13 @@ export function createRealSymmetric ({ config, addScalar, subtract, column, flat
           minE = E[minID]
         }
       }
-      Ef[i] = E.splice(minID, 1)[0]
+      values[i] = E.splice(minID, 1)[0]
       for (let k = 0; k < N; k++) {
-        Sf[k][i] = S[k][minID]
+        vectors[k][i] = S[k][minID]
         S[k].splice(minID, 1)
       }
     }
 
-    const values = clone(Ef)
-    const vectors = []
-
-    for (let i = 0; i < Sf.length; i++) {
-      vectors[i] = flatten(column(Sf, i))
-    }
-
-    // !FIXME vectors are always Array[], never Matrix[]
     return { values, vectors }
   }
 
diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js
index 9bb4439580..84846100e6 100644
--- a/test/unit-tests/function/matrix/eigs.test.js
+++ b/test/unit-tests/function/matrix/eigs.test.js
@@ -2,9 +2,10 @@ import assert from 'assert'
 import math from '../../../../src/defaultInstance.js'
 import approx from '../../../../tools/approx.js'
 const eigs = math.eigs
+const complex = math.complex
 
 describe('eigs', function () {
-  it('should only accept a square matrix', function () {
+  it('only accepts a square matrix', function () {
     assert.throws(function () { eigs(math.matrix([[1, 2, 3], [4, 5, 6]])) }, /Matrix must be square/)
     assert.throws(function () { eigs([[1, 2, 3], [4, 5, 6]]) }, /Matrix must be square/)
     assert.throws(function () { eigs([[1, 2], [4, 5, 6]]) }, /DimensionError: Dimension mismatch/)
@@ -13,7 +14,7 @@ describe('eigs', function () {
     assert.throws(function () { eigs('random') }, /TypeError: Unexpected type of argument/)
   })
 
-  it('should only accept a matrix with valid element type', function () {
+  it('only accepts a matrix with valid element type', function () {
     assert.throws(function () { eigs([['x', 2], [4, 5]]) }, /Cannot convert "x" to a number/)
   })
 
@@ -25,9 +26,12 @@ describe('eigs', function () {
     approx.deepEqual(eigs(
       [[2, 0, 0], [0, 1, 0], [0, 0, 5]]).values, [1, 2, 5]
     )
+    approx.deepEqual(eigs(
+      [[complex(2, 1), 0, 0], [0, 1, 0], [0, 0, complex(0, 5)]]).values, [complex(1, 0), complex(2, 1), complex(0, 5)]
+    )
   })
 
-  it('eigenvalue check for 2x2 simple matrix', function () {
+  it('calculates eigenvalues for 2x2 simple matrix', function () {
     // 2x2 test
     approx.deepEqual(eigs(
       [[1, 0.1], [0.1, 1]]).values, [0.9, 1.1]
@@ -40,7 +44,7 @@ describe('eigs', function () {
     )
   })
 
-  it('eigenvalue check for 3x3 and 4x4 matrix', function () {
+  it('calculates eigenvalues for 3x3 and 4x4 matrix', function () {
     // 3x3 test and 4x4
     approx.deepEqual(eigs(
       [[1.0, 1.0, 1.0],
@@ -75,7 +79,7 @@ describe('eigs', function () {
     approx.deepEqual(Ei, E)
   })
 
-  it('fractions are supported', function () {
+  it('supports fractions', function () {
     const aij = math.fraction('1/2')
     approx.deepEqual(eigs(
       [[aij, aij, aij],
@@ -85,7 +89,7 @@ describe('eigs', function () {
     )
   })
 
-  it('bigNumber diagonalization is supported', function () {
+  it('diagonalizes matrix with bigNumber', function () {
     const x = [[math.bignumber(1), math.bignumber(0)], [math.bignumber(0), math.bignumber(1)]]
     approx.deepEqual(eigs(x).values, [math.bignumber(1), math.bignumber(1)])
     const y = [[math.bignumber(1), math.bignumber(1.0)], [math.bignumber(1.0), math.bignumber(1)]]
@@ -109,7 +113,7 @@ describe('eigs', function () {
     approx.deepEqual(Ei, E)
   })
 
-  it('make sure BigNumbers input is actually calculated with BigNumber precision', function () {
+  it('actually calculates BigNumbers input with BigNumber precision', function () {
     const B = math.bignumber([
       [0, 1],
       [1, 0]

From ab58046b2035c38d38f4ed473e9e89d3d6107fc9 Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Wed, 31 Mar 2021 22:44:28 +0200
Subject: [PATCH 18/19] improved jsdoc description of eigs

---
 src/function/matrix/eigs.js | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index e3a5c7a624..e73b22e4b8 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -13,7 +13,11 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
   const doComplex = createComplex({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, qr, inv, usolveAll, equal, complex, larger, smaller, round, log10, transpose })
 
   /**
-   * Compute eigenvalues and eigenvectors of a matrix.
+   * Compute eigenvalues and eigenvectors of a matrix. The eigenvalues are sorted by their absolute value, ascending.
+   * Each eigenvalue can be listed several times, according to its multiplicity. The eigenvectors are returned as columns
+   * of a matrix – the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`).
+   * If the algorithm fails to converge, it will throw an error – in that case, however, you may still find useful information
+   * in `err.values` and `err.vectors`.
    *
    * Syntax:
    *

From fd697459241f0e887e7970380704847a14bf5f3e Mon Sep 17 00:00:00 2001
From: Michal Grno <m93a.cz@gmail.com>
Date: Wed, 31 Mar 2021 22:57:39 +0200
Subject: [PATCH 19/19] little changes in jsdoc

---
 src/function/matrix/eigs.js | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index e73b22e4b8..f061ee813a 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -14,8 +14,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
 
   /**
    * Compute eigenvalues and eigenvectors of a matrix. The eigenvalues are sorted by their absolute value, ascending.
-   * Each eigenvalue can be listed several times, according to its multiplicity. The eigenvectors are returned as columns
-   * of a matrix – the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`).
+   * An eigenvalue with multiplicity k will be listed k times. The eigenvectors are returned as columns of a matrix –
+   * the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`).
    * If the algorithm fails to converge, it will throw an error – in that case, however, you may still find useful information
    * in `err.values` and `err.vectors`.
    *
@@ -31,7 +31,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
    *     const E = ans.values
    *     const U = ans.vectors
    *     multiply(H, column(U, 0)) // returns multiply(E[0], column(U, 0))
-   *     const UTxHxU = multiply(transpose(U), H, U) // rotates H to the eigen-representation
+   *     const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H
    *     E[0] == UTxHxU[0][0]  // returns true
    *
    * See also:
@@ -41,7 +41,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
    * @param {Array | Matrix} x  Matrix to be diagonalized
    *
    * @param {number | BigNumber} [prec] Precision, default value: 1e-15
-   * @return {{values: Array, vectors: Matrix[]}} Object containing an array of eigenvalues and an array of eigenvectors.
+   * @return {{values: Array, vectors: Matrix}} Object containing an array of eigenvalues and a matrix with eigenvectors as columns.
    *
    */
   return typed('eigs', {