Polynomial
==========

.. cpp:class:: numpy::Polynomial

   Polynomial operations class.

Example
-------

.. code-block:: cpp

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

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

Constructors
------------

.. list-table::
   :widths: 55 25 20
   :header-rows: 1

   * - 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
---------

.. list-table::
   :widths: 40 25 15 20
   :header-rows: 1

   * - 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
--------------

.. list-table::
   :widths: 40 20 15 25
   :header-rows: 1

   * - 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
     - :ref:`View <example-polynomial-identity-0>`
   * - ``std::pair<NDArray<T>,NDArray<T>> linspace(size_tn = 100)``
     - std::pair<NDArray<T>,NDArray<T>>
     - NP_POLYNOMIAL.H:730
     - :ref:`View <example-polynomial-linspace-1>`
   * - ``void trim()``
     - void
     - NP_POLYNOMIAL.H:178
     - :ref:`View <example-polynomial-trim-2>`
   * - ``Polynomial<T> trim_public(Ttol = T(1e-10))``
     - Polynomial<T>
     - NP_POLYNOMIAL.H:824
     - :ref:`View <example-polynomial-trim_public-3>`

Indexing / Selection
--------------------

.. list-table::
   :widths: 40 20 15 25
   :header-rows: 1

   * - Signature
     - Return Type
     - Location
     - Example
   * - ``std::string get_default_printstyle()``
     - std::string
     - NP_POLYNOMIAL.H:327
     - :ref:`View <example-polynomial-get_default_printstyle-4>`
   * - ``PolynomialPrintOptions get_printoptions()``
     - PolynomialPrintOptions
     - NP_POLYNOMIAL.H:335
     - :ref:`View <example-polynomial-get_printoptions-5>`

Math Operations
---------------

.. list-table::
   :widths: 40 20 15 25
   :header-rows: 1

   * - Signature
     - Return Type
     - Location
     - Example
   * - ``Polynomial<T> power(size_tn)``
     - Polynomial<T>
     - NP_POLYNOMIAL.H:466
     - :ref:`View <example-polynomial-power-6>`
   * - ``Polynomial<T> truncate(size_tsize)``
     - Polynomial<T>
     - NP_POLYNOMIAL.H:713
     - :ref:`View <example-polynomial-truncate-7>`

Type Handling
-------------

.. list-table::
   :widths: 40 20 15 25
   :header-rows: 1

   * - 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
     - :ref:`View <example-polynomial-cast-8>`

Other Methods
-------------

.. list-table::
   :widths: 40 20 15 25
   :header-rows: 1

   * - 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
     - :ref:`View <example-polynomial-basis-9>`
   * - ``const std::vector<T>& coefficients()``
     - const std::vector<T>&
     - NP_POLYNOMIAL.H:340
     - :ref:`View <example-polynomial-coefficients-10>`
   * - ``Polynomial<T> cutdeg(size_tdeg)``
     - Polynomial<T>
     - NP_POLYNOMIAL.H:700
     - :ref:`View <example-polynomial-cutdeg-11>`
   * - ``size_t degree()``
     - size_t
     - NP_POLYNOMIAL.H:344
     - :ref:`View <example-polynomial-degree-12>`
   * - ``Polynomial<T> derivative(size_tn = 1)``
     - Polynomial<T>
     - NP_POLYNOMIAL.H:489
     - :ref:`View <example-polynomial-derivative-13>`
   * - ``std::pair<Polynomial<T>,Polynomial<T>> divmod(const Polynomial<T>&other)``
     - std::pair<Polynomial<T>,Polynomial<T>>
     - NP_POLYNOMIAL.H:625
     - :ref:`View <example-polynomial-divmod-14>`
   * - ``const std::array<T,2>& domain()``
     - const std::array<T,2>&
     - NP_POLYNOMIAL.H:341
     - :ref:`View <example-polynomial-domain-15>`
   * - ``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
     - :ref:`View <example-polynomial-fit-16>`
   * - ``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
     - :ref:`View <example-polynomial-fromroots-17>`
   * - ``bool has_samecoef(const Polynomial<T>&other, Ttol = T(1e-10))``
     - bool
     - NP_POLYNOMIAL.H:777
     - :ref:`View <example-polynomial-has_samecoef-18>`
   * - ``bool has_samedomain(const Polynomial<T>&other, Ttol = T(1e-10))``
     - bool
     - NP_POLYNOMIAL.H:803
     - :ref:`View <example-polynomial-has_samedomain-19>`
   * - ``bool has_sametype(const Polynomial<T>&other)``
     - bool
     - NP_POLYNOMIAL.H:816
     - :ref:`View <example-polynomial-has_sametype-20>`
   * - ``bool has_samewindow(const Polynomial<T>&other, Ttol = T(1e-10))``
     - bool
     - NP_POLYNOMIAL.H:809
     - :ref:`View <example-polynomial-has_samewindow-21>`
   * - ``Polynomial<T> integral(size_tn = 1, Tconstant = T(0))``
     - Polynomial<T>
     - NP_POLYNOMIAL.H:522
     - :ref:`View <example-polynomial-integral-22>`
   * - ``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
     - :ref:`View <example-polynomial-pow-23>`
   * - ``std::vector<std::complex<T>> roots()``
     - std::vector<std::complex<T>>
     - NP_POLYNOMIAL.H:555
     - :ref:`View <example-polynomial-roots-24>`
   * - ``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
     - :ref:`View <example-polynomial-set_default_printstyle-25>`
   * - ``void set_precision(intprecision)``
     - void
     - NP_POLYNOMIAL.H:331
     - 
   * - ``size_t size()``
     - size_t
     - NP_POLYNOMIAL.H:345
     - :ref:`View <example-polynomial-size-26>`
   * - ``const std::string & symbol()``
     - const std::string &
     - NP_POLYNOMIAL.H:343
     - 
   * - ``std::string toString()``
     - std::string
     - NP_POLYNOMIAL.H:948
     - :ref:`View <example-polynomial-tostring-27>`
   * - ``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
     - :ref:`View <example-polynomial-window-28>`


Code Examples
-------------

The following examples are extracted from the test suite.

.. _example-polynomial-identity-0:

.. dropdown:: identity (np_test_1_all.cpp:24002)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 23992
      :emphasize-lines: 11

                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);

.. _example-polynomial-linspace-1:

.. dropdown:: linspace (np_test_1_all.cpp:23001)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 22991
      :emphasize-lines: 11

          // 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);

.. _example-polynomial-trim-2:

.. dropdown:: trim (np_test_5_all.cpp:842)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 832
      :emphasize-lines: 11

      
            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");

.. _example-polynomial-trim_public-3:

.. dropdown:: trim_public (np_test_5_all.cpp:7618)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7608
      :emphasize-lines: 11

              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;

.. _example-polynomial-get_default_printstyle-4:

.. dropdown:: get_default_printstyle (np_test_5_all.cpp:8998)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 8988
      :emphasize-lines: 11

      
          // 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 });

.. _example-polynomial-get_printoptions-5:

.. dropdown:: get_printoptions (np_test_2_all.cpp:21259)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 21249
      :emphasize-lines: 11

        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 "

.. _example-polynomial-power-6:

.. dropdown:: power (np_test_1_all.cpp:23257)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 23247
      :emphasize-lines: 11

      
          // 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);

.. _example-polynomial-truncate-7:

.. dropdown:: truncate (np_test_5_all.cpp:7539)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7529
      :emphasize-lines: 11

              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");

.. _example-polynomial-cast-8:

.. dropdown:: cast (np_test_5_all.cpp:7638)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7628
      :emphasize-lines: 11

            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");
            }

.. _example-polynomial-basis-9:

.. dropdown:: basis (np_test_5_all.cpp:7482)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7472
      :emphasize-lines: 11

      #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");

.. _example-polynomial-coefficients-10:

.. dropdown:: coefficients (np_test_1_all.cpp:22719)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 22709
      :emphasize-lines: 11

        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;

.. _example-polynomial-cutdeg-11:

.. dropdown:: cutdeg (np_test_5_all.cpp:7519)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7509
      :emphasize-lines: 11

              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");

.. _example-polynomial-degree-12:

.. dropdown:: degree (np_test_1_all.cpp:22714)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 22704
      :emphasize-lines: 11

      
      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);
          }
      

.. _example-polynomial-derivative-13:

.. dropdown:: derivative (np_test_1_all.cpp:22884)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 22874
      :emphasize-lines: 11

          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}))";

.. _example-polynomial-divmod-14:

.. dropdown:: divmod (np_test_1_all.cpp:21140)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 21130
      :emphasize-lines: 11

          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);

.. _example-polynomial-domain-15:

.. dropdown:: domain (np_test_5_all.cpp:7640)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7630
      :emphasize-lines: 11

      
          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;

.. _example-polynomial-fit-16:

.. dropdown:: fit (np_test_1_all.cpp:23011)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 23001
      :emphasize-lines: 11

          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);

.. _example-polynomial-fromroots-17:

.. dropdown:: fromRoots (np_test_1_all.cpp:22752)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 22742
      :emphasize-lines: 11

              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);
          }

.. _example-polynomial-has_samecoef-18:

.. dropdown:: has_samecoef (np_test_5_all.cpp:7602)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7592
      :emphasize-lines: 11

            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;
          }

.. _example-polynomial-has_samedomain-19:

.. dropdown:: has_samedomain (np_test_5_all.cpp:7603)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7593
      :emphasize-lines: 11

          }
      
          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;
          }
      

.. _example-polynomial-has_sametype-20:

.. dropdown:: has_sametype (np_test_5_all.cpp:7604)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7594
      :emphasize-lines: 11

      
          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() {

.. _example-polynomial-has_samewindow-21:

.. dropdown:: has_samewindow (np_test_5_all.cpp:7603)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7593
      :emphasize-lines: 11

          }
      
          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;
          }
      

.. _example-polynomial-integral-22:

.. dropdown:: integral (np_test_1_all.cpp:10913)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 10903
      :emphasize-lines: 11

      
          // 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;

.. _example-polynomial-pow-23:

.. dropdown:: pow (np_test_1_all.cpp:15230)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 15220
      :emphasize-lines: 11

              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

.. _example-polynomial-roots-24:

.. dropdown:: roots (np_test_1_all.cpp:23043)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 23033
      :emphasize-lines: 11

      
          // 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();
            });
      

.. _example-polynomial-set_default_printstyle-25:

.. dropdown:: set_default_printstyle (np_test_5_all.cpp:8798)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 8788
      :emphasize-lines: 11

          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")) {

.. _example-polynomial-size-26:

.. dropdown:: size (np_test_1_all.cpp:47)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 37
      :emphasize-lines: 11

      
      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);

.. _example-polynomial-tostring-27:

.. dropdown:: toString (np_test_1_all.cpp:6826)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 6816
      :emphasize-lines: 11

      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");

.. _example-polynomial-window-28:

.. dropdown:: window (np_test_5_all.cpp:7642)
   :class-title: example-dropdown

   .. code-block:: cpp
      :linenos:
      :lineno-start: 7632
      :emphasize-lines: 11

            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;
          }