rootfinding.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <utility>
#include <vector>
 
#include "matrix2.hpp"
#include "vector2.hpp"
 
using Vec2 = ginger::Vector2<double>;
using Mat2 = ginger::Matrix2<Vec2>;
 
class Options;
 
extern auto initial_guess(std::vector<double> coeffs) -> std::vector<Vec2>;
 
extern auto pbairstow_even_st(const std::vector<double>& coeffs, std::vector<Vec2>& vrs,
                              const Options& options) -> std::pair<unsigned int, bool>;
 
extern auto pbairstow_even_mt(const std::vector<double>& coeffs, std::vector<Vec2>& vrs,
                              const Options& options) -> std::pair<unsigned int, bool>;
 
inline auto pbairstow_even(const std::vector<double>& coeffs, std::vector<Vec2>& vrs,
                           const Options& options) -> std::pair<unsigned int, bool> {
    return pbairstow_even_st(coeffs, vrs, options);
}
 
extern auto horner(std::vector<double>& coeffs1, std::size_t degree, const Vec2& vr) -> Vec2;
 
extern auto suppress_old(Vec2& vA, Vec2& vA1, const Vec2& vri, const Vec2& vrj) -> void;
 
extern auto suppress(Vec2& vA, Vec2& vA1, const Vec2& vri, const Vec2& vrj) -> void;
 
extern auto suppress2(Vec2& vA, Vec2& vA1, const Vec2& vri, const Vec2& vrj) -> void;
 
inline auto makeadjoint(const Vec2& vr, const Vec2& vp) -> Mat2 {
    auto&& p = vp.x();
    auto&& s = vp.y();
    return {Vec2{s, -p}, Vec2{-p * vr.y(), p * vr.x() + s}};
}
 
inline auto delta_scalar(const Vec2& vA, const Vec2& vr, const Vec2& vp) -> Vec2 {
    const auto r = vr.x();
    const auto q = vr.y();
    const auto p = vp.x();
    const auto s = vp.y();
    // adj = [[s, -p], [-p*q, p*r + s]]
    // det = s*(p*r + s) - p*p*q
    const auto det = s * (p * r + s) - p * p * q;
    // mdot_x = s*vA.x + (-p)*vA.y
    // mdot_y = (-p*q)*vA.x + (p*r + s)*vA.y
    const auto mx = s * vA.x() - p * vA.y();
    const auto my = -p * q * vA.x() + (p * r + s) * vA.y();
    return Vec2{mx / det, my / det};
}
 
inline auto delta(const Vec2& vA, const Vec2& vr, const Vec2& vp) -> Vec2 {
    const auto mp = makeadjoint(vr, vp);  // 2 mul's
    return mp.mdot(vA) / mp.det();        // 6 mul's + 2 div's
}
 
inline auto horner_eval(std::vector<double> coeffs1, std::size_t degree, const double& z)
    -> double {
    for (auto idx = 0U; idx != degree; ++idx) {
        coeffs1[idx + 1] += coeffs1[idx] * z;
    }
    return coeffs1[degree];
}
 
extern auto poly_from_quadratic_factors(const std::vector<Vec2>& vrs) -> std::vector<double>;