Polynomial

class numpy::Polynomial

Polynomial operations class.

Example

#include <numpy/np_ndarray.h>
using namespace numpy;

// Use Polynomial
Polynomial obj;
// ... operations ...

Constructors

Signature	Location	Example
Polynomial(const std::vector<T>&coefficients, const std::array<T,2>&domain = {T(-1),T(1)}, const std::array<T,2>&window = {T(-1),T(1)}, const std::string &symbol = "x")	NP_POLYNOMIAL.H:207
Polynomial(Tconstant)	NP_POLYNOMIAL.H:219
Polynomial(const Polynomial &other)	NP_POLYNOMIAL.H:222
noexcept Polynomial(Polynomial &&other)	NP_POLYNOMIAL.H:226

Operators

Signature	Return Type	Location	Example
Polynomial & operator=(const Polynomial &other)	Polynomial &	NP_POLYNOMIAL.H:223
Polynomial & noexcept operator=(Polynomial &&other)	Polynomial & noexcept	NP_POLYNOMIAL.H:227
T operator()(Tx)	T	NP_POLYNOMIAL.H:353
Polynomial<T> operator+(const Polynomial<T>&other)	Polynomial<T>	NP_POLYNOMIAL.H:378
Polynomial<T> operator-(const Polynomial<T>&other)	Polynomial<T>	NP_POLYNOMIAL.H:392
Polynomial<T> operator\*(const Polynomial<T>&other)	Polynomial<T>	NP_POLYNOMIAL.H:406
Polynomial<T> operator\*(Tscalar)	Polynomial<T>	NP_POLYNOMIAL.H:426
Polynomial<T> operator-()	Polynomial<T>	NP_POLYNOMIAL.H:436
Polynomial<T> operator+()	Polynomial<T>	NP_POLYNOMIAL.H:440
Polynomial<T>& operator+=(const Polynomial<T>&other)	Polynomial<T>&	NP_POLYNOMIAL.H:445
Polynomial<T>& operator-=(const Polynomial<T>&other)	Polynomial<T>&	NP_POLYNOMIAL.H:450
Polynomial<T>& operator\*=(const Polynomial<T>&other)	Polynomial<T>&	NP_POLYNOMIAL.H:455
Polynomial<T>& operator\*=(Tscalar)	Polynomial<T>&	NP_POLYNOMIAL.H:460
Polynomial<T> operator/(const Polynomial<T>&other)	Polynomial<T>	NP_POLYNOMIAL.H:599
Polynomial<T> operator%(const Polynomial<T>&other)	Polynomial<T>	NP_POLYNOMIAL.H:612
bool operator==(const Polynomial<T>&other)	bool	NP_POLYNOMIAL.H:661
bool operator!=(const Polynomial<T>&other)	bool	NP_POLYNOMIAL.H:673

Array Creation

Signature	Return Type	Location	Example
Polynomial<T> identity(const std::array<T,2>&domain = {T(-1),T(1)}, const std::array<T,2>&window = {T(-1),T(1)}, const std::string &symbol = "x")	Polynomial<T>	NP_POLYNOMIAL.H:693	View
std::pair<NDArray<T>,NDArray<T>> linspace(size_tn = 100)	std::pair<NDArray<T>,NDArray<T>>	NP_POLYNOMIAL.H:730	View
void trim()	void	NP_POLYNOMIAL.H:178	View
Polynomial<T> trim_public(Ttol = T(1e-10))	Polynomial<T>	NP_POLYNOMIAL.H:824	View

Indexing / Selection

Signature	Return Type	Location	Example
std::string get_default_printstyle()	std::string	NP_POLYNOMIAL.H:327	View
PolynomialPrintOptions get_printoptions()	PolynomialPrintOptions	NP_POLYNOMIAL.H:335	View

Math Operations

Signature	Return Type	Location	Example
Polynomial<T> power(size_tn)	Polynomial<T>	NP_POLYNOMIAL.H:466	View
Polynomial<T> truncate(size_tsize)	Polynomial<T>	NP_POLYNOMIAL.H:713	View

Type Handling

Signature	Return Type	Location	Example
TargetType cast(const std::array<T,2>&domain = {T(-1),T(1)}, const std::array<T,2>&window = {T(-1),T(1)})	TargetType	NP_POLYNOMIAL.H:842	View

Other Methods

Signature	Return Type	Location	Example
Polynomial<T> basis(size_tdeg, const std::array<T,2>&domain = {T(-1),T(1)}, const std::array<T,2>&window = {T(-1),T(1)}, const std::string &symbol = "x")	Polynomial<T>	NP_POLYNOMIAL.H:683	View
const std::vector<T>& coefficients()	const std::vector<T>&	NP_POLYNOMIAL.H:340	View
Polynomial<T> cutdeg(size_tdeg)	Polynomial<T>	NP_POLYNOMIAL.H:700	View
size_t degree()	size_t	NP_POLYNOMIAL.H:344	View
Polynomial<T> derivative(size_tn = 1)	Polynomial<T>	NP_POLYNOMIAL.H:489	View
std::pair<Polynomial<T>,Polynomial<T>> divmod(const Polynomial<T>&other)	std::pair<Polynomial<T>,Polynomial<T>>	NP_POLYNOMIAL.H:625	View
const std::array<T,2>& domain()	const std::array<T,2>&	NP_POLYNOMIAL.H:341	View
Polynomial<T> fit(const NDArray<T>&x, const NDArray<T>&y, size_tdegree, const std::array<T,2>&domain = {T(-1),T(1)}, const std::array<T,2>&window = {T(-1),T(1)})	Polynomial<T>	NP_POLYNOMIAL.H:249	View
Polynomial<T> fromRoots(const std::vector<T>&roots, const std::array<T,2>&domain = {T(-1),T(1)}, const std::array<T,2>&window = {T(-1),T(1)})	Polynomial<T>	NP_POLYNOMIAL.H:233	View
bool has_samecoef(const Polynomial<T>&other, Ttol = T(1e-10))	bool	NP_POLYNOMIAL.H:777	View
bool has_samedomain(const Polynomial<T>&other, Ttol = T(1e-10))	bool	NP_POLYNOMIAL.H:803	View
bool has_sametype(const Polynomial<T>&other)	bool	NP_POLYNOMIAL.H:816	View
bool has_samewindow(const Polynomial<T>&other, Ttol = T(1e-10))	bool	NP_POLYNOMIAL.H:809	View
Polynomial<T> integral(size_tn = 1, Tconstant = T(0))	Polynomial<T>	NP_POLYNOMIAL.H:522	View
T map_domain_to_window(Tx)	T	NP_POLYNOMIAL.H:188
T map_window_to_domain(Tx)	T	NP_POLYNOMIAL.H:196
Polynomial<T> pow(intexponent)	Polynomial<T>	NP_POLYNOMIAL.H:642	View
std::vector<std::complex<T>> roots()	std::vector<std::complex<T>>	NP_POLYNOMIAL.H:555	View
void setDomain(const std::array<T,2>&domain)	void	NP_POLYNOMIAL.H:348
void setSymbol(const std::string &symbol)	void	NP_POLYNOMIAL.H:350
void setWindow(const std::array<T,2>&window)	void	NP_POLYNOMIAL.H:349
void set_default_printstyle(const std::string &style)	void	NP_POLYNOMIAL.H:323	View
void set_precision(intprecision)	void	NP_POLYNOMIAL.H:331
size_t size()	size_t	NP_POLYNOMIAL.H:345	View
const std::string & symbol()	const std::string &	NP_POLYNOMIAL.H:343
std::string toString()	std::string	NP_POLYNOMIAL.H:948	View
std::string toString_ascii()	std::string	NP_POLYNOMIAL.H:854
std::string toString_unicode()	std::string	NP_POLYNOMIAL.H:901
const std::array<T,2>& window()	const std::array<T,2>&	NP_POLYNOMIAL.H:342	View

Code Examples

The following examples are extracted from the test suite.

identity (np_test_1_all.cpp:24002)

          throw std::runtime_error(description);
      }
      if (!((eye_offset.getShape() == std::vector<size_t>{3, 4}))) {
          std::string description = std::string("testArrayCreationSignatures():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((eye_offset.getShape() == std::vector<size_t>{3, 4}))";
          std::cout << std::string("[FAIL] ") + description << std::endl;
          throw std::runtime_error(description);
      }
      // std::cout << "[OK] eye function with optional parameters works correctly\n";

      // Test identity
      auto identity_result = identity<int>(4);
      if (!((identity_result.getShape() == std::vector<size_t>{4, 4}))) {
          std::string description = std::string("testArrayCreationSignatures():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((identity_result.getShape() == std::vector<size_t>{4, 4}))";
          std::cout << std::string("[FAIL] ") + description << std::endl;
          throw std::runtime_error(description);
      }
      // std::cout << "[OK] identity function has correct signature\n";

      // Test diag with different use cases
      auto vector1d = createInt32Array({ 3 }, 5);
      auto matrix2d = createInt32Array({ 3, 3 }, 0);

linspace (np_test_1_all.cpp:23001)

    // assert(std::abs(p_roots[2] - 1.0) < 1e-10);

    // std::cout << "[OK] Polynomial utilities work correctly" << std::endl;
    std::cout << " -> tests passed" << std::endl;
  }

  void test_polynomial_fitting() {
    std::cout << "========= test_polynomial_fitting =======================";

    // Create test data for y = 1 + 2x + 3x^2
    auto x = numpy::linspace<double>(-1.0, 1.0, 10);
    auto y_true = numpy::NDArray<double>::createZeros({ 10 });

    for (size_t i = 0; i < 10; ++i) {
      double x_val = x.getElementAt({ i });
      double y_val = 1.0 + 2.0 * x_val + 3.0 * x_val * x_val;
      y_true.setElementAt({ i }, y_val);
    }

    // Fit polynomial
    auto fitted_poly = numpy::Polynomial<double>::fit(x, y_true, 2);

trim (np_test_5_all.cpp:842)

      std::cout << " -> tests passed" << std::endl;
    }
#endif

    // ============================================================================
    // I/O UTILITIES TESTS
    // ============================================================================

    void np_test_io_trim() {
      std::cout << "========= trim() utility =======================================";

      // trim() might be in namespace or as a helper function
      // Test basic string trimming
      std::string test1 = "  hello  ";
      std::string test2 = "\t\nworld\t\n";
      std::string test3 = "no_trim";

      // Since trim may not be directly accessible, we'll test the concept
      auto trim = [](const std::string& str) {
        size_t start = str.find_first_not_of(" \t\n\r");

trim_public (np_test_5_all.cpp:7618)

        throw std::runtime_error("np_test_has_same_methods failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }

    void np_test_trim_public() {
      std::cout << "========= trim_public: remove trailing zeros =======================";

      numpy::Polynomial<double> p({ 1.0, 2.0, 0.0, 0.0 });
      auto trimmed = p.trim_public();

      auto coefs = trimmed.coefficients();
      bool passed = (coefs.size() == 2);

      if (!passed) {
        std::cout << "  [FAIL] : trim_public failed";
        throw std::runtime_error("np_test_trim_public failed");
      }

      std::cout << " -> tests passed" << std::endl;

get_default_printstyle (np_test_5_all.cpp:8998)

    // Test 2: Switch to Unicode via static method
    numpy::Polynomial<double>::set_default_printstyle("unicode");
    std::string result2 = p1.toString();
    // std::cout << "Static method Unicode: " << result2 << std::endl;

    if (!verify_contains(result2, "\u00B7", "Static method Unicode")) {
      throw std::runtime_error("Polynomial<T>::set_default_printstyle('unicode') failed");
    }

    // Test 3: Verify Polynomial<T>::get_default_printstyle() static method
    std::string current_style = numpy::Polynomial<double>::get_default_printstyle();
    // std::cout << "Current style via static method: " << current_style << std::endl;

    if (current_style != "unicode") {
      throw std::runtime_error("Polynomial<T>::get_default_printstyle() returned wrong value");
    }

    // Test 4: Verify all polynomial types share the same global print style
    numpy::Chebyshev<double>::set_default_printstyle("ascii");
    numpy::Chebyshev<double> cheb({ 1.0, 2.0 });

get_printoptions (np_test_2_all.cpp:21259)

  namespace numpy_tests_printoptions {

    // ============================================================================
    // BASIC FUNCTIONALITY TESTS
    // ============================================================================

    void np_test_printoptions_basic() {
      std::cout << "========= printoptions: basic get/set ====";

      // Get default options
      auto default_opts = numpy::get_printoptions();

      // Verify defaults
      if (default_opts.precision != 8) {
        std::cout << "  [FAIL] : in np_test_printoptions_basic() : default precision should be 8, got "
          << default_opts.precision << std::endl;
        throw std::runtime_error("np_test_printoptions_basic failed: default precision incorrect");
      }

      if (default_opts.threshold != 1000) {
        std::cout << "  [FAIL] : in np_test_printoptions_basic() : default threshold should be 1000, got "

power (np_test_1_all.cpp:23257)

    // std::cout << "[OK] Laguerre polynomials work correctly" << std::endl;
    std::cout << " -> tests passed" << std::endl;
  }

  void test_polynomial_advanced() {
    std::cout << "========= test_polynomial_advanced =======================";

    // Test polynomial power
    numpy::Polynomial<double> p({ 1, 1 });  // 1 + x
    auto p_squared = p.power(2);  // (1 + x)^2 = 1 + 2x + x^2
    if (!((p_squared.coefficients() == std::vector<double>{1, 2, 1}))) {
        std::string description = std::string("test_polynomial_advanced():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((p_squared.coefficients() == std::vector<double>{1, 2, 1}))";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }

    // Test polynomial division
    std::vector<double> dividend = { -1, 0, 1 };  // x^2 - 1
    std::vector<double> divisor = { -1, 1 };      // x - 1
    auto [quotient, remainder] = numpy::polydiv(dividend, divisor);

truncate (np_test_5_all.cpp:7539)

        throw std::runtime_error("np_test_cutdeg failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }

    void np_test_truncate() {
      std::cout << "========= truncate: truncate to size =======================";

      numpy::Polynomial<double> p({ 1.0, 2.0, 3.0, 4.0, 5.0 });
      auto truncated = p.truncate(3);  // Keep only first 3 coefficients

      auto coefs = truncated.coefficients();
      bool passed = (coefs.size() == 3 &&
        std::abs(coefs[0] - 1.0) < 1e-10 &&
        std::abs(coefs[1] - 2.0) < 1e-10 &&
        std::abs(coefs[2] - 3.0) < 1e-10);

      if (!passed) {
        std::cout << "  [FAIL] : truncate failed";
        throw std::runtime_error("np_test_truncate failed");

cast (np_test_5_all.cpp:7638)

      std::cout << " -> tests passed" << std::endl;
    }

    void np_test_cast() {
      std::cout << "========= cast: convert polynomial type =======================";

      numpy::Polynomial<double> p({ 1.0, 2.0, 3.0 });
      std::array<double, 2> new_domain = { 0.0, 5.0 };
      std::array<double, 2> new_window = { -2.0, 2.0 };

      auto casted = p.cast(new_domain, new_window);

      bool passed = (std::abs(casted.domain()[0] - new_domain[0]) < 1e-10 &&
        std::abs(casted.domain()[1] - new_domain[1]) < 1e-10 &&
        std::abs(casted.window()[0] - new_window[0]) < 1e-10 &&
        std::abs(casted.window()[1] - new_window[1]) < 1e-10);

      if (!passed) {
        std::cout << "  [FAIL] : cast failed";
        throw std::runtime_error("np_test_cast failed");
      }

basis (np_test_5_all.cpp:7482)

#include <cmath>
#include "../numpy/np_polynomial.h"

namespace numpy_tests {

  namespace numpy_tests_polynomial_methods {

    void np_test_basis() {
      std::cout << "========= basis: generate basis polynomial =======================";

      auto b0 = numpy::Polynomial<double>::basis(0);  // Should be [1]
      auto b1 = numpy::Polynomial<double>::basis(1);  // Should be [0, 1]
      auto b2 = numpy::Polynomial<double>::basis(2);  // Should be [0, 0, 1]

      bool passed = (b0.coefficients().size() == 1 && b0.coefficients()[0] == 1.0 &&
        b1.coefficients().size() == 2 && b1.coefficients()[1] == 1.0 &&
        b2.coefficients().size() == 3 && b2.coefficients()[2] == 1.0);

      if (!passed) {
        std::cout << "  [FAIL] : basis polynomials incorrect";
        throw std::runtime_error("np_test_basis failed");

coefficients (np_test_1_all.cpp:22719)

  void test_polynomial_creation() {
    std::cout << "========= test_polynomial_creation =======================";

    // Test default constructor
    numpy::Polynomial<double> p1;
    if (!(p1.degree() == 0)) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(p1.degree() == 0)";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }
    if (!((p1.coefficients() == std::vector<double>{0.0}))) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((p1.coefficients() == std::vector<double>{0.0}))";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }

    // Test coefficient constructor
    numpy::Polynomial<double> p2({ 1, 2, 3 });  // 1 + 2x + 3x^2
    if (!(p2.degree() == 2)) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(p2.degree() == 2)";
        std::cout << std::string("[FAIL] ") + description << std::endl;

cutdeg (np_test_5_all.cpp:7519)

        throw std::runtime_error("np_test_identity failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }

    void np_test_cutdeg() {
      std::cout << "========= cutdeg: truncate to degree =======================";

      numpy::Polynomial<double> p({ 1.0, 2.0, 3.0, 4.0 });  // 1 + 2x + 3x^2 + 4x^3
      auto truncated = p.cutdeg(2);  // Should be 1 + 2x + 3x^2

      auto coefs = truncated.coefficients();
      bool passed = (coefs.size() == 3 &&
        std::abs(coefs[0] - 1.0) < 1e-10 &&
        std::abs(coefs[1] - 2.0) < 1e-10 &&
        std::abs(coefs[2] - 3.0) < 1e-10);

      if (!passed) {
        std::cout << "  [FAIL] : cutdeg failed";
        throw std::runtime_error("np_test_cutdeg failed");

degree (np_test_1_all.cpp:22714)

namespace numpy_tests {

  using namespace numpy;

  void test_polynomial_creation() {
    std::cout << "========= test_polynomial_creation =======================";

    // Test default constructor
    numpy::Polynomial<double> p1;
    if (!(p1.degree() == 0)) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(p1.degree() == 0)";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }
    if (!((p1.coefficients() == std::vector<double>{0.0}))) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((p1.coefficients() == std::vector<double>{0.0}))";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }

derivative (np_test_1_all.cpp:22884)

    std::cout << " -> tests passed" << std::endl;
  }

  void test_polynomial_calculus() {
    std::cout << "========= test_polynomial_calculus =======================";

    // Test derivative: p(x) = 1 + 2x + 3x^2 + 4x^3
    numpy::Polynomial<double> p({ 1, 2, 3, 4 });

    // First derivative: p'(x) = 2 + 6x + 12x^2
    auto p_deriv = p.derivative();
    if (!((p_deriv.coefficients() == std::vector<double>{2, 6, 12}))) {
        std::string description = std::string("test_polynomial_calculus():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((p_deriv.coefficients() == std::vector<double>{2, 6, 12}))";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }

    // Second derivative: p''(x) = 6 + 24x
    auto p_deriv2 = p.derivative(2);
    if (!((p_deriv2.coefficients() == std::vector<double>{6, 24}))) {
        std::string description = std::string("test_polynomial_calculus():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((p_deriv2.coefficients() == std::vector<double>{6, 24}))";

divmod (np_test_1_all.cpp:21140)

    auto a = NDArray<float64>({ 3 });
    a.setElementAt({ 0 }, 7.0);
    a.setElementAt({ 1 }, -7.0);
    a.setElementAt({ 2 }, 7.5);

    auto b = NDArray<float64>({ 3 });
    b.setElementAt({ 0 }, 3.0);
    b.setElementAt({ 1 }, 3.0);
    b.setElementAt({ 2 }, 2.0);

    auto [quotient, remainder] = divmod(a, b);

    if (!(quotient.getElementAt({ 0 }) == 2.0)) {
        std::string description = std::string("testDivmodFunction():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(quotient.getElementAt({ 0 }) == 2.0)";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }
    if (!(std::abs(remainder.getElementAt({ 0 }) - 1.0) < 1e-10)) {
        std::string description = std::string("testDivmodFunction():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(std::abs(remainder.getElementAt({ 0 }) - 1.0) < 1e-10)";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);

domain (np_test_5_all.cpp:7640)

    void np_test_cast() {
      std::cout << "========= cast: convert polynomial type =======================";

      numpy::Polynomial<double> p({ 1.0, 2.0, 3.0 });
      std::array<double, 2> new_domain = { 0.0, 5.0 };
      std::array<double, 2> new_window = { -2.0, 2.0 };

      auto casted = p.cast(new_domain, new_window);

      bool passed = (std::abs(casted.domain()[0] - new_domain[0]) < 1e-10 &&
        std::abs(casted.domain()[1] - new_domain[1]) < 1e-10 &&
        std::abs(casted.window()[0] - new_window[0]) < 1e-10 &&
        std::abs(casted.window()[1] - new_window[1]) < 1e-10);

      if (!passed) {
        std::cout << "  [FAIL] : cast failed";
        throw std::runtime_error("np_test_cast failed");
      }

      std::cout << " -> tests passed" << std::endl;

fit (np_test_1_all.cpp:23011)

    auto x = numpy::linspace<double>(-1.0, 1.0, 10);
    auto y_true = numpy::NDArray<double>::createZeros({ 10 });

    for (size_t i = 0; i < 10; ++i) {
      double x_val = x.getElementAt({ i });
      double y_val = 1.0 + 2.0 * x_val + 3.0 * x_val * x_val;
      y_true.setElementAt({ i }, y_val);
    }

    // Fit polynomial
    auto fitted_poly = numpy::Polynomial<double>::fit(x, y_true, 2);

    // Debug output
    // std::cout << "Fitted coefficients: ["
      // << fitted_poly.coefficients()[0] << ", "
      // << fitted_poly.coefficients()[1] << ", "
      // << fitted_poly.coefficients()[2] << "]" << std::endl;

    // Check that fitted coefficients are close to original
    // TODO: Fix fitting algorithm - temporarily disabled
    // assert(std::abs(fitted_poly.coefficients()[0] - 1.0) < 1e-10);

fromRoots (np_test_1_all.cpp:22752)

        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }
    if (!((p3.coefficients() == std::vector<double>{5.0}))) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !((p3.coefficients() == std::vector<double>{5.0}))";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }

    // Test polynomial from roots
    auto p4 = numpy::Polynomial<double>::fromRoots({ 1.0, -1.0 });  // (x-1)(x+1) = x^2 - 1
    if (!(p4.degree() == 2)) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(p4.degree() == 2)";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }
    if (!(std::abs(p4.coefficients()[0] + 1.0) < 1e-10)) {
        std::string description = std::string("test_polynomial_creation():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(std::abs(p4.coefficients()[0] + 1.0) < 1e-10)";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }

has_samecoef (np_test_5_all.cpp:7602)

      std::cout << " -> tests passed" << std::endl;
    }

    void np_test_has_same_methods() {
      std::cout << "========= has_same* methods =======================";

      numpy::Polynomial<double> p1({ 1.0, 2.0 });
      numpy::Polynomial<double> p2({ 1.0, 2.0 });
      numpy::Polynomial<double> p3({ 1.0, 3.0 });

      bool passed = (p1.has_samecoef(p2) && !p1.has_samecoef(p3) &&
        p1.has_samedomain(p2) && p1.has_samewindow(p2) &&
        p1.has_sametype(p2) && p1.has_sametype(p3));

      if (!passed) {
        std::cout << "  [FAIL] : has_same methods failed";
        throw std::runtime_error("np_test_has_same_methods failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }

has_samedomain (np_test_5_all.cpp:7603)

    }

    void np_test_has_same_methods() {
      std::cout << "========= has_same* methods =======================";

      numpy::Polynomial<double> p1({ 1.0, 2.0 });
      numpy::Polynomial<double> p2({ 1.0, 2.0 });
      numpy::Polynomial<double> p3({ 1.0, 3.0 });

      bool passed = (p1.has_samecoef(p2) && !p1.has_samecoef(p3) &&
        p1.has_samedomain(p2) && p1.has_samewindow(p2) &&
        p1.has_sametype(p2) && p1.has_sametype(p3));

      if (!passed) {
        std::cout << "  [FAIL] : has_same methods failed";
        throw std::runtime_error("np_test_has_same_methods failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }

has_sametype (np_test_5_all.cpp:7604)

    void np_test_has_same_methods() {
      std::cout << "========= has_same* methods =======================";

      numpy::Polynomial<double> p1({ 1.0, 2.0 });
      numpy::Polynomial<double> p2({ 1.0, 2.0 });
      numpy::Polynomial<double> p3({ 1.0, 3.0 });

      bool passed = (p1.has_samecoef(p2) && !p1.has_samecoef(p3) &&
        p1.has_samedomain(p2) && p1.has_samewindow(p2) &&
        p1.has_sametype(p2) && p1.has_sametype(p3));

      if (!passed) {
        std::cout << "  [FAIL] : has_same methods failed";
        throw std::runtime_error("np_test_has_same_methods failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }

    void np_test_trim_public() {

has_samewindow (np_test_5_all.cpp:7603)

    }

    void np_test_has_same_methods() {
      std::cout << "========= has_same* methods =======================";

      numpy::Polynomial<double> p1({ 1.0, 2.0 });
      numpy::Polynomial<double> p2({ 1.0, 2.0 });
      numpy::Polynomial<double> p3({ 1.0, 3.0 });

      bool passed = (p1.has_samecoef(p2) && !p1.has_samecoef(p3) &&
        p1.has_samedomain(p2) && p1.has_samewindow(p2) &&
        p1.has_sametype(p2) && p1.has_sametype(p3));

      if (!passed) {
        std::cout << "  [FAIL] : has_same methods failed";
        throw std::runtime_error("np_test_has_same_methods failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }

integral (np_test_1_all.cpp:10913)

    // Si(π/2) ≈ 1.851937 (approximately)
    double si_pi2 = si(M_PI / 2.0);
    if (!(!std::isnan(si_pi2) && std::isfinite(si_pi2))) {
        std::string description = std::string("testTrigonometricIntegrals():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(!std::isnan(si_pi2) && std::isfinite(si_pi2))";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }
    // std::cout << "[OK] Sine integral basic functionality\n";

    // Test cosine integral (undefined at x = 0)
    if (!(std::isnan(ci(0.0)))) {
        std::string description = std::string("testTrigonometricIntegrals():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(std::isnan(ci(0.0)))";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }

    double ci_1 = ci(1.0);
    if (!(!std::isnan(ci_1) && std::isfinite(ci_1))) {
        std::string description = std::string("testTrigonometricIntegrals():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(!std::isnan(ci_1) && std::isfinite(ci_1))";
        std::cout << std::string("[FAIL] ") + description << std::endl;

pow (np_test_1_all.cpp:15230)

        std::cout << "[FAIL] Broadcasting validation failed" << std::endl;
        errors++;
      }

      // Test 11: Performance with caching
      std::cout << "========= test_caching_performance =======================" ;

      int call_count = 0;
      auto expensive_func = [&call_count](double x) {
        call_count++;
        return std::pow(x, 3) + std::sin(x);  // Simulate expensive computation
        };

      auto cached_vectorized = numpy::vectorize(expensive_func, true, 100);

      numpy::NDArray<double> test_array({ 3 });
      test_array.setElementAt({ 0 }, 1.0);
      test_array.setElementAt({ 1 }, 2.0);
      test_array.setElementAt({ 2 }, 3.0);

      // First call

roots (np_test_1_all.cpp:23043)

    // std::cout << "[OK] Polynomial fitting works correctly" << std::endl;
    std::cout << " -> tests passed" << std::endl;
  }

  void test_polynomial_roots() {
    std::cout << "========= test_polynomial_roots =======================";

    // Test quadratic: x^2 - 1 = (x-1)(x+1), roots should be ±1
    numpy::Polynomial<double> quad({ -1, 0, 1 });
    auto roots = quad.roots();

    // TODO: Fix root finding algorithm - temporarily disabled
    // assert(roots.size() == 2);

    // Sort roots by real part
    std::sort(roots.begin(), roots.end(),
      [](const std::complex<double>& a, const std::complex<double>& b) {
        return a.real() < b.real();
      });

set_default_printstyle (np_test_5_all.cpp:8798)

    std::cout << " -> tests passed" << std::endl;
  }

  // ===========================
  // Polynomial Tests
  // ===========================

  void np_test_toString_polynomial_basic() {
    std::cout << "========= Polynomial toString() basic ============================";

    numpy::set_default_printstyle("ascii");
    numpy::set_polynomial_precision(2);

    numpy::Polynomial<double> p({ 1.0, 2.0, 3.0 });
    std::string result = p.toString();

    // std::cout << "Polynomial: " << result << std::endl;

    verify_not_empty(result, "Polynomial toString()");

    if (!verify_contains(result, "x", "Polynomial")) {

size (np_test_1_all.cpp:47)

using namespace numpy;
using namespace numpy::benchmark;
using namespace numpy::char_;

// Helper functions for array comparison tests
namespace {
    const double TOLERANCE = 1e-10;

    bool allTrue(const NDArray<bool>& arr) {
        std::vector<size_t> indices(arr.getShape().size(), 0);
        do {
            if (!arr.getElementAt(indices)) {
                return false;
            }
        } while (incrementIndices(indices, arr.getShape()));
        return true;
    }

    bool allFalse(const NDArray<bool>& arr) {
        std::vector<size_t> indices(arr.getShape().size(), 0);

toString (np_test_1_all.cpp:6826)

void testDatetime64CreationDateTime() {
    std::cout << "========= testDatetime64CreationDateTime =======================";

    // Test default constructor
    datetime64 default_dt;
    if (!(default_dt.isNaT())) {
        std::string description = std::string("testDatetime64CreationDateTime():") + __FILE__ + ":" + std::to_string(__LINE__) + ": !(default_dt.isNaT())";
        std::cout << std::string("[FAIL] ") + description << std::endl;
        throw std::runtime_error(description);
    }
    // std::cout << "Default datetime64 is NaT: " << default_dt.toString() << std::endl;

    // Test value constructor
    datetime64 epoch_day(0, DateTimeUnit::Day);
    // std::cout << "Epoch day: " << epoch_day.toString() << std::endl;

    // Test string constructor
    datetime64 date_from_string("2024-01-15");
    // std::cout << "Date from string '2024-01-15': " << date_from_string.toString() << std::endl;

    datetime64 datetime_from_string("2024-01-15T10:30:45");

window (np_test_5_all.cpp:7642)

      std::cout << "========= cast: convert polynomial type =======================";

      numpy::Polynomial<double> p({ 1.0, 2.0, 3.0 });
      std::array<double, 2> new_domain = { 0.0, 5.0 };
      std::array<double, 2> new_window = { -2.0, 2.0 };

      auto casted = p.cast(new_domain, new_window);

      bool passed = (std::abs(casted.domain()[0] - new_domain[0]) < 1e-10 &&
        std::abs(casted.domain()[1] - new_domain[1]) < 1e-10 &&
        std::abs(casted.window()[0] - new_window[0]) < 1e-10 &&
        std::abs(casted.window()[1] - new_window[1]) < 1e-10);

      if (!passed) {
        std::cout << "  [FAIL] : cast failed";
        throw std::runtime_error("np_test_cast failed");
      }

      std::cout << " -> tests passed" << std::endl;
    }