std::complex<double>
poly_function(
conststd::valarray<double> &coeffs,
55std::complex<double> x) {
56 doublereal = 0.f, imag = 0.f;
62 for(n = 0; n < coeffs.size(); n++) {
63std::complex<double> tmp =
64coeffs[n] * std::pow(x, coeffs.size() - n - 1);
69 returnstd::complex<double>(real, imag);
78#define MAX_BUFF_SIZE 50 79 static charmsg[MAX_BUFF_SIZE];
81std::snprintf(msg, MAX_BUFF_SIZE,
"% 7.04g%+7.04gj", x.real(), x.imag());
93 static long doublepast_delta = INFINITY;
111 conststd::valarray<double> &coeffs,
112std::valarray<std::complex<double>> *roots,
boolwrite_log =
false) {
113 long doubletol_condition = 1;
116std::ofstream log_file;
122log_file.open(
"durand_kerner.log.csv");
123 if(!log_file.is_open()) {
124perror(
"Unable to create a storage log file!");
125std::exit(EXIT_FAILURE);
127log_file <<
"iter#,";
129 for(n = 0; n < roots->size(); n++) log_file <<
"root_"<< n <<
",";
131log_file <<
"avg. correction";
133 for(n = 0; n < roots->size(); n++)
137 boolbreak_loop =
false;
144 if(log_file.is_open())
145log_file <<
"\n"<< iter <<
",";
148#pragma omp parallel for shared(break_loop, tol_condition) 150 for(n = 0; n < roots->size(); n++) {
154std::complex<double> numerator, denominator;
157 for(
inti = 0; i < roots->size(); i++)
159denominator *= (*roots)[n] - (*roots)[i];
161std::complex<long double> delta = numerator / denominator;
163 if(std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
164std::cerr <<
"\n\nOverflow/underrun error - got value = " 165<< std::abs(delta) <<
"\n";
170(*roots)[n] -= delta;
175tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
182 if(log_file.is_open()) {
183 for(n = 0; n < roots->size(); n++)
187#if defined(DEBUG) || !defined(NDEBUG) 188 if(iter % 500 == 0) {
189std::cout <<
"Iter: "<< iter <<
"\t";
190 for(n = 0; n < roots->size(); n++)
192std::cout <<
"\t\tabsolute average change: "<< tol_condition
197 if(log_file.is_open())
198log_file << tol_condition;
201 returnstd::pair<uint32_t, long double>(iter, tol_condition);
209 conststd::valarray<double> coeffs = {1, 0, 4};
210std::valarray<std::complex<double>> roots(2);
211std::valarray<std::complex<double>> expected = {
212std::complex<double>(0., 2.),
213std::complex<double>(0., -2.)
217 for(
intn = 0; n < roots.size(); n++) {
218roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
225 for(
inti = 0; i < roots.size(); i++) {
229 for(
intj = 0; j < roots.size(); j++)
230err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
234std::cout <<
"Test 1 passed! - "<< result.first <<
" iterations, " 235<< result.second <<
" accuracy" 244 conststd::valarray<double> coeffs = {
2451. / 64., 0., 0., -1.};
246std::valarray<std::complex<double>> roots(3);
247 conststd::valarray<std::complex<double>> expected = {
248std::complex<double>(4., 0.), std::complex<double>(-2., 3.46410162),
249std::complex<double>(-2., -3.46410162)
253 for(
intn = 0; n < roots.size(); n++) {
254roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
261 for(
inti = 0; i < roots.size(); i++) {
265 for(
intj = 0; j < roots.size(); j++)
266err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
270std::cout <<
"Test 2 passed! - "<< result.first <<
" iterations, " 271<< result.second <<
" accuracy" 284int main(
intargc,
char**argv) {
286std::srand(std::time(
nullptr));
291std::cout <<
"Please pass the coefficients of the polynomial as " 297 intn, degree = argc - 1;
298std::valarray<double> coeffs(degree);
301std::valarray<std::complex<double>> s0(degree - 1);
303std::cout <<
"Computing the roots for:\n\t";
304 for(n = 0; n < degree; n++) {
305coeffs[n] = strtod(argv[n + 1],
nullptr);
306 if(n < degree - 1 && coeffs[n] != 0)
307std::cout <<
"("<< coeffs[n] <<
") x^"<< degree - n - 1 <<
" + ";
308 else if(coeffs[n] != 0)
309std::cout <<
"("<< coeffs[n] <<
") x^"<< degree - n - 1
313 if(n < degree - 1) {
314s0[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
323 doubletmp = coeffs[0];
327clock_t end_time, start_time = clock();
331std::cout <<
"\nIterations: "<<
result.first <<
"\n";
332 for(n = 0; n < degree - 1; n++)
334std::cout <<
"absolute average change: "<<
result.second <<
"\n";
335std::cout <<
"Time taken: " 336<<
static_cast<double>(end_time - start_time) / CLOCKS_PER_SEC
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)
const char * complex_str(const std::complex< double > &x)
std::complex< double > poly_function(const std::valarray< double > &coeffs, std::complex< double > x)
uint64_t result(uint64_t n)
int main()
Main function.
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4