#pragma once
#include <memory>
#include <limits>


template <typename T>
struct Complex {
    T re {0}, im {0};
    inline Complex<T> operator+(const T& other) const {
        return Complex<T>{other + re, im};
    }
    inline Complex<T> operator*(const T& other) const {
        return Complex<T>{other*re, other*im};
    }
    inline Complex<T> operator*(const Complex<T>& other) const {
        return Complex<T>{re*other.re - im*other.im, re*other.im + im*other.re};
    }
};

template <typename T>
struct RPoly {
    /*  static constexpr <-- qdlib not constexpr friendly */ T eta {std::numeric_limits<T>::epsilon()};
    const int degree;
private:
    T sr,si,u,v,a,b,c,d,a1,a2;
    T a3,a6,a7,e,f,g,h,szr,szi,lzr,lzi;
    T are,mre;
    int n,nn,nmi,zerok;

    std::unique_ptr<T[]> buffer;
    T * const __restrict temp;
    T * const __restrict pt;
    T * const __restrict p;
    T * const __restrict qp;
    T * const __restrict k;
    T * const __restrict qk;
    T * const __restrict svk;

public:
    RPoly(int degree);
    int findRoots(const T * const coeffs, T * const zeror, T * const zeroi);
    T refineRealRoot(const T * const coeffs, const T& rr) const;
    T evalPoly(const T * const coeffs, T x) const;
    Complex<T> evalPoly(const T * const coeffs, T re, T im) const;
private:
    void quad(T a,T b1,T c,T *sr,T *si,
            T *lr,T *li);
    void fxshfr(int l2, int *nz);
    void quadit(T *uu,T *vv,int *nz);
    void realit(T *sss, int *nz, int *iflag);
    void calcsc(int *type);
    void nextk(int *type);
    void newest(int type,T *uu,T *vv);
    void quadsd(int n,T *u,T *v,T *p,T *q, T *a,T *b);
};