feat: add BisectMethod impl

Author: tiankaima

CustomMath_lib/BisectMethod/BisectMethod.cpp

@@ -0,0 +1,104 @@
+//
+// Created by TianKai Ma on 2023/12/19.
+//
+
+#include "BisectMethod.h"
+
+ull CalculateSignChange(const Vector &x, const Vector &y, lld mu) {
+#ifdef DEBUG
+    if (x.size != y.size + 1) {
+        throw std::invalid_argument("x.size != y.size + 1");
+    }
+#endif
+
+    auto n = x.size;
+    auto s = 0;
+    auto q = x.array[0] - mu;
+    for (ull k = 0; k < n; k++) {
+        if (q < 0) {
+            s++;
+        }
+        if (k != n - 1) {
+            if (q == 0) {
+                q = std::abs(y.array[k]) * 1e-10;
+            }
+
+            q = x.array[k + 1] - mu - y.array[k] * y.array[k] / q;
+        }
+    }
+    return s;
+}
+
+std::vector<lld> BisectMethodCall(const Vector &x, const Vector &y, lld start, lld end, lld precision) {
+#ifdef DEBUG
+    if (start > end) {
+        throw std::invalid_argument("start > end");
+    }
+    if (precision < 0) {
+        throw std::invalid_argument("precision < 0");
+    }
+#endif
+
+    if (end - start < precision) {
+        auto l = CalculateSignChange(x, y, start);
+        auto r = CalculateSignChange(x, y, end);
+
+        if (l > r) {
+            return {};
+        }
+
+        auto result = std::vector<lld>(r - l, (start + end) / 2);
+        return result;
+    }
+
+    auto mid = (start + end) / 2;
+    auto l = CalculateSignChange(x, y, start);
+    auto r = CalculateSignChange(x, y, end);
+
+    if (l == r) {
+        return {};
+    }
+
+    auto m = CalculateSignChange(x, y, mid);
+    auto result = std::vector<lld>();
+    if (l != m) {
+        auto left = BisectMethodCall(x, y, start, mid, precision);
+        result.insert(result.end(), left.begin(), left.end());
+    }
+    if (m != r) {
+        auto right = BisectMethodCall(x, y, mid, end, precision);
+        result.insert(result.end(), right.begin(), right.end());
+    }
+    return result;
+}
+
+Vector BisectMethod(const Vector &x, const Vector &y, lld precision) {
+#ifdef DEBUG
+    if (x.size != y.size + 1) {
+        throw std::invalid_argument("x.size != y.size + 1");
+    }
+    if (precision < 0) {
+        throw std::invalid_argument("precision < 0");
+    }
+#endif
+
+    // determine Infinity norm of T: x_n + y_n + y_(n-1)
+    auto n = x.size;
+    auto max = x.array[0];
+    for (ull i = 1; i < n; i++) {
+        lld temp;
+        if (i == n - 1) {
+            temp = std::abs(x.array[i]) + std::abs(y.array[i]);
+        } else {
+            temp = std::abs(x.array[i]) + std::abs(y.array[i - 1]) + std::abs(y.array[i]);
+        }
+        if (temp > max) {
+            max = temp;
+        }
+    }
+
+    auto l = -max;
+    auto r = max;
+    auto result = BisectMethodCall(x, y, l, r, precision);
+    return Vector(result);
+}

CustomMath_lib/BisectMethod/BisectMethod.h

@@ -0,0 +1,14 @@
+//
+// Created by TianKai Ma on 2023/12/19.
+//
+
+#ifndef NUMERICAL_ALGEBRA_BISECTMETHOD_H
+#define NUMERICAL_ALGEBRA_BISECTMETHOD_H
+
+#include "CustomMath_lib.h"
+
+ull CalculateSignChange(const Vector &x, const Vector &y, lld mu);
+
+Vector BisectMethod(const Vector &x, const Vector &y, lld precision = 1e-10);
+
+#endif //NUMERICAL_ALGEBRA_BISECTMETHOD_H

CustomMath_lib/CMakeLists.txt

@@ -43,6 +43,8 @@ set(SOURCE_FILES
         Matrix/MatrixOperations.cpp
         JacobiMethod/JacobiMethod.cpp
         JacobiMethod/JacobiMethod.h
+        BisectMethod/BisectMethod.cpp
+        BisectMethod/BisectMethod.h
 )
 
 add_library(CustomMath_lib STATIC ${SOURCE_FILES} ${HEADER_FILES})

CustomMath_lib/CustomMath_lib.h

@@ -10,6 +10,7 @@
 #include "Vector/Vector.h"
 #include "Matrix/Matrix.h"
 
+#include "BisectMethod/BisectMethod.h"
 #include "CholeskyMethod/CholeskyMethod.h"
 #include "GaussMethod/GaussMethod.h"
 #include "HessenbergMethod/HessenbergMethod.h"

main.cpp

@@ -1,16 +1,38 @@
 #include "CustomMath_lib.h"
 
 int main() {
-    auto m = Matrix("[5 1 -2; 1 2 0; -2 0 -10]");
-    auto I = Matrix::identity(m.rows);
-    auto lambda = -10.263471;
-    auto A = m - I * lambda;
-    auto x = Vector("[1 1 1]");
+    auto n = 3;
+    auto x = Vector(n, 2);
+    auto y = Vector(n - 1, -1);
 
-    auto result = RevPowerIteration(PowerIterationInput{A, x, 1000});
+    auto A = Matrix(n, n);
+    for(ull i = 0; i < n; i++) {
+        for(ull j = 0; j < n; j++) {
+            if(i == j) {
+                A.matrix[i][j] = 2;
+            } else if(i == j + 1 || i == j - 1) {
+                A.matrix[i][j] = -1;
+            }
+        }
+    }
+    auto I = Matrix::identity(n);
+    auto x_default = Vector(n, 1);
 
-    std::cout << "Eigenvector: " << std::endl;
-    result.result.print();
+//    for (lld k = -4.0; k <= 4.0; k += 0.1) {
+//        auto r = CalculateSignChange(x, y, k);
+//        std::cout << k << " " << r << std::endl;
+//    }
+
+    auto r = BisectMethod(x, y);
+    r.print();
+
+    for (auto &lambda: r.array) {
+        auto B = A - I * lambda;
+        auto input = PowerIterationInput{B, x_default, 1000,};
+        auto k = RevPowerIteration(input);
+        k.result.print();
+        (B * k.result).print();
+    }
 
     return 0;
 }