aberth.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <complex>
#include <lds/lds.hpp>
#include <utility>
#include <vector>
 
namespace ginger {
    class Options;
}
 
template <unsigned long N, unsigned long Base = 2> constexpr auto make_vdc_table()
    -> std::array<double, N> {
    std::array<double, N> table{};
    lds::VdCorput<Base> gen;
    for (unsigned long i = 0; i < N; ++i) {
        table[i] = gen.pop();
    }
    return table;
}
 
constexpr const auto VDC_TABLE_SIZE = 1000UL;
 
constexpr std::array<double, VDC_TABLE_SIZE> VDC_TABLE_2 = make_vdc_table<VDC_TABLE_SIZE, 2>();
 
inline double vdc2_table(unsigned long index) { return VDC_TABLE_2[index]; }
 
inline auto make_circle_table() noexcept -> std::array<std::array<double, 2>, VDC_TABLE_SIZE> {
    std::array<std::array<double, 2>, VDC_TABLE_SIZE> table{};
    for (unsigned long i = 0; i < VDC_TABLE_SIZE; ++i) {
        auto theta = VDC_TABLE_2[i] * lds::TWO_PI;
        table[i] = {std::cos(theta), std::sin(theta)};
    }
    return table;
}
 
static const auto CIRCLE_TABLE_2 = make_circle_table();
 
inline auto make_cos_pi_vdc2_table() noexcept -> std::array<double, VDC_TABLE_SIZE> {
    std::array<double, VDC_TABLE_SIZE> table{};
    for (unsigned long i = 0; i < VDC_TABLE_SIZE; ++i) {
        table[i] = std::cos(lds::TWO_PI / 2.0 * VDC_TABLE_2[i]);
    }
    return table;
}
 
static const auto COS_PI_VDC2_TABLE = make_cos_pi_vdc2_table();
 
constexpr double cos_pi_vdc2(unsigned long index) { return COS_PI_VDC2_TABLE[index]; }
 
constexpr double circle2_table_x(unsigned long index) { return CIRCLE_TABLE_2[index][0]; }
 
constexpr double circle2_table_y(unsigned long index) { return CIRCLE_TABLE_2[index][1]; }
 
inline auto horner_eval_c(const std::vector<double>& coeffs, const std::complex<double>& zval)
    -> std::complex<double> {
    std::complex<double> result(0.0, 0.0);
    for (auto coeff : coeffs) {
        result = result * zval + coeff;
    }
    return result;
}
 
inline auto horner_eval_f(const std::vector<double>& coeffs, const double& zval) -> double {
    double result(0.0);
    for (auto coeff : coeffs) {
        result = result * zval + coeff;
    }
    return result;
}
 
extern auto initial_aberth(const std::vector<double>& coeffs) -> std::vector<std::complex<double>>;
 
extern auto aberth(const std::vector<double>& coeffs, std::vector<std::complex<double>>& zs,
                   const ginger::Options& options) -> std::pair<unsigned int, bool>;
 
extern auto initial_aberth_autocorr(const std::vector<double>& coeffs)
    -> std::vector<std::complex<double>>;
 
extern auto aberth_autocorr(const std::vector<double>& coeffs,
                            std::vector<std::complex<double>>& zs, const ginger::Options& options)
    -> std::pair<unsigned int, bool>;
 
extern auto poly_from_roots(const std::vector<std::complex<double>>& zs) -> std::vector<double>;
 
extern auto leja_order(const std::vector<std::complex<double>>& points)
    -> std::vector<std::complex<double>>;
 
extern auto poly_from_autocorr_roots(const std::vector<std::complex<double>>& zs)
    -> std::vector<double>;