58 lines
1.5 KiB
C++
58 lines
1.5 KiB
C++
#include "rooting.hpp"
|
|
#include <cmath>
|
|
#include <stdexcept>
|
|
#include <limits>
|
|
|
|
double rooting::nth_root(long double number, double n) {
|
|
// Special case handling
|
|
if (n <= 0) {
|
|
throw std::invalid_argument("The root must be positive and non-zero");
|
|
}
|
|
|
|
if (number < 0 && std::fmod(n, 2.0) == 0) {
|
|
throw std::invalid_argument("Cannot compute even root of a negative number");
|
|
}
|
|
|
|
if (number == 0) {
|
|
return 0.0;
|
|
}
|
|
|
|
if (number == 1 || n == 1.0) {
|
|
return number;
|
|
}
|
|
|
|
// Handle negative numbers for odd roots
|
|
bool negative = number < 0;
|
|
double abs_number = negative ? -number : number;
|
|
|
|
// Use logarithm method for initial guess
|
|
double log_result = std::log(abs_number) / n;
|
|
double x = std::exp(log_result);
|
|
|
|
// Constants
|
|
const double epsilon = 1e-15;
|
|
const int max_iterations = 100;
|
|
|
|
// Newton-Raphson method
|
|
for (int i = 0; i < max_iterations; i++) {
|
|
double x_prev = x;
|
|
|
|
// Calculate new approximation, handling potential overflow
|
|
double x_pow = std::pow(x, n - 1);
|
|
|
|
// Prevent division by zero
|
|
if (x_pow < std::numeric_limits<double>::min()) {
|
|
break;
|
|
}
|
|
|
|
x = ((n - 1) * x + abs_number / x_pow) / n;
|
|
|
|
// Check for convergence using relative error
|
|
if (std::fabs(x - x_prev) <= epsilon * std::fabs(x)) {
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Return correct sign for odd roots of negative numbers
|
|
return negative ? -x : x;
|
|
} |