From 5624e6938f59cbcfb6f1cdb7873b64ce6d31ad99 Mon Sep 17 00:00:00 2001
From: jos <wjosdejong@gmail.com>
Date: Sun, 29 Mar 2020 16:17:19 +0200
Subject: [PATCH] Fix #1789: Function `eigs` not calculating with BigNumber
 precision when input contains BigNumbers

---
 HISTORY.md                                   |  2 ++
 src/function/matrix/eigs.js                  | 34 +++++++++-----------
 test/unit-tests/function/matrix/eigs.test.js | 13 ++++++++
 3 files changed, 30 insertions(+), 19 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index ac07ec355d..82fea999ed 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -3,6 +3,8 @@
 
 # not yet published, version 6.6.2
 
+- Fix #1789: Function `eigs` not calculating with BigNumber precision 
+  when input contains BigNumbers.
 - Run the build script during npm `prepare`, so you can use the library
   directly when installing directly from git. Thanks @cinderblock.
 
diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index 304da41e18..a4a7274561 100644
--- a/src/function/matrix/eigs.js
+++ b/src/function/matrix/eigs.js
@@ -3,9 +3,9 @@ import { factory } from '../../utils/factory'
 import { format } from '../../utils/string'
 
 const name = 'eigs'
-const dependencies = ['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']
 
-export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ 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 }) => {
   /**
    * Compute eigenvalue and eigenvector of a real symmetric matrix.
    * Only applicable to two dimensional symmetric matrices. Uses Jacobi
@@ -172,26 +172,22 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
 
   // 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
     } else {
-      th = 0.5 * Math.atan(2.0 * aij / (ajj - aii))
+      return 0.5 * Math.atan(2 * 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, aij, inv(denom))))
     }
-    return th
   }
 
   // update eigvec
@@ -216,8 +212,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     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]))
@@ -236,10 +232,10 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     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]))
@@ -251,8 +247,8 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     // 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) {
diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js
index 8a0dddc20e..483628753e 100644
--- a/test/unit-tests/function/matrix/eigs.test.js
+++ b/test/unit-tests/function/matrix/eigs.test.js
@@ -105,4 +105,17 @@ describe('eigs', function () {
     }
     approx.deepEqual(Ei, E)
   })
+
+  it('make sure BigNumbers input is actually calculated with BigNumber precision', function () {
+    const B = math.bignumber([
+      [0, 1],
+      [1, 0]
+    ])
+    const eig = math.eigs(B)
+
+    assert.strictEqual(eig.values[0].toString(),
+      '-0.9999999999999999999999999999999999999999999999999999999999999999')
+    assert.strictEqual(eig.values[1].toString(),
+      '0.9999999999999999999999999999999999999999999999999999999999999999')
+  })
 })