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/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) 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/numerical_methods/durand_kerner_roots.cpp b/numerical_methods/durand_kerner_roots.cpp new file mode 100644 index 000000000..e92c8d6b9 --- /dev/null +++ b/numerical_methods/durand_kerner_roots.cpp @@ -0,0 +1,336 @@ +/** + * @file + * \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 + * + * 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 + * ``` + * 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 +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef _OPENMP +#include +#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] x point at which to evaluate the polynomial + * \returns \f$f(x)\f$ + **/ +std::complex poly_function(const std::valarray &coeffs, + std::complex 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 tmp = + coeffs[n] * std::pow(x, coeffs.size() - n - 1); + real += tmp.real(); + imag += tmp.imag(); + } + + return std::complex(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 &x) { +#define MAX_BUFF_SIZE 50 + static char msg[MAX_BUFF_SIZE]; + + std::snprintf(msg, MAX_BUFF_SIZE, "% 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; +} + +/** + * 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 durand_kerner_algo( + const std::valarray &coeffs, + std::valarray> *roots, bool write_log = false) { + long double tol_condition = 1; + uint32_t iter = 0; + int n; + 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 < roots->size(); n++) log_file << "root_" << n << ","; + + log_file << "avg. correction"; + log_file << "\n0,"; + for (n = 0; n < roots->size(); n++) + log_file << complex_str((*roots)[n]) << ","; + } + + bool break_loop = false; + while (!check_termination(tol_condition) && iter < INT16_MAX && + !break_loop) { + std::complex delta = 0; + tol_condition = 0; + iter++; + break_loop = false; + + if (log_file.is_open()) + log_file << "\n" << iter << ","; + +#ifdef _OPENMP +#pragma omp parallel for shared(break_loop, tol_condition) +#endif + for (n = 0; n < roots->size(); n++) { + if (break_loop) + continue; + + std::complex numerator, denominator; + numerator = poly_function(coeffs, (*roots)[n]); + denominator = 1.0; + for (int i = 0; i < roots->size(); i++) + if (i != n) + denominator *= (*roots)[n] - (*roots)[i]; + + std::complex 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) << "\n"; + // return std::pair(iter, tol_condition); + break_loop = true; + } + + (*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); + + if (break_loop) + break; + +#if defined(DEBUG) || !defined(NDEBUG) + if (iter % 500 == 0) { + std::cout << "Iter: " << iter << "\t"; + for (n = 0; n < roots->size(); 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(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 coeffs = {1, 0, 4}; // x^2 - 2 = 0 + std::valarray> roots(2); + std::valarray> expected = { + std::complex(0., 2.), + std::complex(0., -2.) // known expected roots + }; + + /* initialize root approximations with random values */ + for (int n = 0; n < roots.size(); n++) { + roots[n] = std::complex(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++) { + // 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 coeffs = {// 0.015625 x^3 - 1 = 0 + 1. / 64., 0., 0., -1.}; + std::valarray> roots(3); + const std::valarray> expected = { + std::complex(4., 0.), std::complex(-2., 3.46410162), + std::complex(-2., -3.46410162) // known expected roots + }; + + /* initialize root approximations with random values */ + for (int n = 0; n < roots.size(); n++) { + roots[n] = std::complex(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++) { + // 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 + * ```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) { + /* 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 + std::cout << "Please pass the coefficients of the polynomial as " + "commandline " + "arguments.\n"; + return 0; + } + + int n, degree = argc - 1; // detected polynomial degree + std::valarray coeffs(degree); // create coefficiencts array + + // number of roots = degree - 1 + std::valarray> 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(std::rand() % 100, std::rand() % 100); + s0[n] -= 50.f; + s0[n] /= 50.f; + } + } + + // numerical errors less when the first coefficient is "1" + // hence, we normalize the first coefficient + { + double tmp = coeffs[0]; + coeffs /= tmp; + } + + clock_t end_time, start_time = clock(); + auto result = durand_kerner_algo(coeffs, &s0, true); + end_time = clock(); + + 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: " << result.second << "\n"; + std::cout << "Time taken: " + << static_cast(end_time - start_time) / CLOCKS_PER_SEC + << " sec\n"; + + return 0; +} 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