tpl_polynomial.H

Go to the documentation of this file.

 
/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
#ifndef TPL_POLYNOMIAL_H
#define TPL_POLYNOMIAL_H
 
#include <cmath>
#include <initializer_list>
#include <limits>
#include <sstream>
#include <string>
#include <type_traits>
#include <utility>
#include <numeric>  // std::gcd
#include <random>   // std::mt19937_64
 
#include <ah-errors.H>
#include <tpl_dynMapTree.H>
#include <tpl_array.H>
#include <htlist.H>
#include <modular_arithmetic.H>  // mod_inv, mod_mul
 
namespace Aleph {
 
namespace polynomial_detail {
 
template <typename C>
C power(C base, size_t exp)
{
  C result = C(1);
  while (exp > 0)
    {
      if (exp & 1)
        result = result * base;
      base = base * base;
      exp >>= 1;
    }
  return result;
}
 
template <typename Int>
  requires std::is_integral_v<Int>
[[nodiscard]] uint64_t abs_to_u64(Int value) noexcept
{
  if constexpr (std::is_signed_v<Int>)
    {
      using Unsigned = std::make_unsigned_t<Int>;
      const Unsigned magnitude = value < Int{} ? Unsigned{} - static_cast<Unsigned>(value)
                                               : static_cast<Unsigned>(value);
      return static_cast<uint64_t>(magnitude);
    }
  else
    return static_cast<uint64_t>(value);
}
 
template <typename Int>
  requires std::is_integral_v<Int>
[[nodiscard]] uint64_t normalize_mod_u64(Int value, const uint64_t modulus) noexcept
{
  if (modulus == 0)
    return 0;
 
  if constexpr (std::is_signed_v<Int>)
    {
      const auto mod_i64 = static_cast<int64_t>(modulus);
      auto rem = static_cast<int64_t>(value % static_cast<Int>(mod_i64));
      if (rem < 0)
        rem += mod_i64;
      return static_cast<uint64_t>(rem);
    }
  else
    return static_cast<uint64_t>(value % static_cast<Int>(modulus));
}
 
template <typename Int>
  requires std::is_integral_v<Int>
[[nodiscard]] Int centered_from_mod_u64(const uint64_t value, const uint64_t modulus) noexcept
{
  if constexpr (std::is_signed_v<Int>)
    {
      const uint64_t canonical = modulus == 0 ? 0 : value % modulus;
      const uint64_t half = modulus / 2;
      if (canonical > half)
        return static_cast<Int>(static_cast<int64_t>(canonical) - static_cast<int64_t>(modulus));
      return static_cast<Int>(canonical);
    }
  else
    return static_cast<Int>(modulus == 0 ? 0 : value % modulus);
}
 
}  // end namespace polynomial_detail
 
template <typename Coefficient = double>
class Gen_Polynomial
{
public:
  using Coeff = Coefficient;
 
private:
  using Term = std::pair<size_t, Coefficient>;
  DynMapTree<size_t, Coefficient> coeffs;
 
  // --- Epsilon-based zero detection for floating-point types ----------
 
  static constexpr Coefficient epsilon() noexcept
  {
    if constexpr (std::is_floating_point_v<Coefficient>)
      return Coefficient(128) * std::numeric_limits<Coefficient>::epsilon();
    else
      return Coefficient{};
  }
 
  static bool coeff_is_zero(const Coefficient &c) noexcept
  {
    if constexpr (std::is_floating_point_v<Coefficient>)
      return std::abs(c) <= epsilon();
    else
      return c == Coefficient{};
  }
 
  void remove_zeros()
  {
    DynList<size_t> zeros;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        if (coeff_is_zero(p.second))
          zeros.append(p.first);
      }
    zeros.for_each([this](size_t e)
    {
      coeffs.remove_key(e);
    });
  }
 
  Coefficient coeff_at(size_t exp) const
  {
    auto *p = coeffs.search(exp);
    return p != nullptr ? p->second : Coefficient{};
  }
 
  Coefficient *coeff_ptr(size_t exp) noexcept
  {
    auto *p = coeffs.search(exp);
    return p != nullptr ? &p->second : nullptr;
  }
 
  const Coefficient *coeff_ptr(size_t exp) const noexcept
  {
    auto *p = coeffs.search(exp);
    return p != nullptr ? &p->second : nullptr;
  }
 
  void add_to_coeff(size_t exp, const Coefficient &delta)
  {
    if (coeff_is_zero(delta))
      return;
 
    if (auto *slot = coeff_ptr(exp); slot != nullptr)
      {
        *slot = *slot + delta;
        if (coeff_is_zero(*slot))
          coeffs.remove_key(exp);
      }
    else
      coeffs.insert(exp, delta);
  }
 
  void add_scaled_shifted(const Gen_Polynomial &src, const Coefficient &scale, size_t shift = 0)
  {
    if (coeff_is_zero(scale) or src.is_zero())
      return;
 
    src.for_each_term([this, &scale, shift](size_t exp, const Coefficient &c)
    {
      add_to_coeff(exp + shift, c * scale);
    });
  }
 
  void scale_inplace(const Coefficient &s)
  {
    if (coeffs.is_empty())
      return;
 
    if (coeff_is_zero(s))
      {
        coeffs = DynMapTree<size_t, Coefficient>();
        return;
      }
 
    if (s == Coefficient(1))
      return;
 
    DynList<size_t> zeros;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        p.second = p.second * s;
        if (coeff_is_zero(p.second))
          zeros.append(p.first);
      }
 
    zeros.for_each([this](size_t exp)
    {
      coeffs.remove_key(exp);
    });
  }
 
  void divide_scalar_inplace(const Coefficient &s)
  {
    ah_domain_error_if(coeff_is_zero(s)) << "polynomial division by zero scalar";
 
    if (coeffs.is_empty())
      return;
 
    if (s == Coefficient(1))
      return;
 
    DynList<size_t> zeros;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        p.second = p.second / s;
        if (coeff_is_zero(p.second))
          zeros.append(p.first);
      }
 
    zeros.for_each([this](size_t exp)
    {
      coeffs.remove_key(exp);
    });
  }
 
  template <class Op>
  void for_each_term_desc(Op &&op) const
  {
    Array<const Term *> desc_terms(coeffs.size());
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        const auto &p = it.get_curr();
        desc_terms.append(&p);
      }
 
    for (size_t i = desc_terms.size(); i > 0; --i)
      {
        const Term &term = *desc_terms(i - 1);
        op(term.first, term.second);
      }
  }
 
  [[nodiscard]] Gen_Polynomial multiply_by_monic_linear(const Coefficient &c) const
  {
    if (is_zero())
      return Gen_Polynomial();
 
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result.add_to_coeff(p.first, p.second * c);
        result.add_to_coeff(p.first + 1, p.second);
      }
    return result;
  }
 
public:
  // =================================================================
  //  Construction
  // =================================================================
 
  Gen_Polynomial() = default;
 
  explicit Gen_Polynomial(const Coefficient &c)
  {
    if (not coeff_is_zero(c))
      coeffs[0] = c;
  }
 
  Gen_Polynomial(const Coefficient &c, size_t exp)
  {
    if (not coeff_is_zero(c))
      coeffs[exp] = c;
  }
 
  Gen_Polynomial(std::initializer_list<Coefficient> il)
  {
    size_t exp = 0;
    for (auto &c : il)
      {
        if (not coeff_is_zero(c))
          coeffs[exp] = c;
        ++exp;
      }
  }
 
  explicit Gen_Polynomial(const DynList<Coefficient> &l)
  {
    size_t exp = 0;
    l.for_each([this, &exp](const Coefficient &c)
    {
      if (not coeff_is_zero(c))
        coeffs[exp] = c;
      ++exp;
    });
  }
 
  explicit Gen_Polynomial(const DynList<std::pair<size_t, Coefficient>> &term_list)
  {
    term_list.for_each([this](const Term &t)
    {
      if (not coeff_is_zero(t.second))
        add_to_coeff(t.first, t.second);
    });
  }
 
  // =================================================================
  //  Static factory methods
  // =================================================================
 
  [[nodiscard]] static Gen_Polynomial zero()
  {
    return Gen_Polynomial();
  }
 
  [[nodiscard]] static Gen_Polynomial one()
  {
    return Gen_Polynomial(Coefficient(1));
  }
 
  [[nodiscard]] static Gen_Polynomial x_to(size_t n)
  {
    return Gen_Polynomial(Coefficient(1), n);
  }
 
  [[nodiscard]] static Gen_Polynomial from_roots(const DynList<Coefficient> &roots)
  {
    Gen_Polynomial result(Coefficient(1));  // start with constant 1
    roots.for_each([&result](const Coefficient &r)
    {
      result = result.multiply_by_monic_linear(-r);  // (x - r)
    });
    return result;
  }
 
  [[nodiscard]] static Gen_Polynomial interpolate(
    const DynList<std::pair<Coefficient, Coefficient>> &points)
  {
    ah_domain_error_if(points.is_empty()) << "interpolate requires at least one point";
 
    const size_t n = points.size();
    Array<Coefficient> xs;
    Array<Coefficient> dd;
    xs.putn(n);
    dd.putn(n);
 
    size_t idx = 0;
    points.for_each([&](const std::pair<Coefficient, Coefficient> &pt)
    {
      xs[idx] = pt.first;
      dd[idx] = pt.second;
      ++idx;
    });
 
    for (size_t i = 0; i < n; ++i)
      for (size_t j = i + 1; j < n; ++j)
        ah_domain_error_if(coeff_is_zero(xs[i] - xs[j])) << "interpolate: duplicate x-values";
 
    for (size_t order = 1; order < n; ++order)
      for (size_t i = n - 1; i >= order; --i)
        {
          Coefficient denom = xs[i] - xs[i - order];
          dd[i] = (dd[i] - dd[i - 1]) / denom;
          if (i == order)
            break;
        }
 
    Gen_Polynomial result(dd[0]);
    Gen_Polynomial basis(Coefficient(1));
    for (size_t i = 1; i < n; ++i)
      {
        basis = basis.multiply_by_monic_linear(-xs[i - 1]);
        if (not coeff_is_zero(dd[i]))
          result.add_scaled_shifted(basis, dd[i]);
      }
 
    return result;
  }
 
  // =================================================================
  //  Properties
  // =================================================================
 
  [[nodiscard]] bool is_zero() const noexcept
  {
    return coeffs.is_empty();
  }
 
  [[nodiscard]] size_t degree() const noexcept
  {
    if (coeffs.is_empty())
      return 0;
    return coeffs.max().first;
  }
 
  [[nodiscard]] size_t num_terms() const noexcept
  {
    return coeffs.size();
  }
 
  [[nodiscard]] bool is_constant() const noexcept
  {
    return coeffs.is_empty() or (coeffs.size() == 1 and coeffs.min().first == 0);
  }
 
  [[nodiscard]] bool is_monomial() const noexcept
  {
    return coeffs.size() == 1;
  }
 
  [[nodiscard]] bool is_monic() const noexcept
  {
    return not coeffs.is_empty() and coeff_is_zero(coeffs.max().second - Coefficient(1));
  }
 
  // =================================================================
  //  Coefficient access
  // =================================================================
 
  [[nodiscard]] Coefficient operator[](size_t exp) const
  {
    return coeff_at(exp);
  }
 
  [[nodiscard]] Coefficient leading_coeff() const noexcept
  {
    if (coeffs.is_empty())
      return Coefficient{};
    return coeffs.max().second;
  }
 
  [[nodiscard]] bool has_term(size_t exp) const noexcept
  {
    return coeffs.has(exp);
  }
 
  // =================================================================
  //  Evaluation
  // =================================================================
 
  [[nodiscard]] Coefficient horner_eval(const Coefficient &x) const
  {
    if (coeffs.is_empty())
      return Coefficient{};
 
    if (x == Coefficient{})
      return coeff_at(0);
 
    Coefficient result = Coefficient{};
    size_t prev_exp = 0;
    bool first = true;
 
    for_each_term_desc([&](size_t exp, const Coefficient &c)
    {
      if (first)
        {
          result = c;
          prev_exp = exp;
          first = false;
          return;
        }
 
      result = result * polynomial_detail::power(x, prev_exp - exp) + c;
      prev_exp = exp;
    });
 
    if (prev_exp > 0)
      result = result * polynomial_detail::power(x, prev_exp);
 
    return result;
  }
 
  [[nodiscard]] Coefficient sparse_eval(const Coefficient &x) const
  {
    Coefficient result = Coefficient{};
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result = result + p.second * polynomial_detail::power(x, p.first);
      }
    return result;
  }
 
  [[nodiscard]] Coefficient eval(const Coefficient &x) const
  {
    if (coeffs.is_empty())
      return Coefficient{};
 
    if (x == Coefficient{})
      return coeff_at(0);
 
    if (num_terms() <= 2)
      return sparse_eval(x);
 
    return horner_eval(x);
  }
 
  [[nodiscard]] Coefficient operator()(const Coefficient &x) const
  {
    return eval(x);
  }
 
  [[nodiscard]] DynList<Coefficient> multi_eval(const DynList<Coefficient> &points) const
  {
    DynList<Coefficient> results;
    points.for_each([this, &results](const Coefficient &x)
    {
      results.append(eval(x));
    });
    return results;
  }
 
  // =================================================================
  //  Arithmetic operators
  // =================================================================
 
  [[nodiscard]] Gen_Polynomial operator-() const
  {
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result.coeffs[p.first] = -p.second;
      }
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial operator+(const Gen_Polynomial &q) const
  {
    Gen_Polynomial result = *this;
    result += q;
    return result;
  }
 
  Gen_Polynomial &operator+=(const Gen_Polynomial &q)
  {
    if (this == &q)
      {
        scale_inplace(Coefficient(1) + Coefficient(1));
        return *this;
      }
 
    for (auto it = q.coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        add_to_coeff(p.first, p.second);
      }
    return *this;
  }
 
  [[nodiscard]] Gen_Polynomial operator+(const Coefficient &c) const
  {
    Gen_Polynomial result = *this;
    result += c;
    return result;
  }
 
  Gen_Polynomial &operator+=(const Coefficient &c)
  {
    add_to_coeff(0, c);
    return *this;
  }
 
  [[nodiscard]] Gen_Polynomial operator-(const Gen_Polynomial &q) const
  {
    Gen_Polynomial result = *this;
    result -= q;
    return result;
  }
 
  Gen_Polynomial &operator-=(const Gen_Polynomial &q)
  {
    if (this == &q)
      {
        coeffs = DynMapTree<size_t, Coefficient>();
        return *this;
      }
 
    for (auto it = q.coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        add_to_coeff(p.first, -p.second);
      }
    return *this;
  }
 
  [[nodiscard]] Gen_Polynomial operator-(const Coefficient &c) const
  {
    Gen_Polynomial result = *this;
    result -= c;
    return result;
  }
 
  Gen_Polynomial &operator-=(const Coefficient &c)
  {
    add_to_coeff(0, -c);
    return *this;
  }
 
  [[nodiscard]] Gen_Polynomial operator*(const Gen_Polynomial &q) const
  {
    if (is_zero() or q.is_zero())
      return Gen_Polynomial();
 
    Gen_Polynomial result;
    const Gen_Polynomial *outer = this;
    const Gen_Polynomial *inner = &q;
 
    if (q.num_terms() < num_terms())
      {
        outer = &q;
        inner = this;
      }
 
    for (auto it1 = outer->coeffs.get_it(); it1.has_curr(); it1.next_ne())
      {
        auto &t1 = it1.get_curr();
        for (auto it2 = inner->coeffs.get_it(); it2.has_curr(); it2.next_ne())
          {
            auto &t2 = it2.get_curr();
            result.add_to_coeff(t1.first + t2.first, t1.second * t2.second);
          }
      }
    return result;
  }
 
  Gen_Polynomial &operator*=(const Gen_Polynomial &q)
  {
    *this = *this * q;
    return *this;
  }
 
  [[nodiscard]] Gen_Polynomial operator*(const Coefficient &s) const
  {
    if (coeff_is_zero(s))
      return Gen_Polynomial();
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        Coefficient prod = p.second * s;
        if (not coeff_is_zero(prod))
          result.coeffs[p.first] = prod;
      }
    return result;
  }
 
  Gen_Polynomial &operator*=(const Coefficient &s)
  {
    scale_inplace(s);
    return *this;
  }
 
  [[nodiscard]] Gen_Polynomial operator/(const Coefficient &s) const
  {
    ah_domain_error_if(coeff_is_zero(s)) << "polynomial division by zero scalar";
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        Coefficient q = p.second / s;
        if (not coeff_is_zero(q))
          result.coeffs[p.first] = q;
      }
    return result;
  }
 
  Gen_Polynomial &operator/=(const Coefficient &s)
  {
    divide_scalar_inplace(s);
    return *this;
  }
 
  // =================================================================
  //  Polynomial division
  // =================================================================
 
  [[nodiscard]] std::pair<Gen_Polynomial, Gen_Polynomial> divmod(const Gen_Polynomial &d) const
  {
    ah_domain_error_if(d.is_zero()) << "polynomial division by zero";
 
    if (is_zero() or degree() < d.degree())
      return {Gen_Polynomial(), *this};
 
    Gen_Polynomial q;          // quotient
    Gen_Polynomial r = *this;  // remainder
    const size_t d_degree = d.degree();
    const Coefficient d_lc = d.leading_coeff();
 
    while (not r.is_zero())
      {
        const size_t r_degree = r.degree();
        if (r_degree < d_degree)
          break;
 
        const Coefficient scale = r.leading_coeff() / d_lc;
        ah_domain_error_if(coeff_is_zero(scale))
          << "divmod: inexact coefficient division (leading_coeff=" << r.leading_coeff()
          << " / divisor_lc=" << d_lc << " = 0). "
          << "For integral/non-field coefficients, use pseudo_divmod() instead.";
        const size_t shift = r_degree - d_degree;
 
        q.add_to_coeff(shift, scale);
        r.add_scaled_shifted(d, -scale, shift);
      }
 
    return {q, r};
  }
 
  [[nodiscard]] Gen_Polynomial operator/(const Gen_Polynomial &d) const
  {
    return divmod(d).first;
  }
 
  [[nodiscard]] Gen_Polynomial operator%(const Gen_Polynomial &d) const
  {
    return divmod(d).second;
  }
 
  Gen_Polynomial &operator/=(const Gen_Polynomial &d)
  {
    *this = divmod(d).first;
    return *this;
  }
 
  Gen_Polynomial &operator%=(const Gen_Polynomial &d)
  {
    *this = divmod(d).second;
    return *this;
  }
 
  // =================================================================
  //  Calculus
  // =================================================================
 
  [[nodiscard]] Gen_Polynomial derivative() const
  {
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        if (p.first == 0)
          continue;  // derivative of constant is 0
        Coefficient dc = p.second * Coefficient(p.first);
        if (not coeff_is_zero(dc))
          result.coeffs[p.first - 1] = dc;
      }
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial nth_derivative(size_t n) const
  {
    Gen_Polynomial result = *this;
    for (size_t i = 0; i < n and not result.is_zero(); ++i)
      result = result.derivative();
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial integral(const Coefficient &C = Coefficient{}) const
  {
    Gen_Polynomial result;
    if (not coeff_is_zero(C))
      result.coeffs[0] = C;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        Coefficient ic = p.second / Coefficient(p.first + 1);
        if (not coeff_is_zero(ic))
          result.coeffs[p.first + 1] = ic;
      }
    return result;
  }
 
  [[nodiscard]] Coefficient definite_integral(const Coefficient &a, const Coefficient &b) const
  {
    Coefficient result = Coefficient{};
 
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        const size_t next_exp = p.first + 1;
        const Coefficient scale = p.second / Coefficient(next_exp);
        if (coeff_is_zero(scale))
          continue;
 
        result
          = result
          + scale * (polynomial_detail::power(b, next_exp) - polynomial_detail::power(a, next_exp));
      }
 
    return result;
  }
 
  // =================================================================
  //  Composition and power
  // =================================================================
 
  [[nodiscard]] Gen_Polynomial compose(const Gen_Polynomial &q) const
  {
    if (is_zero())
      return Gen_Polynomial();
 
    if (q.is_constant())
      return Gen_Polynomial(eval(q[0]));
 
    Gen_Polynomial result;
    size_t prev_exp = 0;
    bool first = true;
 
    for_each_term_desc([&](size_t exp, const Coefficient &c)
    {
      if (first)
        {
          result = Gen_Polynomial(c);
          prev_exp = exp;
          first = false;
          return;
        }
 
      const size_t gap = prev_exp - exp;
      if (gap == 1)
        result *= q;
      else if (gap > 1)
        result *= q.pow(gap);
 
      result.add_to_coeff(0, c);
      prev_exp = exp;
    });
 
    if (prev_exp == 1)
      result *= q;
    else if (prev_exp > 1)
      result *= q.pow(prev_exp);
 
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial pow(size_t n) const
  {
    Gen_Polynomial result(Coefficient(1));
    Gen_Polynomial base = *this;
    while (n > 0)
      {
        if (n & 1)
          result *= base;
        base *= base;
        n >>= 1;
      }
    return result;
  }
 
  // =================================================================
  //  GCD
  // =================================================================
 
  [[nodiscard]] static Gen_Polynomial gcd(Gen_Polynomial a, Gen_Polynomial b)
  {
    while (not b.is_zero())
      {
        Gen_Polynomial r = a % b;
        a = b;
        b = r;
      }
    // Make monic
    if (not a.is_zero())
      a /= a.leading_coeff();
    return a;
  }
 
  struct Xgcd_Result
  {
    Gen_Polynomial g;  
    Gen_Polynomial s;  
    Gen_Polynomial t;  
  };
 
  [[nodiscard]] static Xgcd_Result xgcd(Gen_Polynomial a, Gen_Polynomial b)
  {
    Gen_Polynomial s0(Coefficient(1)), s1;
    Gen_Polynomial t0, t1(Coefficient(1));
 
    while (not b.is_zero())
      {
        auto [q, r] = a.divmod(b);
        a = b;
        b = r;
        Gen_Polynomial s2 = s0 - q * s1;
        Gen_Polynomial t2 = t0 - q * t1;
        s0 = s1;
        s1 = s2;
        t0 = t1;
        t1 = t2;
      }
 
    if (not a.is_zero())
      {
        Coefficient lc = a.leading_coeff();
        a /= lc;
        s0 /= lc;
        t0 /= lc;
      }
 
    return {a, s0, t0};
  }
 
  [[nodiscard]] static Gen_Polynomial lcm(const Gen_Polynomial &a, const Gen_Polynomial &b)
  {
    if (a.is_zero() or b.is_zero())
      return Gen_Polynomial();
    Gen_Polynomial g = gcd(a, b);
    return a / g * b;
  }
 
  [[nodiscard]] std::pair<Gen_Polynomial, Gen_Polynomial> pseudo_divmod(const Gen_Polynomial &d) const
  {
    ah_domain_error_if(d.is_zero()) << "pseudo-division by zero polynomial";
 
    if (is_zero() or degree() < d.degree())
      return {Gen_Polynomial(), *this};
 
    const Coefficient lc_d = d.leading_coeff();
    const size_t d_degree = d.degree();
    size_t delta = degree() - d_degree + 1;
 
    Gen_Polynomial q;
    Gen_Polynomial r = *this;
 
    while (not r.is_zero())
      {
        const size_t r_degree = r.degree();
        if (r_degree < d_degree)
          break;
 
        const Coefficient lc_r = r.leading_coeff();
        const size_t exp = r_degree - d_degree;
 
        // Multiply q and r by lc(d) to avoid division
        q.scale_inplace(lc_d);
        r.scale_inplace(lc_d);
        --delta;
 
        q.add_to_coeff(exp, lc_r);
        r.add_scaled_shifted(d, -lc_r, exp);
      }
 
    // Multiply by remaining powers of lc(d)
    if (delta > 0)
      {
        const Coefficient scale = polynomial_detail::power(lc_d, delta);
        q.scale_inplace(scale);
        r.scale_inplace(scale);
      }
 
    return {q, r};
  }
 
  // =================================================================
  //  Algebraic transformations
  // =================================================================
 
  [[nodiscard]] Gen_Polynomial to_monic() const
  {
    ah_domain_error_if(is_zero()) << "cannot make zero polynomial monic";
    return *this / leading_coeff();
  }
 
  [[nodiscard]] Gen_Polynomial reverse() const
  {
    if (is_zero())
      return Gen_Polynomial();
    size_t d = degree();
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result.coeffs[d - p.first] = p.second;
      }
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial negate_arg() const
  {
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result.coeffs[p.first] = (p.first % 2 == 0) ? p.second : -p.second;
      }
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial shift(const Coefficient &k) const
  {
    return compose(Gen_Polynomial({k, Coefficient(1)}));
  }
 
  [[nodiscard]] Gen_Polynomial shift_up(size_t k) const
  {
    if (k == 0)
      return *this;
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result.coeffs[p.first + k] = p.second;
      }
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial shift_down(size_t k) const
  {
    if (k == 0)
      return *this;
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        if (p.first >= k)
          result.coeffs[p.first - k] = p.second;
      }
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial truncate(size_t n) const
  {
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        if (p.first < n)
          result.coeffs[p.first] = p.second;
      }
    return result;
  }
 
  [[nodiscard]] Array<Coefficient> to_dense() const
  {
    if (is_zero())
      return Array<Coefficient>();
    size_t d = degree();
    Array<Coefficient> result(d + 1, Coefficient{});
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result(p.first) = p.second;
      }
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial square_free() const
  {
    if (is_zero() or is_constant())
      return *this;
    Gen_Polynomial g = gcd(*this, derivative());
    return *this / g;
  }
 
  // =================================================================
  //  Root analysis (floating-point coefficients)
  // =================================================================
 
  [[nodiscard]] Coefficient cauchy_bound() const
  {
    ah_domain_error_if(is_zero()) << "cauchy_bound: zero polynomial has no roots";
    if (is_constant())
      return Coefficient{};
 
    const size_t d = degree();
    const Coefficient lc = leading_coeff();
    if constexpr (std::is_floating_point_v<Coefficient>)
      {
        Coefficient max_ratio = Coefficient{};
        for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
          {
            auto &p = it.get_curr();
            if (p.first == d)
              continue;
            const Coefficient ratio = std::abs(p.second / lc);
            if (ratio > max_ratio)
              max_ratio = ratio;
          }
        return Coefficient(1) + max_ratio;
      }
    else
      {
        long double max_ratio_ld = 0.0L;
        for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
          {
            auto &p = it.get_curr();
            if (p.first == d)
              continue;
            const long double ratio
              = std::fabs(static_cast<long double>(p.second) / static_cast<long double>(lc));
            if (ratio > max_ratio_ld)
              max_ratio_ld = ratio;
          }
        return Coefficient(1) + static_cast<Coefficient>(std::ceil(max_ratio_ld));
      }
  }
 
  [[nodiscard]] size_t sign_variations() const
  {
    if (coeffs.size() < 2)
      return 0;
 
    size_t changes = 0;
    bool have_prev = false;
    bool prev_positive = false;
 
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        bool positive;
        if constexpr (std::is_floating_point_v<Coefficient>)
          positive = p.second > epsilon();
        else
          positive = p.second > Coefficient{};
 
        bool is_neg;
        if constexpr (std::is_floating_point_v<Coefficient>)
          is_neg = p.second < -epsilon();
        else
          is_neg = p.second < Coefficient{};
 
        if (not positive and not is_neg)
          continue;  // skip zero coefficients
 
        if (have_prev and positive != prev_positive)
          ++changes;
        prev_positive = positive;
        have_prev = true;
      }
    return changes;
  }
 
  [[nodiscard]] DynList<Gen_Polynomial> sturm_chain() const
  {
    DynList<Gen_Polynomial> chain;
    chain.append(*this);
    if (is_zero() or is_constant())
      return chain;
 
    Gen_Polynomial p0 = *this;
    Gen_Polynomial p1 = derivative();
    chain.append(p1);
 
    while (not p1.is_zero())
      {
        Gen_Polynomial r = -(p0 % p1);
        if (r.is_zero())
          break;
        chain.append(r);
        p0 = p1;
        p1 = r;
      }
 
    return chain;
  }
 
  [[nodiscard]] static size_t sturm_sign_changes(const DynList<Gen_Polynomial> &chain,
                                                 const Coefficient &x)
  {
    size_t changes = 0;
    bool have_prev = false;
    bool prev_positive = false;
 
    chain.for_each([&](const Gen_Polynomial &p)
    {
      Coefficient val = p(x);
      if (coeff_is_zero(val))
        return;  // skip zeros in Sturm evaluation
 
      bool positive = val > Coefficient{};
      if (have_prev and positive != prev_positive)
        ++changes;
      prev_positive = positive;
      have_prev = true;
    });
 
    return changes;
  }
 
  [[nodiscard]] size_t count_real_roots(const Coefficient &a, const Coefficient &b) const
  {
    if (b < a)
      return count_real_roots(b, a);
 
    if (is_zero())
      return 0;
 
    auto chain = sturm_chain();
    size_t sa = sturm_sign_changes(chain, a);
    size_t sb = sturm_sign_changes(chain, b);
    size_t count = sa > sb ? sa - sb : 0;
 
    // Sturm's theorem counts roots in (a, b]; add root at a if present
    if (coeff_is_zero((*this)(a)))
      ++count;
 
    return count;
  }
 
  [[nodiscard]] size_t count_all_real_roots() const
  {
    if (is_zero() or is_constant())
      return 0;
    Coefficient bound = cauchy_bound();
    return count_real_roots(-bound, bound);
  }
 
  [[nodiscard]] Coefficient bisect_root(Coefficient a,
                                        Coefficient b,
                                        Coefficient tol = Coefficient(1e-12),
                                        size_t max_iter = 200) const
  {
    if (b < a)
      std::swap(a, b);
 
    Coefficient fa = eval(a), fb = eval(b);
    if (coeff_is_zero(fa))
      return a;
    if (coeff_is_zero(fb))
      return b;
 
    ah_domain_error_if((fa > Coefficient{}) == (fb > Coefficient{}) and not coeff_is_zero(fa)
                       and not coeff_is_zero(fb))
      << "bisect_root: f(a) and f(b) must have opposite signs";
 
    for (size_t i = 0; i < max_iter; ++i)
      {
        Coefficient mid = (a + b) / Coefficient(2);
        if (b - a < tol)
          return mid;
        Coefficient fm = eval(mid);
        if (coeff_is_zero(fm))
          return mid;
        if (fa > Coefficient{} == fm > Coefficient{})
          {
            a = mid;
            fa = fm;
          }
        else
          {
            b = mid;
            fb = fm;
          }
      }
    return (a + b) / Coefficient(2);
  }
 
  [[nodiscard]] Coefficient newton_root(Coefficient x0,
                                        Coefficient tol = Coefficient(1e-12),
                                        size_t max_iter = 100) const
    requires std::is_floating_point_v<Coefficient>
  {
    Gen_Polynomial dp = derivative();
 
    for (size_t i = 0; i < max_iter; ++i)
      {
        Coefficient fx = eval(x0);
        if (coeff_is_zero(fx))
          return x0;
        Coefficient dfx = dp.eval(x0);
        ah_domain_error_if(coeff_is_zero(dfx)) << "newton_root: derivative is zero at x = " << x0;
        Coefficient x1 = x0 - fx / dfx;
        Coefficient delta = x1 - x0;
        if constexpr (std::is_floating_point_v<Coefficient>)
          {
            if (std::abs(delta) < tol)
              return x1;
          }
        else
          {
            if (delta == Coefficient{})
              return x1;
          }
        x0 = x1;
      }
    return x0;
  }
 
  // =================================================================
  //  Comparison
  // =================================================================
 
  [[nodiscard]] bool operator==(const Gen_Polynomial &q) const noexcept
  {
    if (this == &q)
      return true;
 
    auto it1 = coeffs.get_it();
    auto it2 = q.coeffs.get_it();
 
    while (it1.has_curr() and it2.has_curr())
      {
        auto &lhs = it1.get_curr();
        auto &rhs = it2.get_curr();
 
        if (lhs.first == rhs.first)
          {
            if (not coeff_is_zero(lhs.second - rhs.second))
              return false;
            it1.next_ne();
            it2.next_ne();
            continue;
          }
 
        if (lhs.first < rhs.first)
          {
            if (not coeff_is_zero(lhs.second))
              return false;
            it1.next_ne();
            continue;
          }
 
        if (not coeff_is_zero(rhs.second))
          return false;
        it2.next_ne();
      }
 
    while (it1.has_curr())
      {
        if (not coeff_is_zero(it1.get_curr().second))
          return false;
        it1.next_ne();
      }
 
    while (it2.has_curr())
      {
        if (not coeff_is_zero(it2.get_curr().second))
          return false;
        it2.next_ne();
      }
 
    return true;
  }
 
  [[nodiscard]] bool operator!=(const Gen_Polynomial &q) const noexcept
  {
    return not(*this == q);
  }
 
  // =================================================================
  //  Iteration
  // =================================================================
 
  template <class Op>
  void for_each_term(Op &&op) const
  {
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        op(p.first, p.second);
      }
  }
 
  [[nodiscard]] DynList<std::pair<size_t, Coefficient>> terms_list() const
  {
    DynList<std::pair<size_t, Coefficient>> result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result.append({p.first, p.second});
      }
    return result;
  }
 
  [[nodiscard]] DynList<size_t> exponents() const
  {
    return coeffs.keys();
  }
 
  [[nodiscard]] DynList<Coefficient> coefficients() const
  {
    return coeffs.values();
  }
 
  // =================================================================
  //  String representation
  // =================================================================
 
  [[nodiscard]] std::string to_str(const std::string &var = "x") const
  {
    if (coeffs.is_empty())
      return "0";
 
    std::ostringstream oss;
    bool first = true;
 
    for_each_term_desc([&](size_t exp, const Coefficient &c)
    {
      bool negative = false;
      if constexpr (requires { c < Coefficient{}; })
        negative = c < Coefficient{};
 
      Coefficient ac = negative ? -c : c;
 
      // Sign
      if (first)
        {
          if (negative)
            oss << "-";
        }
      else
        oss << (negative ? " - " : " + ");
 
      // Coefficient and variable
      if (exp == 0)
        oss << ac;
      else
        {
          if (not(ac == Coefficient(1)))
            oss << ac;
          oss << var;
          if (exp > 1)
            oss << "^" << exp;
        }
 
      first = false;
    });
 
    return oss.str();
  }
 
  // =================================================================
  //  Layer 5 — Exact Univariate Factorization
  // =================================================================
 
  struct SfdTerm
  {
    Gen_Polynomial factor;
    size_t multiplicity;
  };
 
  [[nodiscard]] Coefficient get_coeff(size_t exp) const noexcept
  {
    return coeff_at(exp);
  }
 
  void set_coeff(size_t exp, const Coefficient &c)
  {
    if (auto *slot = coeff_ptr(exp); slot != nullptr)
      {
        if (coeff_is_zero(c))
          coeffs.remove_key(exp);
        else
          *slot = c;
      }
    else if (not coeff_is_zero(c))
      coeffs.insert(exp, c);
  }
 
  [[nodiscard]] Coefficient content() const
    requires std::is_integral_v<Coefficient>
  {
    if (coeffs.is_empty())
      return Coefficient{};
 
    Coefficient result = Coefficient{};
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result = std::gcd(result, p.second);
      }
 
    // Ensure sign follows leading coefficient
    const Coefficient lc = leading_coeff();
    if (lc < Coefficient{} and result > Coefficient{})
      result = -result;
 
    return result;
  }
 
  [[nodiscard]] Gen_Polynomial primitive_part() const
    requires std::is_integral_v<Coefficient>
  {
    if (is_zero())
      return Gen_Polynomial();
 
    Coefficient c = content();
    if (c == Coefficient{} or c == Coefficient(1) or c == Coefficient(-1))
      {
        // Already primitive; ensure leading coefficient is positive
        Gen_Polynomial result = *this;
        if (leading_coeff() < Coefficient{})
          result.scale_inplace(Coefficient(-1));
        return result;
      }
 
    Gen_Polynomial result;
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        result.coeffs.insert(p.first, p.second / c);
      }
 
    // Ensure leading coefficient is positive
    if (result.leading_coeff() < Coefficient{})
      result.scale_inplace(Coefficient(-1));
 
    return result;
  }
 
  static Gen_Polynomial integer_gcd(Gen_Polynomial a, Gen_Polynomial b)
    requires std::is_integral_v<Coefficient>
  {
    if (a.is_zero())
      return b.primitive_part();
    if (b.is_zero())
      return a.primitive_part();
 
    a = a.primitive_part();
    b = b.primitive_part();
 
    while (not b.is_zero())
      {
        auto [quo, rem] = a.pseudo_divmod(b);
        a = b;
        b = rem.is_zero() ? rem : rem.primitive_part();
      }
 
    return a.primitive_part();
  }
 
  [[nodiscard]] DynList<SfdTerm> yun_sfd() const
    requires std::is_integral_v<Coefficient>
  {
    DynList<SfdTerm> result;
 
    if (is_constant())
      return result;
 
    Gen_Polynomial f = primitive_part();
    Gen_Polynomial fp = f.derivative();
 
    Gen_Polynomial g = integer_gcd(f, fp);
    Gen_Polynomial h = _integer_exact_quot(f, g);
 
    size_t multiplicity = 1;
 
    while (not h.is_constant())
      {
        Gen_Polynomial h_prev = h;
        h = integer_gcd(g, h);
        Gen_Polynomial w = _integer_exact_quot(h_prev, h);
 
        if (not w.is_constant())
          result.append(SfdTerm{w.primitive_part(), multiplicity});
 
        g = _integer_exact_quot(g, h);
        multiplicity++;
      }
 
    // Final factor
    if (not g.is_constant())
      result.append(SfdTerm{g.primitive_part(), multiplicity});
 
    return result;
  }
 
  static Gen_Polynomial _integer_exact_quot(const Gen_Polynomial &a, const Gen_Polynomial &b)
    requires std::is_integral_v<Coefficient>
  {
    ah_domain_error_if(b.is_zero()) << "_integer_exact_quot: zero divisor";
    if (b.is_constant())
      {
        Coefficient d = b.leading_coeff();
        Gen_Polynomial r;
        for (auto it = a.coeffs.get_it(); it.has_curr(); it.next_ne())
          {
            auto &p = it.get_curr();
            ah_domain_error_if(not coeff_is_zero(p.second % d))
              << "_integer_exact_quot: inexact constant division for coefficient "
              << p.second << " and divisor " << d;
            r.coeffs.insert(p.first, p.second / d);
          }
        return r;
      }
    auto [q, rem] = a.pseudo_divmod(b);
    ah_domain_error_if(not rem.is_zero())
      << "_integer_exact_quot: divisor does not divide dividend exactly";
    return q.is_zero() ? q : q.primitive_part();
  }
 
  [[nodiscard]] static Gen_Polynomial reduce_coeffs_mod(const Gen_Polynomial &f, Coefficient mod)
    requires std::is_integral_v<Coefficient>
  {
    Gen_Polynomial reduced;
    f.for_each_term([&reduced, mod](size_t exp, const Coefficient &coeff)
    {
      Coefficient c = coeff % mod;
      if (c < Coefficient{})
        c = c + mod;
      if (c != Coefficient{})
        reduced.set_coeff(exp, c);
    });
    return reduced;
  }
 
  [[nodiscard]] static uint64_t eval_mod_u64(const Gen_Polynomial &f, uint64_t x, uint64_t mod)
    requires std::is_integral_v<Coefficient>
  {
    if (mod == 0 or mod == 1)
      return 0;
 
    uint64_t result = 0;
    f.for_each_term([&result, x, mod](size_t exp, const Coefficient &coeff)
    {
      const uint64_t c = polynomial_detail::normalize_mod_u64(coeff, mod);
      const uint64_t xpow = mod_exp(x, exp, mod);
      const uint64_t term = mod_mul(c, xpow, mod);
      result = (result + term) % mod;
    });
    return result;
  }
 
  [[nodiscard]] static Gen_Polynomial divide_by_monic_linear_mod(const Gen_Polynomial &f,
                                                                 Coefficient root,
                                                                 Coefficient mod)
    requires std::is_integral_v<Coefficient>
  {
    ah_domain_error_if(f.degree() == 0)
      << "divide_by_monic_linear_mod: constant polynomial is not divisible by x - r";
 
    const size_t degree = f.degree();
    Array<Coefficient> synthetic(degree + 1, Coefficient{});
 
    synthetic(degree) = ((f.get_coeff(degree) % mod) + mod) % mod;
 
    for (size_t idx = degree; idx-- > 0;)
      {
        Coefficient term = (f.get_coeff(idx) % mod + mod) % mod;
        Coefficient next = (term + (synthetic(idx + 1) * root) % mod) % mod;
        synthetic(idx) = next;
      }
 
    ah_domain_error_if(synthetic(0) != Coefficient{})
      << "divide_by_monic_linear_mod: x - r does not divide polynomial modulo p";
 
    Gen_Polynomial quotient;
    for (size_t exp = 0; exp < degree; ++exp)
      if (synthetic(exp + 1) != Coefficient{})
        quotient.set_coeff(exp, synthetic(exp + 1));
 
    return quotient;
  }
 
  [[nodiscard]] static DynList<Coefficient> positive_divisors(Coefficient value)
    requires std::is_integral_v<Coefficient>
  {
    DynList<Coefficient> result;
 
    const uint64_t abs_value = polynomial_detail::abs_to_u64(value);
    if (abs_value == 0)
      return result;
 
    for (uint64_t d = 1; d * d <= abs_value; ++d)
      if (abs_value % d == 0)
        {
          result.append(static_cast<Coefficient>(d));
          const uint64_t other = abs_value / d;
          if (other != d)
            result.append(static_cast<Coefficient>(other));
        }
 
    return result;
  }
 
  [[nodiscard]] static bool try_exact_candidate_division(const Gen_Polynomial &f,
                                                         const Gen_Polynomial &candidate,
                                                         Gen_Polynomial &quotient)
    requires std::is_integral_v<Coefficient>
  {
    try
      {
        auto [q, r] = f.divmod(candidate);
        if (not r.is_zero())
          return false;
        quotient = std::move(q);
        return true;
      }
    catch (const std::domain_error &)
      {
        return false;
      }
  }
 
  [[nodiscard]] static bool try_divide_by_linear_factor(const Gen_Polynomial &f,
                                                        Coefficient a,
                                                        Coefficient b,
                                                        Gen_Polynomial &quotient)
    requires std::is_integral_v<Coefficient>
  {
    if (f.is_zero() or f.degree() == 0)
      return false;
 
    if (a == Coefficient{})
      return false;
 
    const size_t degree = f.degree();
    Array<Coefficient> q_coeffs(degree, Coefficient{});
 
    const Coefficient leading = f.get_coeff(degree);
    if (leading % a != Coefficient{})
      return false;
 
    q_coeffs(degree - 1) = leading / a;
 
    for (size_t k = degree - 1; k > 0; --k)
      {
        const Coefficient numerator = f.get_coeff(k) - b * q_coeffs(k);
        if (numerator % a != Coefficient{})
          return false;
        q_coeffs(k - 1) = numerator / a;
      }
 
    if (f.get_coeff(0) - b * q_coeffs(0) != Coefficient{})
      return false;
 
    quotient = Gen_Polynomial();
    for (size_t exp = 0; exp < degree; ++exp)
      if (q_coeffs(exp) != Coefficient{})
        quotient.set_coeff(exp, q_coeffs(exp));
 
    return true;
  }
 
  [[nodiscard]] static bool try_extract_rational_linear_factor(const Gen_Polynomial &f,
                                                               Gen_Polynomial &factor,
                                                               Gen_Polynomial &quotient)
    requires std::is_integral_v<Coefficient>
  {
    if (f.is_zero() or f.degree() == 0)
      return false;
 
    if (f.get_coeff(0) == Coefficient{})
      {
        if (try_divide_by_linear_factor(f, Coefficient(1), Coefficient{}, quotient))
          {
            factor = Gen_Polynomial(Coefficient(1), 1);
            return true;
          }
        return false;
      }
 
    DynList<Coefficient> numerator_divs = positive_divisors(f.get_coeff(0));
    DynList<Coefficient> denominator_divs = positive_divisors(f.leading_coeff());
 
    for (auto den_it = denominator_divs.get_it(); den_it.has_curr(); den_it.next_ne())
      {
        const Coefficient den = den_it.get_curr();
 
        for (auto num_it = numerator_divs.get_it(); num_it.has_curr(); num_it.next_ne())
          {
            const Coefficient num_abs = num_it.get_curr();
 
            for (Coefficient sign : {Coefficient(1), Coefficient(-1)})
              {
                const Coefficient num = sign * num_abs;
                if (std::gcd(polynomial_detail::abs_to_u64(num), polynomial_detail::abs_to_u64(den))
                    != 1)
                  continue;
 
                if (try_divide_by_linear_factor(f, den, -num, quotient))
                  {
                    factor = Gen_Polynomial();
                    factor.set_coeff(0, -num);
                    factor.set_coeff(1, den);
                    return true;
                  }
              }
          }
      }
 
    return false;
  }
 
  [[nodiscard]] static bool try_extract_primitive_quadratic_factor(const Gen_Polynomial &f,
                                                                   Gen_Polynomial &factor,
                                                                   Gen_Polynomial &quotient)
    requires std::is_integral_v<Coefficient>
  {
    if (f.is_zero() or f.degree() < 2)
      return false;
 
    if (f.get_coeff(0) == Coefficient{})
      return false;
 
    DynList<Coefficient> leading_divs = positive_divisors(f.leading_coeff());
    DynList<Coefficient> constant_divs = positive_divisors(f.get_coeff(0));
    if (leading_divs.is_empty() or constant_divs.is_empty())
      return false;
 
    const double abs_lc = std::abs(static_cast<double>(f.leading_coeff()));
    double max_ratio = 0.0;
    for (auto it = f.coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        const auto &term = it.get_curr();
        if (term.first == f.degree())
          continue;
        const double ratio = std::abs(static_cast<double>(term.second)) / abs_lc;
        if (ratio > max_ratio)
          max_ratio = ratio;
      }
 
    const double root_bound = 1.0 + max_ratio;
 
    for (auto a_it = leading_divs.get_it(); a_it.has_curr(); a_it.next_ne())
      {
        const Coefficient a = a_it.get_curr();
        const int64_t b_bound
          = static_cast<int64_t>(std::ceil(2.0 * std::abs(static_cast<double>(a)) * root_bound));
 
        for (auto c_it = constant_divs.get_it(); c_it.has_curr(); c_it.next_ne())
          {
            const Coefficient c_abs = c_it.get_curr();
 
            for (Coefficient sign : {Coefficient(1), Coefficient(-1)})
              {
                const Coefficient c = sign * c_abs;
 
                for (int64_t b_i = -b_bound; b_i <= b_bound; ++b_i)
                  {
                    const Coefficient b = static_cast<Coefficient>(b_i);
                    const uint64_t gcd_ab
                      = std::gcd(polynomial_detail::abs_to_u64(a), polynomial_detail::abs_to_u64(b));
                    if (std::gcd(gcd_ab, polynomial_detail::abs_to_u64(c)) != 1)
                      continue;
 
                    Gen_Polynomial candidate;
                    candidate.set_coeff(0, c);
                    if (b != Coefficient{})
                      candidate.set_coeff(1, b);
                    candidate.set_coeff(2, a);
 
                    Gen_Polynomial q;
                    if (try_exact_candidate_division(f, candidate, q))
                      {
                        factor = candidate.primitive_part();
                        quotient = q.is_zero() ? q : q.primitive_part();
                        return true;
                      }
                  }
              }
          }
      }
 
    return false;
  }
 
  [[nodiscard]] static bool try_extract_primitive_cubic_factor(const Gen_Polynomial &f,
                                                               Gen_Polynomial &factor,
                                                               Gen_Polynomial &quotient)
    requires std::is_integral_v<Coefficient>
  {
    if (f.is_zero() or f.degree() < 3)
      return false;
 
    if (f.get_coeff(0) == Coefficient{})
      return false;
 
    DynList<Coefficient> leading_divs = positive_divisors(f.leading_coeff());
    DynList<Coefficient> constant_divs = positive_divisors(f.get_coeff(0));
    if (leading_divs.is_empty() or constant_divs.is_empty())
      return false;
 
    const double abs_lc = std::abs(static_cast<double>(f.leading_coeff()));
    double l2_norm_sq = 0.0;
    for (auto it = f.coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        const auto &term = it.get_curr();
        const double coeff_abs = std::abs(static_cast<double>(term.second));
        l2_norm_sq += coeff_abs * coeff_abs;
      }
 
    const int64_t coeff_bound = static_cast<int64_t>(std::ceil(8.0 * std::sqrt(l2_norm_sq) / abs_lc));
    if (coeff_bound <= 0)
      return false;
 
    for (auto a_it = leading_divs.get_it(); a_it.has_curr(); a_it.next_ne())
      {
        const Coefficient a = a_it.get_curr();
        if (polynomial_detail::abs_to_u64(a) > static_cast<uint64_t>(coeff_bound))
          continue;
 
        for (auto d_it = constant_divs.get_it(); d_it.has_curr(); d_it.next_ne())
          {
            const Coefficient d_abs = d_it.get_curr();
            if (polynomial_detail::abs_to_u64(d_abs) > static_cast<uint64_t>(coeff_bound))
              continue;
 
            for (Coefficient sign : {Coefficient(1), Coefficient(-1)})
              {
                const Coefficient d = sign * d_abs;
 
                for (int64_t b_i = -coeff_bound; b_i <= coeff_bound; ++b_i)
                  {
                    const Coefficient b = static_cast<Coefficient>(b_i);
                    const uint64_t gcd_ab
                      = std::gcd(polynomial_detail::abs_to_u64(a), polynomial_detail::abs_to_u64(b));
 
                    for (int64_t c_i = -coeff_bound; c_i <= coeff_bound; ++c_i)
                      {
                        const Coefficient c = static_cast<Coefficient>(c_i);
                        const uint64_t gcd_abc = std::gcd(gcd_ab, polynomial_detail::abs_to_u64(c));
                        if (std::gcd(gcd_abc, polynomial_detail::abs_to_u64(d)) != 1)
                          continue;
 
                        Gen_Polynomial candidate;
                        candidate.set_coeff(0, d);
                        if (c != Coefficient{})
                          candidate.set_coeff(1, c);
                        if (b != Coefficient{})
                          candidate.set_coeff(2, b);
                        candidate.set_coeff(3, a);
 
                        Gen_Polynomial q;
                        if (try_exact_candidate_division(f, candidate, q))
                          {
                            factor = candidate.primitive_part();
                            quotient = q.is_zero() ? q : q.primitive_part();
                            return true;
                          }
                      }
                  }
              }
          }
      }
 
    return false;
  }
 
  static DynList<Gen_Polynomial> factor_mod_p(const Gen_Polynomial &f, Coefficient p)
    requires std::is_integral_v<Coefficient>
  {
    ah_domain_error_if(p <= Coefficient(1)) << "factor_mod_p: modulus p must be > 1, got " << p;
 
    DynList<Gen_Polynomial> result;
 
    if (f.is_zero() or f.is_constant())
      return result;
 
    Gen_Polynomial remaining = reduce_coeffs_mod(f, p);
 
    for (Coefficient r = Coefficient{}; r < p; r = r + Coefficient(1))
      {
        while (
          not remaining.is_zero() and remaining.degree() > 0
          and eval_mod_u64(remaining,
                           polynomial_detail::normalize_mod_u64(r, polynomial_detail::abs_to_u64(p)),
                           polynomial_detail::abs_to_u64(p))
                == 0)
          {
            Gen_Polynomial factor;
            factor.set_coeff(0, ((-r) % p + p) % p);
            factor.set_coeff(1, Coefficient(1));
            result.append(factor);
            remaining = divide_by_monic_linear_mod(remaining, r, p);
          }
      }
 
    if (not remaining.is_zero() and not remaining.is_constant())
      result.append(remaining);
 
    return result;
  }
 
  [[nodiscard]] double mignotte_bound() const
    requires std::is_integral_v<Coefficient>
  {
    if (is_zero())
      return 0.0;
 
    const size_t d = degree();
    if (d == 0)
      return 0.0;
 
    Coefficient max_coeff = Coefficient{};
    for (auto it = coeffs.get_it(); it.has_curr(); it.next_ne())
      {
        auto &p = it.get_curr();
        Coefficient abs_c = p.second < Coefficient{} ? -p.second : p.second;
        if (abs_c > max_coeff)
          max_coeff = abs_c;
      }
 
    const auto d_real = static_cast<double>(d);
    const auto max_real = static_cast<double>(max_coeff);
    const double two_pow_d = polynomial_detail::power(2.0, d);
 
    return std::sqrt(d_real + 1.0) * two_pow_d * max_real;
  }
 
  static DynList<Gen_Polynomial> hensel_lift(const Gen_Polynomial &f,
                                             DynList<Gen_Polynomial> factors,
                                             Coefficient p,
                                             size_t levels)
    requires std::is_integral_v<Coefficient>
  {
    ah_domain_error_if(p <= Coefficient(1)) << "hensel_lift: modulus p must be > 1, got " << p;
    ah_domain_error_if(levels == 0)
      << "hensel_lift: levels must be >= 1, got " << levels;
    ah_domain_error_if(levels > 30)
      << "hensel_lift: levels must be <= 30 to avoid overflow, got " << levels;
 
    if (factors.is_empty())
      return factors;
 
    const size_t target_steps = static_cast<size_t>(1) << levels;
    const uint64_t p_u64 = polynomial_detail::abs_to_u64(p);
    const auto checked_mul_u64 = [](uint64_t lhs, uint64_t rhs, const char *name) -> uint64_t
    {
      ah_domain_error_if(rhs != 0 and lhs > std::numeric_limits<uint64_t>::max() / rhs)
        << "hensel_lift: overflow computing " << name << " (" << lhs << " * " << rhs << ")";
      return lhs * rhs;
    };
 
    uint64_t target_u64 = p_u64;
    for (size_t i = 1; i < target_steps; ++i)
      target_u64 = checked_mul_u64(target_u64, p_u64, "target_u64");
 
    ah_domain_error_if(target_u64 > static_cast<uint64_t>(std::numeric_limits<Coefficient>::max()))
      << "hensel_lift: target modulus " << target_u64
      << " does not fit in coefficient type";
    const Coefficient p_power = static_cast<Coefficient>(target_u64);
    Gen_Polynomial derivative = f.derivative();
 
    DynList<Gen_Polynomial> result;
    for (auto it = factors.get_it(); it.has_curr(); it.next_ne())
      {
        const Gen_Polynomial &factor = it.get_curr();
        bool lifted_linear_factor = false;
 
        if (factor.degree() == 1 and factor.get_coeff(1) == Coefficient(1))
          {
            uint64_t root = polynomial_detail::normalize_mod_u64(-factor.get_coeff(0), p_u64);
            uint64_t modulus = p_u64;
            bool simple = true;
 
            for (size_t step = 1; step < target_steps; ++step)
              {
                const uint64_t next_mod = checked_mul_u64(modulus, p_u64, "next_mod");
                const uint64_t f_val = eval_mod_u64(f, root, next_mod);
                const uint64_t deriv_modp = eval_mod_u64(derivative, root % p_u64, p_u64);
 
                if (f_val % modulus != 0 or std::gcd(deriv_modp, p_u64) != 1)
                  {
                    simple = false;
                    break;
                  }
 
                const uint64_t e = (f_val / modulus) % p_u64;
                const uint64_t inv = mod_inv(deriv_modp, p_u64);
                const uint64_t t = (p_u64 - mod_mul(e, inv, p_u64)) % p_u64;
 
                const uint64_t increment = checked_mul_u64(modulus, t, "root increment");
                root = static_cast<uint64_t>((static_cast<unsigned __int128>(root) + increment)
                                             % next_mod);
                modulus = next_mod;
              }
 
            if (simple)
              {
                Gen_Polynomial lifted;
                lifted.set_coeff(0,
                                 polynomial_detail::centered_from_mod_u64<Coefficient>(
                                   (target_u64 - (root % target_u64)) % target_u64, target_u64));
                lifted.set_coeff(1, Coefficient(1));
                result.append(lifted);
                lifted_linear_factor = true;
              }
          }
 
        if (lifted_linear_factor)
          continue;
 
        Gen_Polynomial normalized;
        factor.for_each_term([&normalized, p_power](size_t exp, const Coefficient &coeff)
        {
          Coefficient canonical = coeff % p_power;
          if (canonical < Coefficient{})
            canonical += p_power;
          if (canonical != Coefficient{})
            normalized.set_coeff(exp, canonical);
        });
        result.append(normalized);
      }
 
    return result;
  }
 
  [[nodiscard]] DynList<SfdTerm> factorize() const
    requires std::is_integral_v<Coefficient>
  {
    DynList<SfdTerm> result;
 
    if (is_zero() or is_constant())
      return result;
 
    const Coefficient c = content();
    if (c != Coefficient(1))
      result.append(SfdTerm{Gen_Polynomial(c), 1});
 
    // Step 1: Get primitive part
    Gen_Polynomial prim = primitive_part();
 
    // Step 2: Compute square-free decomposition
    DynList<SfdTerm> sfd_terms = prim.yun_sfd();
 
    if (sfd_terms.is_empty())
      return result;
 
    // Step 3: For each square-free group, attempt factorization
    for (auto it = sfd_terms.get_it(); it.has_curr(); it.next_ne())
      {
        auto &sfd_term = it.get_curr();
        Gen_Polynomial factor = sfd_term.factor;
        size_t mult = sfd_term.multiplicity;
        Gen_Polynomial remaining = factor.primitive_part();
 
        while (remaining.degree() > 0)
          {
            Gen_Polynomial exact_linear, quotient;
            if (not try_extract_rational_linear_factor(remaining, exact_linear, quotient))
              break;
 
            result.append(SfdTerm{exact_linear.primitive_part(), mult});
            remaining = quotient.is_zero() ? quotient : quotient.primitive_part();
          }
 
        while (remaining.degree() > 1)
          {
            Gen_Polynomial exact_quadratic, quotient;
            if (not try_extract_primitive_quadratic_factor(remaining, exact_quadratic, quotient))
              break;
 
            result.append(SfdTerm{exact_quadratic.primitive_part(), mult});
            remaining = quotient.is_zero() ? quotient : quotient.primitive_part();
          }
 
        while (remaining.degree() > 2)
          {
            Gen_Polynomial exact_cubic, quotient;
            if (not try_extract_primitive_cubic_factor(remaining, exact_cubic, quotient))
              break;
 
            result.append(SfdTerm{exact_cubic.primitive_part(), mult});
            remaining = quotient.is_zero() ? quotient : quotient.primitive_part();
          }
 
        if (remaining.is_constant() or remaining.is_zero())
          continue;
 
        // Try modular factorization with a few small primes and exact redivision.
        for (Coefficient p : {Coefficient(3), Coefficient(5), Coefficient(7), Coefficient(11)})
          {
            if (remaining.leading_coeff() % p == Coefficient{})
              continue;
 
            auto mod_factors = factor_mod_p(remaining, p);
            if (mod_factors.size() <= 1)
              continue;
 
            auto lifted = hensel_lift(remaining, mod_factors, p, 2);
 
            bool changed = false;
            for (auto jt = lifted.get_it(); jt.has_curr(); jt.next_ne())
              {
                Gen_Polynomial cand = jt.get_curr().primitive_part();
                if (cand.is_zero() or cand.is_constant() or cand.degree() != 1)
                  continue;
 
                Gen_Polynomial quotient;
                if (not try_divide_by_linear_factor(
                      remaining, cand.get_coeff(1), cand.get_coeff(0), quotient))
                  continue;
 
                result.append(SfdTerm{cand, mult});
                remaining = quotient.is_zero() ? quotient : quotient.primitive_part();
                changed = true;
 
                if (remaining.is_constant())
                  break;
              }
 
            if (changed)
              {
                if (remaining.is_constant())
                  break;
              }
          }
 
        if (not remaining.is_constant() and not remaining.is_zero())
          result.append(SfdTerm{remaining.primitive_part(), mult});
      }
 
    return result;
  }
};
 
// =================================================================
//  Free functions
// =================================================================
 
template <typename C>
[[nodiscard]] Gen_Polynomial<C> operator*(const C &s, const Gen_Polynomial<C> &p)
{
  return p * s;
}
 
template <typename C>
[[nodiscard]] Gen_Polynomial<C> operator+(const C &s, const Gen_Polynomial<C> &p)
{
  return p + s;
}
 
template <typename C>
[[nodiscard]] Gen_Polynomial<C> operator-(const C &s, const Gen_Polynomial<C> &p)
{
  return Gen_Polynomial<C>(s) - p;
}
 
template <typename C>
std::ostream &operator<<(std::ostream &out, const Gen_Polynomial<C> &p)
{
  return out << p.to_str();
}
 
// =================================================================
//  Convenience type aliases
// =================================================================
 
using Polynomial = Gen_Polynomial<double>;
 
using PolynomialF = Gen_Polynomial<float>;
 
using PolynomialLD = Gen_Polynomial<long double>;
 
using IntPolynomial = Gen_Polynomial<long long>;
 
}  // end namespace Aleph
 
#endif  // TPL_POLYNOMIAL_H