From 1116ab2babd2ff5fecfe449d4b50d8f932476707 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 07:41:08 -0400
Subject: [PATCH 01/18] rename computer.... folder to numerical methods

---
 CMakeLists.txt                                                  | 2 +-
 .../CMakeLists.txt                                              | 2 +-
 .../bisection_method.cpp                                        | 0
 .../false_position.cpp                                          | 0
 .../gaussian_elimination.cpp                                    | 0
 .../newton_raphson_method.cpp                                   | 0
 .../ordinary_least_squares_regressor.cpp                        | 0
 .../successive_approximation.cpp                                | 0
 8 files changed, 2 insertions(+), 2 deletions(-)
 rename {computer_oriented_statistical_methods => numerical_methods}/CMakeLists.txt (91%)
 rename {computer_oriented_statistical_methods => numerical_methods}/bisection_method.cpp (100%)
 rename {computer_oriented_statistical_methods => numerical_methods}/false_position.cpp (100%)
 rename {computer_oriented_statistical_methods => numerical_methods}/gaussian_elimination.cpp (100%)
 rename {computer_oriented_statistical_methods => numerical_methods}/newton_raphson_method.cpp (100%)
 rename {computer_oriented_statistical_methods => numerical_methods}/ordinary_least_squares_regressor.cpp (100%)
 rename {computer_oriented_statistical_methods => numerical_methods}/successive_approximation.cpp (100%)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index d304c776f..430be40d6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -70,7 +70,7 @@ add_subdirectory(sorting)
 add_subdirectory(geometry)
 add_subdirectory(probability)
 add_subdirectory(machine_learning)
-add_subdirectory(computer_oriented_statistical_methods)
+add_subdirectory(numerical_methods)
 
 set(CPACK_PROJECT_NAME ${PROJECT_NAME})
 set(CPACK_PROJECT_VERSION ${PROJECT_VERSION})
diff --git a/computer_oriented_statistical_methods/CMakeLists.txt b/numerical_methods/CMakeLists.txt
similarity index 91%
rename from computer_oriented_statistical_methods/CMakeLists.txt
rename to numerical_methods/CMakeLists.txt
index acff1ed1c..fc4c343de 100644
--- a/computer_oriented_statistical_methods/CMakeLists.txt
+++ b/numerical_methods/CMakeLists.txt
@@ -13,6 +13,6 @@ foreach( testsourcefile ${APP_SOURCES} )
     if(OpenMP_CXX_FOUND)
         target_link_libraries(${testname} OpenMP::OpenMP_CXX)
     endif()
-    install(TARGETS ${testname} DESTINATION "bin/stats")
+    install(TARGETS ${testname} DESTINATION "bin/numerical_methods")
 
 endforeach( testsourcefile ${APP_SOURCES} )
diff --git a/computer_oriented_statistical_methods/bisection_method.cpp b/numerical_methods/bisection_method.cpp
similarity index 100%
rename from computer_oriented_statistical_methods/bisection_method.cpp
rename to numerical_methods/bisection_method.cpp
diff --git a/computer_oriented_statistical_methods/false_position.cpp b/numerical_methods/false_position.cpp
similarity index 100%
rename from computer_oriented_statistical_methods/false_position.cpp
rename to numerical_methods/false_position.cpp
diff --git a/computer_oriented_statistical_methods/gaussian_elimination.cpp b/numerical_methods/gaussian_elimination.cpp
similarity index 100%
rename from computer_oriented_statistical_methods/gaussian_elimination.cpp
rename to numerical_methods/gaussian_elimination.cpp
diff --git a/computer_oriented_statistical_methods/newton_raphson_method.cpp b/numerical_methods/newton_raphson_method.cpp
similarity index 100%
rename from computer_oriented_statistical_methods/newton_raphson_method.cpp
rename to numerical_methods/newton_raphson_method.cpp
diff --git a/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp b/numerical_methods/ordinary_least_squares_regressor.cpp
similarity index 100%
rename from computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp
rename to numerical_methods/ordinary_least_squares_regressor.cpp
diff --git a/computer_oriented_statistical_methods/successive_approximation.cpp b/numerical_methods/successive_approximation.cpp
similarity index 100%
rename from computer_oriented_statistical_methods/successive_approximation.cpp
rename to numerical_methods/successive_approximation.cpp

From eb977ba9d2220d94344d75e3ece98f8e746c35e2 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 11:32:57 -0400
Subject: [PATCH 02/18] added durand kerner method for root computation for
 arbitrarily large polynomials

---
 numerical_methods/durand_kerner_roots.cpp | 226 ++++++++++++++++++++++
 1 file changed, 226 insertions(+)
 create mode 100644 numerical_methods/durand_kerner_roots.cpp

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
new file mode 100644
index 000000000..42ed474e0
--- /dev/null
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -0,0 +1,226 @@
+/**
+ * @file
+ * Compute riots of any given polynomial.
+ *
+ * Test the algorithm online:
+ * https://gist.github.com/kvedala/27f1b0b6502af935f6917673ec43bcd7
+ *
+ * Try the highly unstable Wilkinson's polynomial:
+ * ```
+ * ./numerical_methods/durand_kerner_roots 1 -210 20615 -1256850 53327946
+ * -1672280820 40171771630 -756111184500 11310276995381 -135585182899530
+ * 1307535010540395 -10142299865511450 63030812099294896 -311333643161390640
+ * 1206647803780373360 -3599979517947607200 8037811822645051776
+ * -12870931245150988800 13803759753640704000 -8752948036761600000
+ * 2432902008176640000
+ * ```
+ */
+
+#include <algorithm>
+#include <climits>
+#include <cmath>
+#include <complex>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <valarray>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#define ACCURACY 1e-10 /**< maximum accuracy limit */
+
+/**
+ * Evaluate the value of a polynomial with given coefficients
+ * \param[in] coeffs coefficients of the polynomial
+ * \param[in] degree degree of polynomial
+ * \param[in] x point at which to evaluate the polynomial
+ * \returns \f$f(x)\f$
+ **/
+std::complex<double> poly_function(const std::valarray<double> &coeffs,
+                                   std::complex<double> x) {
+    double real = 0.f, imag = 0.f;
+    int n;
+
+#ifdef _OPENMP
+#pragma omp parallel for reduction(+ : real, imag)
+#endif
+    for (n = 0; n < coeffs.size(); n++) {
+        std::complex<double> tmp =
+            coeffs[n] * std::pow(x, coeffs.size() - n - 1);
+        real += tmp.real();
+        imag += tmp.imag();
+    }
+
+    return std::complex<double>(real, imag);
+}
+
+/**
+ * create a textual form of complex number
+ * \param[in] x point at which to evaluate the polynomial
+ * \returns pointer to converted string
+ */
+const char *complex_str(const std::complex<double> &x) {
+    static char msg[50];
+
+    sprintf(msg, "% 7.04g%+7.04gj", x.real(), x.imag());
+
+    return msg;
+}
+
+/**
+ * check for termination condition
+ * \param[in] delta point at which to evaluate the polynomial
+ * \returns `false` if termination not reached
+ * \returns `true` if termination reached
+ */
+bool check_termination(long double delta) {
+    static long double past_delta = INFINITY;
+    if (std::abs(past_delta - delta) <= ACCURACY || delta < ACCURACY)
+        return true;
+    past_delta = delta;
+    return false;
+}
+
+/***
+ * Main function.
+ * The comandline input arguments are taken as coeffiecients of a polynomial.
+ * For example, this command
+ * ```sh
+ * ./durand_kerner_roots 1 0 -4
+ * ```
+ * will find roots of the polynomial \f$1\cdot x^2 + 0\cdot x^1 + (-4)=0\f$
+ **/
+int main(int argc, char **argv) {
+    unsigned int degree = 0;
+    unsigned int n, i;
+
+    if (argc < 2) {
+        std::cout
+            << "Please pass the coefficients of the polynomial as commandline "
+               "arguments.\n";
+        return 0;
+    }
+
+    degree = argc - 1;                     // detected polynomial degree
+    std::valarray<double> coeffs(degree);  // create coefficiencts array
+
+    /* initialize random seed: */
+    std::srand(std::time(nullptr));
+
+#if defined(DEBUG) || !defined(NDEBUG)
+    /*
+     * store intermediate values to a CSV file
+     */
+    std::ofstream log_file;
+    log_file.open("durand_kerner.log.csv");
+    if (!log_file.is_open()) {
+        perror("Unable to create a storage log file!");
+        return EXIT_FAILURE;
+    }
+    log_file << "iter#,";
+#endif
+
+    // number of roots = degree - 1
+    std::valarray<std::complex<double>> s0(degree - 1);
+
+    std::cout << "Computing the roots for:\n\t";
+    for (n = 0; n < degree; n++) {
+        coeffs[n] = strtod(argv[n + 1], nullptr);
+        if (n < degree - 1 && coeffs[n] != 0)
+            std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1 << " + ,";
+        else if (coeffs[n] != 0)
+            std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1
+                      << " = 0\n";
+
+        /* initialize root approximations with random values */
+        if (n < degree - 1) {
+            s0[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
+            s0[n] -= 50.f;
+            s0[n] /= 50.f;
+#if defined(DEBUG) || !defined(NDEBUG)
+            log_file << "root_" << n << ",";
+#endif
+        }
+    }
+
+    // numerical errors less when the first coefficient is "1"
+    // hence, we normalize the first coefficient
+    coeffs /= coeffs[0];
+
+#if defined(DEBUG) || !defined(NDEBUG)
+    log_file << "avg. correction";
+    log_file << "\n0,";
+    for (n = 0; n < degree - 1; n++) log_file << complex_str(s0[n]) << ",";
+#endif
+
+    double tol_condition = 1;
+    unsigned long iter = 0;
+
+    clock_t end_time, start_time = clock();
+    while (!check_termination(tol_condition) && iter < INT_MAX) {
+        std::complex<double> delta = 0;
+        tol_condition = 0;
+        iter++;
+
+#if defined(DEBUG) || !defined(NDEBUG)
+        log_file << "\n" << iter << ",";
+#endif
+
+        for (n = 0; n < degree - 1; n++) {
+            std::complex<double> numerator, denominator;
+            numerator = poly_function(coeffs, s0[n]);
+            denominator = 1.0;
+            for (i = 0; i < degree - 1; i++)
+                if (i != n)
+                    denominator *= s0[n] - s0[i];
+
+            delta = numerator / denominator;
+
+            if (std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
+                std::cerr << "\n\nOverflow/underrun error - got value = "
+                          << std::abs(delta);
+                goto end;
+            }
+
+            s0[n] -= delta;
+
+            tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
+
+#if defined(DEBUG) || !defined(NDEBUG)
+            log_file << complex_str(s0[n]) << ",";
+#endif
+        }
+        // tol_condition /= (degree - 1);
+
+#if defined(DEBUG) || !defined(NDEBUG)
+        if (iter % 500 == 0) {
+            std::cout << "Iter: " << iter << "\t";
+            for (n = 0; n < degree - 1; n++)
+                std::cout << "\t" << complex_str(s0[n]);
+            std::cout << "\t\tabsolute average change: " << tol_condition
+                      << "\n";
+        }
+
+        log_file << tol_condition;
+#endif
+    }
+end:
+
+    end_time = clock();
+
+#if defined(DEBUG) || !defined(NDEBUG)
+    log_file.close();
+#endif
+
+    std::cout << "\nIterations: " << iter << "\n";
+    for (n = 0; n < degree - 1; n++)
+        std::cout << "\t" << complex_str(s0[n]) << "\n";
+    std::cout << "absolute average change: " << tol_condition << "\n";
+    std::cout << "Time taken: "
+              << (end_time - start_time) / (double)CLOCKS_PER_SEC << " sec\n";
+
+    return 0;
+}
\ No newline at end of file

From bbd578d4a15c818b35c5d901433cc72c7baaa2a1 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 11:38:34 -0400
Subject: [PATCH 03/18] fixed additional comma

---
 numerical_methods/durand_kerner_roots.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 42ed474e0..dd746ad7c 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -130,7 +130,7 @@ int main(int argc, char **argv) {
     for (n = 0; n < degree; n++) {
         coeffs[n] = strtod(argv[n + 1], nullptr);
         if (n < degree - 1 && coeffs[n] != 0)
-            std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1 << " + ,";
+            std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1 << " + ";
         else if (coeffs[n] != 0)
             std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1
                       << " = 0\n";

From 355d84dba91c17d453f54e4c78ae65c8ac2c1d15 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 11:46:30 -0400
Subject: [PATCH 04/18] fix cpplint errors

---
 numerical_methods/durand_kerner_roots.cpp | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index dd746ad7c..030364d27 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -63,9 +63,10 @@ std::complex<double> poly_function(const std::valarray<double> &coeffs,
  * \returns pointer to converted string
  */
 const char *complex_str(const std::complex<double> &x) {
-    static char msg[50];
+#define MAX_BUFF_SIZE 50
+    static char msg[MAX_BUFF_SIZE];
 
-    sprintf(msg, "% 7.04g%+7.04gj", x.real(), x.imag());
+    std::snprintf(msg, MAX_BUFF_SIZE, "% 7.04g%+7.04gj", x.real(), x.imag());
 
     return msg;
 }
@@ -148,7 +149,10 @@ int main(int argc, char **argv) {
 
     // numerical errors less when the first coefficient is "1"
     // hence, we normalize the first coefficient
-    coeffs /= coeffs[0];
+    {
+        double tmp = coeffs[0];
+        coeffs /= tmp;
+    }
 
 #if defined(DEBUG) || !defined(NDEBUG)
     log_file << "avg. correction";
@@ -157,7 +161,7 @@ int main(int argc, char **argv) {
 #endif
 
     double tol_condition = 1;
-    unsigned long iter = 0;
+    uint32_t iter = 0;
 
     clock_t end_time, start_time = clock();
     while (!check_termination(tol_condition) && iter < INT_MAX) {
@@ -220,7 +224,8 @@ end:
         std::cout << "\t" << complex_str(s0[n]) << "\n";
     std::cout << "absolute average change: " << tol_condition << "\n";
     std::cout << "Time taken: "
-              << (end_time - start_time) / (double)CLOCKS_PER_SEC << " sec\n";
+              << static_cast<double>(end_time - start_time) / CLOCKS_PER_SEC
+              << " sec\n";
 
     return 0;
-}
\ No newline at end of file
+}

From 7a4ca60754d102ceaba06de43273cc7f7ca7680f Mon Sep 17 00:00:00 2001
From: github-actions <${GITHUB_ACTOR}@users.noreply.github.com>
Date: Thu, 4 Jun 2020 15:47:37 +0000
Subject: [PATCH 05/18] updating DIRECTORY.md

---
 DIRECTORY.md | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/DIRECTORY.md b/DIRECTORY.md
index 5507717ef..3bca579f6 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -10,14 +10,6 @@
   * [Rat Maze](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/backtracking/rat_maze.cpp)
   * [Sudoku Solve](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/backtracking/sudoku_solve.cpp)
 
-## Computer Oriented Statistical Methods
-  * [Bisection Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/bisection_method.cpp)
-  * [False Position](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/false_position.cpp)
-  * [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/gaussian_elimination.cpp)
-  * [Newton Raphson Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/newton_raphson_method.cpp)
-  * [Ordinary Least Squares Regressor](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp)
-  * [Successive Approximation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/successive_approximation.cpp)
-
 ## Data Structure
   * [Avltree](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/data_structure/avltree.cpp)
   * [Binary Search Tree](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/data_structure/binary_search_tree.cpp)
@@ -134,6 +126,15 @@
   * [Sqrt Double](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/sqrt_double.cpp)
   * [String Fibonacci](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/string_fibonacci.cpp)
 
+## Numerical Methods
+  * [Bisection Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/bisection_method.cpp)
+  * [Durand Kerner Roots](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/durand_kerner_roots.cpp)
+  * [False Position](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/false_position.cpp)
+  * [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/gaussian_elimination.cpp)
+  * [Newton Raphson Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/newton_raphson_method.cpp)
+  * [Ordinary Least Squares Regressor](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ordinary_least_squares_regressor.cpp)
+  * [Successive Approximation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/successive_approximation.cpp)
+
 ## Operations On Datastructures
   * [Array Left Rotation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/operations_on_datastructures/array_left_rotation.cpp)
   * [Array Right Rotation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/operations_on_datastructures/array_right_rotation.cpp)

From 4d548fe29eb2cd37e40926b1d296aad1f447712b Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 12:14:26 -0400
Subject: [PATCH 06/18] convert to function module

---
 numerical_methods/durand_kerner_roots.cpp | 162 +++++++++++-----------
 1 file changed, 81 insertions(+), 81 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 030364d27..e407a5734 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -85,6 +85,83 @@ bool check_termination(long double delta) {
     return false;
 }
 
+std::pair<uint32_t, double> durand_kerner_algo(
+    const std::valarray<double> &coeffs,
+    std::valarray<std::complex<double>> *roots, bool write_log = false) {
+    double tol_condition = 1;
+    uint32_t iter = 0;
+    int i, n, degree = coeffs.size();
+    std::ofstream log_file;
+
+    if (write_log) {
+        /*
+         * store intermediate values to a CSV file
+         */
+        log_file.open("durand_kerner.log.csv");
+        if (!log_file.is_open()) {
+            perror("Unable to create a storage log file!");
+            std::exit(EXIT_FAILURE);
+        }
+        log_file << "iter#,";
+
+        for (n = 0; n < degree - 1; n++) log_file << "root_" << n << ",";
+
+        log_file << "avg. correction";
+        log_file << "\n0,";
+        for (n = 0; n < degree - 1; n++)
+            log_file << complex_str((*roots)[n]) << ",";
+    }
+
+    while (!check_termination(tol_condition) && iter < INT_MAX) {
+        std::complex<double> delta = 0;
+        tol_condition = 0;
+        iter++;
+
+        if (log_file.is_open())
+            log_file << "\n" << iter << ",";
+
+        for (n = 0; n < degree - 1; n++) {
+            std::complex<double> numerator, denominator;
+            numerator = poly_function(coeffs, (*roots)[n]);
+            denominator = 1.0;
+            for (i = 0; i < degree - 1; i++)
+                if (i != n)
+                    denominator *= (*roots)[n] - (*roots)[i];
+
+            delta = numerator / denominator;
+
+            if (std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
+                std::cerr << "\n\nOverflow/underrun error - got value = "
+                          << std::abs(delta);
+                return std::pair<uint32_t, double>(iter, tol_condition);
+            }
+
+            (*roots)[n] -= delta;
+
+            tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
+
+            if (log_file.is_open())
+                log_file << complex_str((*roots)[n]) << ",";
+        }
+        // tol_condition /= (degree - 1);
+
+#if defined(DEBUG) || !defined(NDEBUG)
+        if (iter % 500 == 0) {
+            std::cout << "Iter: " << iter << "\t";
+            for (n = 0; n < degree - 1; n++)
+                std::cout << "\t" << complex_str((*roots)[n]);
+            std::cout << "\t\tabsolute average change: " << tol_condition
+                      << "\n";
+        }
+#endif
+
+        if (log_file.is_open())
+            log_file << tol_condition;
+    }
+
+    return std::pair<uint32_t, double>(iter, tol_condition);
+}
+
 /***
  * Main function.
  * The comandline input arguments are taken as coeffiecients of a polynomial.
@@ -96,7 +173,7 @@ bool check_termination(long double delta) {
  **/
 int main(int argc, char **argv) {
     unsigned int degree = 0;
-    unsigned int n, i;
+    unsigned int n;
 
     if (argc < 2) {
         std::cout
@@ -111,19 +188,6 @@ int main(int argc, char **argv) {
     /* initialize random seed: */
     std::srand(std::time(nullptr));
 
-#if defined(DEBUG) || !defined(NDEBUG)
-    /*
-     * store intermediate values to a CSV file
-     */
-    std::ofstream log_file;
-    log_file.open("durand_kerner.log.csv");
-    if (!log_file.is_open()) {
-        perror("Unable to create a storage log file!");
-        return EXIT_FAILURE;
-    }
-    log_file << "iter#,";
-#endif
-
     // number of roots = degree - 1
     std::valarray<std::complex<double>> s0(degree - 1);
 
@@ -141,9 +205,6 @@ int main(int argc, char **argv) {
             s0[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
             s0[n] -= 50.f;
             s0[n] /= 50.f;
-#if defined(DEBUG) || !defined(NDEBUG)
-            log_file << "root_" << n << ",";
-#endif
         }
     }
 
@@ -154,75 +215,14 @@ int main(int argc, char **argv) {
         coeffs /= tmp;
     }
 
-#if defined(DEBUG) || !defined(NDEBUG)
-    log_file << "avg. correction";
-    log_file << "\n0,";
-    for (n = 0; n < degree - 1; n++) log_file << complex_str(s0[n]) << ",";
-#endif
-
-    double tol_condition = 1;
-    uint32_t iter = 0;
-
     clock_t end_time, start_time = clock();
-    while (!check_termination(tol_condition) && iter < INT_MAX) {
-        std::complex<double> delta = 0;
-        tol_condition = 0;
-        iter++;
-
-#if defined(DEBUG) || !defined(NDEBUG)
-        log_file << "\n" << iter << ",";
-#endif
-
-        for (n = 0; n < degree - 1; n++) {
-            std::complex<double> numerator, denominator;
-            numerator = poly_function(coeffs, s0[n]);
-            denominator = 1.0;
-            for (i = 0; i < degree - 1; i++)
-                if (i != n)
-                    denominator *= s0[n] - s0[i];
-
-            delta = numerator / denominator;
-
-            if (std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
-                std::cerr << "\n\nOverflow/underrun error - got value = "
-                          << std::abs(delta);
-                goto end;
-            }
-
-            s0[n] -= delta;
-
-            tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
-
-#if defined(DEBUG) || !defined(NDEBUG)
-            log_file << complex_str(s0[n]) << ",";
-#endif
-        }
-        // tol_condition /= (degree - 1);
-
-#if defined(DEBUG) || !defined(NDEBUG)
-        if (iter % 500 == 0) {
-            std::cout << "Iter: " << iter << "\t";
-            for (n = 0; n < degree - 1; n++)
-                std::cout << "\t" << complex_str(s0[n]);
-            std::cout << "\t\tabsolute average change: " << tol_condition
-                      << "\n";
-        }
-
-        log_file << tol_condition;
-#endif
-    }
-end:
-
+    auto result = durand_kerner_algo(coeffs, &s0, true);
     end_time = clock();
 
-#if defined(DEBUG) || !defined(NDEBUG)
-    log_file.close();
-#endif
-
-    std::cout << "\nIterations: " << iter << "\n";
+    std::cout << "\nIterations: " << result.first << "\n";
     for (n = 0; n < degree - 1; n++)
         std::cout << "\t" << complex_str(s0[n]) << "\n";
-    std::cout << "absolute average change: " << tol_condition << "\n";
+    std::cout << "absolute average change: " << result.second << "\n";
     std::cout << "Time taken: "
               << static_cast<double>(end_time - start_time) / CLOCKS_PER_SEC
               << " sec\n";

From 4b68a888eb2aa38a241a208d13b25169a2933b90 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 12:34:52 -0400
Subject: [PATCH 07/18] update documentation

---
 numerical_methods/durand_kerner_roots.cpp | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index e407a5734..424319329 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -1,6 +1,8 @@
 /**
  * @file
- * Compute riots of any given polynomial.
+ * \brief Compute all possible approximate roots of any given polynomial using
+ * [Durand Kerner
+ * algorithm](https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method)
  *
  * Test the algorithm online:
  * https://gist.github.com/kvedala/27f1b0b6502af935f6917673ec43bcd7
@@ -85,6 +87,16 @@ bool check_termination(long double delta) {
     return false;
 }
 
+/**
+ * Implements Durand Kerner iterative algorithm to compute all roots of a
+ * polynomial.
+ *
+ * \param[in] coeffs coefficients of the polynomial
+ * \param[out] roots the computed roots of the polynomial
+ * \param[in] write_log flag whether to save the log file (default = `false`)
+ * \returns pair of values - number of iterations taken and final accuracy
+ * achieved
+ */
 std::pair<uint32_t, double> durand_kerner_algo(
     const std::valarray<double> &coeffs,
     std::valarray<std::complex<double>> *roots, bool write_log = false) {

From 4e75c6d1eac51fe2731a59eb11e099b9a8e9eb99 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 12:35:13 -0400
Subject: [PATCH 08/18] move openmp to main loop

---
 numerical_methods/durand_kerner_roots.cpp | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 424319329..a714e8f1b 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -46,9 +46,9 @@ std::complex<double> poly_function(const std::valarray<double> &coeffs,
     double real = 0.f, imag = 0.f;
     int n;
 
-#ifdef _OPENMP
-#pragma omp parallel for reduction(+ : real, imag)
-#endif
+    // #ifdef _OPENMP
+    // #pragma omp parallel for reduction(+ : real, imag)
+    // #endif
     for (n = 0; n < coeffs.size(); n++) {
         std::complex<double> tmp =
             coeffs[n] * std::pow(x, coeffs.size() - n - 1);
@@ -132,6 +132,9 @@ std::pair<uint32_t, double> durand_kerner_algo(
         if (log_file.is_open())
             log_file << "\n" << iter << ",";
 
+#ifdef _OPENMP
+#pragma omp for
+#endif
         for (n = 0; n < degree - 1; n++) {
             std::complex<double> numerator, denominator;
             numerator = poly_function(coeffs, (*roots)[n]);

From 29ed37d67e65ccf4e02b7b48bd4b41c02a04f800 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 13:15:28 -0400
Subject: [PATCH 09/18] added two test cases

---
 numerical_methods/durand_kerner_roots.cpp | 82 ++++++++++++++++++++---
 1 file changed, 71 insertions(+), 11 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index a714e8f1b..b49a2cef6 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -19,14 +19,13 @@
  */
 
 #include <algorithm>
-#include <climits>
+#include <cassert>
 #include <cmath>
 #include <complex>
 #include <cstdlib>
 #include <ctime>
 #include <fstream>
 #include <iostream>
-#include <string>
 #include <valarray>
 #ifdef _OPENMP
 #include <omp.h>
@@ -37,7 +36,6 @@
 /**
  * Evaluate the value of a polynomial with given coefficients
  * \param[in] coeffs coefficients of the polynomial
- * \param[in] degree degree of polynomial
  * \param[in] x point at which to evaluate the polynomial
  * \returns \f$f(x)\f$
  **/
@@ -177,27 +175,89 @@ std::pair<uint32_t, double> durand_kerner_algo(
     return std::pair<uint32_t, double>(iter, tol_condition);
 }
 
+/**
+ * Self test the algorithm by checking the roots for \f$x^2+4=0\f$ to which the
+ * roots are \f$0 \pm 2i\f$
+ */
+void test1() {
+    const std::valarray<double> coeffs = {1, 0, 4};  // x^2 - 2 = 0
+    std::valarray<std::complex<double>> roots = {
+        std::complex<double>(-100., 100.),
+        std::complex<double>(100., -100.)  // initial approximate roots
+    };
+    std::valarray<std::complex<double>> expected = {
+        std::complex<double>(0., 2.),
+        std::complex<double>(0., -2.)  // known expected roots
+    };
+
+    auto result = durand_kerner_algo(coeffs, &roots, false);
+
+    for (int i = 0; i < roots.size(); i++) {
+        // check if approximations are have < 0.1% error with one of the
+        // expected roots
+        bool err1 = false;
+        for (int j = 0; j < roots.size(); j++)
+            err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
+        assert(err1);
+    }
+
+    std::cout << "Test 1 passed! - " << result.first << " iterations, "
+              << result.second << " accuracy"
+              << "\n";
+}
+
+/**
+ * Self test the algorithm by checking the roots for \f$0.015625x^3-1=0\f$ to
+ * which the roots are \f$(4+0i),\,(-2\pm3.464i)\f$
+ */
+void test2() {
+    const std::valarray<double> coeffs = {// 0.015625 x^3 - 1 = 0
+                                          1. / 64., 0., 0., -1.};
+    std::valarray<std::complex<double>> roots = {
+        std::complex<double>(-100., 100.), std::complex<double>(-101., 101.),
+        std::complex<double>(2., -2.)  // initial approximate roots
+    };
+    std::valarray<std::complex<double>> expected = {
+        std::complex<double>(4., 0.), std::complex<double>(-2., 3.46410162),
+        std::complex<double>(-2., -3.46410162)  // known expected roots
+    };
+
+    auto result = durand_kerner_algo(coeffs, &roots, false);
+
+    for (int i = 0; i < roots.size(); i++) {
+        // check if approximations are have < 0.1% error with one of the
+        // expected roots
+        bool err1 = false;
+        for (int j = 0; j < roots.size(); j++)
+            err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
+        assert(err1);
+    }
+    std::cout << "Test 2 passed! - " << result.first << " iterations, "
+              << result.second << " accuracy"
+              << "\n";
+}
+
 /***
  * Main function.
- * The comandline input arguments are taken as coeffiecients of a polynomial.
- * For example, this command
+ * The comandline input arguments are taken as coeffiecients of a
+ *polynomial. For example, this command
  * ```sh
  * ./durand_kerner_roots 1 0 -4
  * ```
  * will find roots of the polynomial \f$1\cdot x^2 + 0\cdot x^1 + (-4)=0\f$
  **/
 int main(int argc, char **argv) {
-    unsigned int degree = 0;
-    unsigned int n;
+    test1();
+    test2();
 
     if (argc < 2) {
-        std::cout
-            << "Please pass the coefficients of the polynomial as commandline "
-               "arguments.\n";
+        std::cout << "Please pass the coefficients of the polynomial as "
+                     "commandline "
+                     "arguments.\n";
         return 0;
     }
 
-    degree = argc - 1;                     // detected polynomial degree
+    int n, degree = argc - 1;              // detected polynomial degree
     std::valarray<double> coeffs(degree);  // create coefficiencts array
 
     /* initialize random seed: */

From 3b37d06d2a3ed215ff896da75a5204cf5dfa9ddf Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 13:37:34 -0400
Subject: [PATCH 10/18] use INT16_MAX

---
 numerical_methods/durand_kerner_roots.cpp | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index b49a2cef6..90789fd4d 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -100,7 +100,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
     std::valarray<std::complex<double>> *roots, bool write_log = false) {
     double tol_condition = 1;
     uint32_t iter = 0;
-    int i, n, degree = coeffs.size();
+    int i, n;
     std::ofstream log_file;
 
     if (write_log) {
@@ -114,15 +114,15 @@ std::pair<uint32_t, double> durand_kerner_algo(
         }
         log_file << "iter#,";
 
-        for (n = 0; n < degree - 1; n++) log_file << "root_" << n << ",";
+        for (n = 0; n < roots->size(); n++) log_file << "root_" << n << ",";
 
         log_file << "avg. correction";
         log_file << "\n0,";
-        for (n = 0; n < degree - 1; n++)
+        for (n = 0; n < roots->size(); n++)
             log_file << complex_str((*roots)[n]) << ",";
     }
 
-    while (!check_termination(tol_condition) && iter < INT_MAX) {
+    while (!check_termination(tol_condition) && iter < INT16_MAX) {
         std::complex<double> delta = 0;
         tol_condition = 0;
         iter++;
@@ -133,11 +133,11 @@ std::pair<uint32_t, double> durand_kerner_algo(
 #ifdef _OPENMP
 #pragma omp for
 #endif
-        for (n = 0; n < degree - 1; n++) {
+        for (n = 0; n < roots->size(); n++) {
             std::complex<double> numerator, denominator;
             numerator = poly_function(coeffs, (*roots)[n]);
             denominator = 1.0;
-            for (i = 0; i < degree - 1; i++)
+            for (i = 0; i < roots->size(); i++)
                 if (i != n)
                     denominator *= (*roots)[n] - (*roots)[i];
 
@@ -161,7 +161,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
 #if defined(DEBUG) || !defined(NDEBUG)
         if (iter % 500 == 0) {
             std::cout << "Iter: " << iter << "\t";
-            for (n = 0; n < degree - 1; n++)
+            for (n = 0; n < roots->size(); n++)
                 std::cout << "\t" << complex_str((*roots)[n]);
             std::cout << "\t\tabsolute average change: " << tol_condition
                       << "\n";

From 5e34779d9506abf0ac3ac1cae44c586b7e8906b3 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 13:57:57 -0400
Subject: [PATCH 11/18] remove return statement from omp-for loop and use
 "break"

---
 numerical_methods/durand_kerner_roots.cpp | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 90789fd4d..a970829ca 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -122,10 +122,13 @@ std::pair<uint32_t, double> durand_kerner_algo(
             log_file << complex_str((*roots)[n]) << ",";
     }
 
-    while (!check_termination(tol_condition) && iter < INT16_MAX) {
+    bool break_loop = false;
+    while (!check_termination(tol_condition) && iter < INT16_MAX &&
+           !break_loop) {
         std::complex<double> delta = 0;
         tol_condition = 0;
         iter++;
+        break_loop = false;
 
         if (log_file.is_open())
             log_file << "\n" << iter << ",";
@@ -146,9 +149,13 @@ std::pair<uint32_t, double> durand_kerner_algo(
             if (std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
                 std::cerr << "\n\nOverflow/underrun error - got value = "
                           << std::abs(delta);
-                return std::pair<uint32_t, double>(iter, tol_condition);
+                // return std::pair<uint32_t, double>(iter, tol_condition);
+                break_loop = true;
             }
 
+            if (break_loop)
+                break;
+
             (*roots)[n] -= delta;
 
             tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
@@ -158,6 +165,9 @@ std::pair<uint32_t, double> durand_kerner_algo(
         }
         // tol_condition /= (degree - 1);
 
+        if (break_loop)
+            break;
+
 #if defined(DEBUG) || !defined(NDEBUG)
         if (iter % 500 == 0) {
             std::cout << "Iter: " << iter << "\t";

From 586ac1aee190bb21d2115811fefc9f9c867f2c46 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 13:59:24 -0400
Subject: [PATCH 12/18] run tests when no input is provided and skip tests when
 input polynomial is provided

---
 numerical_methods/durand_kerner_roots.cpp | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index a970829ca..13e79dbd9 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -257,10 +257,9 @@ void test2() {
  * will find roots of the polynomial \f$1\cdot x^2 + 0\cdot x^1 + (-4)=0\f$
  **/
 int main(int argc, char **argv) {
-    test1();
-    test2();
-
     if (argc < 2) {
+        test1();  // run tests when no input is provided
+        test2();  // and skip tests when input polynomial is provided
         std::cout << "Please pass the coefficients of the polynomial as "
                      "commandline "
                      "arguments.\n";

From 0e5d841e0a18c897685a383a6e1b31bf7f56832c Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 14:02:18 -0400
Subject: [PATCH 13/18] while loop cannot have break - replaced with continue
 and check is present in the main while condition

---
 numerical_methods/durand_kerner_roots.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 13e79dbd9..dd6a10f22 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -166,7 +166,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
         // tol_condition /= (degree - 1);
 
         if (break_loop)
-            break;
+            continue;
 
 #if defined(DEBUG) || !defined(NDEBUG)
         if (iter % 500 == 0) {

From 89f33a49ba2bf70e4189c484a60eea59a2d5b286 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 14:14:46 -0400
Subject: [PATCH 14/18] (1) break while loop (2) skip runs on break_loop
 instead of hard-break

---
 numerical_methods/durand_kerner_roots.cpp | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index dd6a10f22..01fd4ebaf 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -134,9 +134,12 @@ std::pair<uint32_t, double> durand_kerner_algo(
             log_file << "\n" << iter << ",";
 
 #ifdef _OPENMP
-#pragma omp for
+#pragma omp parallel for shared(break_loop)
 #endif
         for (n = 0; n < roots->size(); n++) {
+            if (break_loop)
+                continue;
+
             std::complex<double> numerator, denominator;
             numerator = poly_function(coeffs, (*roots)[n]);
             denominator = 1.0;
@@ -153,9 +156,6 @@ std::pair<uint32_t, double> durand_kerner_algo(
                 break_loop = true;
             }
 
-            if (break_loop)
-                break;
-
             (*roots)[n] -= delta;
 
             tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
@@ -166,7 +166,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
         // tol_condition /= (degree - 1);
 
         if (break_loop)
-            continue;
+            break;
 
 #if defined(DEBUG) || !defined(NDEBUG)
         if (iter % 500 == 0) {

From b7d05eab9e99988e2a85612c27d87d0d8740fd1a Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 14:14:59 -0400
Subject: [PATCH 15/18] add documentation images

---
 numerical_methods/durand_kerner_roots.cpp | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 01fd4ebaf..2a8d8b5e7 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -16,6 +16,12 @@
  * -12870931245150988800 13803759753640704000 -8752948036761600000
  * 2432902008176640000
  * ```
+ * Sample implementation results to compute approximate roots of the equation
+ * \f$x^4-1=0\f$:\n
+ * ![Error
+ * evolution](https://raw.githubusercontent.com/kvedala/C-Plus-Plus/docs/images/numerical_methods/durand_kerner_error.svg)
+ * ![Roots
+ * evolution](https://raw.githubusercontent.com/kvedala/C-Plus-Plus/docs/images/numerical_methods/durand_kerner_roots.svg)
  */
 
 #include <algorithm>

From 42f78bbe1998bac1f9b5bd8554462ba647a17141 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 15:20:31 -0400
Subject: [PATCH 16/18] use long double for errors and tolerance checks

---
 numerical_methods/durand_kerner_roots.cpp | 44 ++++++++++++++---------
 1 file changed, 28 insertions(+), 16 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 2a8d8b5e7..5d5b65bf5 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -104,7 +104,7 @@ bool check_termination(long double delta) {
 std::pair<uint32_t, double> durand_kerner_algo(
     const std::valarray<double> &coeffs,
     std::valarray<std::complex<double>> *roots, bool write_log = false) {
-    double tol_condition = 1;
+    long double tol_condition = 1;
     uint32_t iter = 0;
     int i, n;
     std::ofstream log_file;
@@ -153,11 +153,11 @@ std::pair<uint32_t, double> durand_kerner_algo(
                 if (i != n)
                     denominator *= (*roots)[n] - (*roots)[i];
 
-            delta = numerator / denominator;
+            std::complex<long double> delta = numerator / denominator;
 
             if (std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
                 std::cerr << "\n\nOverflow/underrun error - got value = "
-                          << std::abs(delta);
+                          << std::abs(delta) << "\n";
                 // return std::pair<uint32_t, double>(iter, tol_condition);
                 break_loop = true;
             }
@@ -188,7 +188,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
             log_file << tol_condition;
     }
 
-    return std::pair<uint32_t, double>(iter, tol_condition);
+    return std::pair<uint32_t, long double>(iter, tol_condition);
 }
 
 /**
@@ -197,15 +197,19 @@ std::pair<uint32_t, double> durand_kerner_algo(
  */
 void test1() {
     const std::valarray<double> coeffs = {1, 0, 4};  // x^2 - 2 = 0
-    std::valarray<std::complex<double>> roots = {
-        std::complex<double>(-100., 100.),
-        std::complex<double>(100., -100.)  // initial approximate roots
-    };
+    std::valarray<std::complex<double>> roots(2);
     std::valarray<std::complex<double>> expected = {
         std::complex<double>(0., 2.),
         std::complex<double>(0., -2.)  // known expected roots
     };
 
+    /* initialize root approximations with random values */
+    for (int n = 0; n < roots.size(); n++) {
+        roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
+        roots[n] -= 50.f;
+        roots[n] /= 25.f;
+    }
+
     auto result = durand_kerner_algo(coeffs, &roots, false);
 
     for (int i = 0; i < roots.size(); i++) {
@@ -229,17 +233,24 @@ void test1() {
 void test2() {
     const std::valarray<double> coeffs = {// 0.015625 x^3 - 1 = 0
                                           1. / 64., 0., 0., -1.};
-    std::valarray<std::complex<double>> roots = {
-        std::complex<double>(-100., 100.), std::complex<double>(-101., 101.),
-        std::complex<double>(2., -2.)  // initial approximate roots
-    };
-    std::valarray<std::complex<double>> expected = {
+    std::valarray<std::complex<double>> roots(3);
+    const std::valarray<std::complex<double>> expected = {
         std::complex<double>(4., 0.), std::complex<double>(-2., 3.46410162),
         std::complex<double>(-2., -3.46410162)  // known expected roots
     };
 
+    /* initialize root approximations with random values */
+    for (int n = 0; n < roots.size(); n++) {
+        roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
+        roots[n] -= 50.f;
+        roots[n] /= 25.f;
+    }
+
     auto result = durand_kerner_algo(coeffs, &roots, false);
 
+    for (int n = 0; n < roots.size(); n++)
+        std::cout << "\t" << complex_str(roots[n]) << "\n";
+
     for (int i = 0; i < roots.size(); i++) {
         // check if approximations are have < 0.1% error with one of the
         // expected roots
@@ -248,6 +259,7 @@ void test2() {
             err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
         assert(err1);
     }
+
     std::cout << "Test 2 passed! - " << result.first << " iterations, "
               << result.second << " accuracy"
               << "\n";
@@ -263,6 +275,9 @@ void test2() {
  * will find roots of the polynomial \f$1\cdot x^2 + 0\cdot x^1 + (-4)=0\f$
  **/
 int main(int argc, char **argv) {
+    /* initialize random seed: */
+    std::srand(std::time(nullptr));
+
     if (argc < 2) {
         test1();  // run tests when no input is provided
         test2();  // and skip tests when input polynomial is provided
@@ -275,9 +290,6 @@ int main(int argc, char **argv) {
     int n, degree = argc - 1;              // detected polynomial degree
     std::valarray<double> coeffs(degree);  // create coefficiencts array
 
-    /* initialize random seed: */
-    std::srand(std::time(nullptr));
-
     // number of roots = degree - 1
     std::valarray<std::complex<double>> s0(degree - 1);
 

From 645321c38fa65bf98a0e9d52f14b5f2834ed7f87 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 15:21:09 -0400
Subject: [PATCH 17/18] make iterator variable i local to threads

---
 numerical_methods/durand_kerner_roots.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 5d5b65bf5..395797bfc 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -106,7 +106,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
     std::valarray<std::complex<double>> *roots, bool write_log = false) {
     long double tol_condition = 1;
     uint32_t iter = 0;
-    int i, n;
+    int n;
     std::ofstream log_file;
 
     if (write_log) {
@@ -149,7 +149,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
             std::complex<double> numerator, denominator;
             numerator = poly_function(coeffs, (*roots)[n]);
             denominator = 1.0;
-            for (i = 0; i < roots->size(); i++)
+            for (int i = 0; i < roots->size(); i++)
                 if (i != n)
                     denominator *= (*roots)[n] - (*roots)[i];
 

From e4f53974989b7317ade3e85c2c4d9d6b0e62c203 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 4 Jun 2020 15:21:45 -0400
Subject: [PATCH 18/18] add critical secions to omp threads

---
 numerical_methods/durand_kerner_roots.cpp | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp
index 395797bfc..e92c8d6b9 100644
--- a/numerical_methods/durand_kerner_roots.cpp
+++ b/numerical_methods/durand_kerner_roots.cpp
@@ -140,7 +140,7 @@ std::pair<uint32_t, double> durand_kerner_algo(
             log_file << "\n" << iter << ",";
 
 #ifdef _OPENMP
-#pragma omp parallel for shared(break_loop)
+#pragma omp parallel for shared(break_loop, tol_condition)
 #endif
         for (n = 0; n < roots->size(); n++) {
             if (break_loop)
@@ -164,9 +164,15 @@ std::pair<uint32_t, double> durand_kerner_algo(
 
             (*roots)[n] -= delta;
 
+#ifdef _OPENMP
+#pragma omp critical
+#endif
             tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
 
             if (log_file.is_open())
+#ifdef _OPENMP
+#pragma omp critical
+#endif
                 log_file << complex_str((*roots)[n]) << ",";
         }
         // tol_condition /= (degree - 1);
@@ -248,9 +254,6 @@ void test2() {
 
     auto result = durand_kerner_algo(coeffs, &roots, false);
 
-    for (int n = 0; n < roots.size(); n++)
-        std::cout << "\t" << complex_str(roots[n]) << "\n";
-
     for (int i = 0; i < roots.size(); i++) {
         // check if approximations are have < 0.1% error with one of the
         // expected roots