chore: update project configuration in pyproject.toml

- Drop support for Python 3.8.
- Bump minimum dependencies to `numpy>=2.0` and `numba>=0.60`.
- Replace legacy NumPy polynomial functions (polyder, polyadd, etc.) with
  the modern `numpy.polynomial.Polynomial` class for stable calculus and
  arithmetic operations.
- Remove `_strip_leading_zeros` as the modern API handles degree alignment
  and zero-trimming natively.
- Retain highest-to-lowest `self.coefficients` array to maintain compatibility
  with Numba and CUDA GPU kernels evaluating via Horner's method.
- Fix `TypeError` in `quadratic_solve` by evaluating `b.real >= 0` to safely
  support complex coefficients.
- Update `quadratic_solve` return type hint to `np.ndarray`.

.gitea/workflows/publish-python.yml

@@ -16,7 +16,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v5
       with:
-        python-version: '3.12'
+        python-version: '3.14'
 
     - name: Install dependencies
       run: |

.gitea/workflows/test.yml

@@ -6,7 +6,7 @@ jobs:
     runs-on: nvidia-cuda
     strategy:
       matrix:
-        python-version: ["3.8", "3.10", "3.12"]
+        python-version: ["3.10", "3.12", "3.14"]
 
     steps:
     - name: Checkout repository

CONTRIBUTING.md

@@ -81,7 +81,7 @@ You can replicate our CI environment exactly to minimize "it works on my machine
 * **CUDA Toolkit:** 12.5.1
 * **Base Docker Image:** `nvidia/cuda:12.5.1-devel-ubuntu24.04`
 * **Node.js Version:** 20.x
-* **Python Versions Tested:** 3.8, 3.10, 3.12
+* **Python Versions Tested:** 3.8, 3.10, 3.12, 3.14
 * **CI Docker Image:** You can pull the exact image used by our runners from Docker Hub:
 
 ```bash

README.md

@@ -51,6 +51,8 @@ from polysolve import Function, GA_Options
 #    Coefficients can be integers or floats.
 f1 = Function(largest_exponent=2)
 f1.set_coeffs([2, -3, -5])
+f2 = Function(largest_exponent=2)
+f2.set_coeffs([3, -4, -5])
 
 print(f"Function f1: {f1}")
 # > Function f1: 2x^2 - 3x - 5
@@ -70,20 +72,30 @@ ddf1 = f1.nth_derivative(2)
 print(f"2nd Derivative of f1: {ddf1}")
 # > Derivative of f1: 4
 
-# 5. Find roots analytically using the quadratic formula
+# 5. Find the integral: x^3 - 2x^2 - 5x + 0
+igl = f2.integral()
+print(f"Integral of f2: {igl}")
+# > integral of f2: x^3 - 2x^2 - 5x + 0
+
+# 6. Find the integral with custom constant of integration: x^3 - 2x^2 - 5x + 2
+igl = f2.integral(constant_of_integration=2)
+print(f"Integral of f2: {igl}")
+# > integral of f2: x^3 - 2x^2 - 5x + 2
+
+# 7. Find roots analytically using the quadratic formula
 #    This is exact and fast for degree-2 polynomials.
 roots_analytic = f1.quadratic_solve()
 print(f"Analytic roots: {sorted(roots_analytic)}")
 # > Analytic roots: [-1.0, 2.5]
 
-# 6. Find REAL roots with the genetic algorithm (Numba CPU)
+# 8. Find REAL roots with the genetic algorithm (Numba CPU)
 #    This is the default, JIT-compiled CPU solver.
 ga_opts = GA_Options(num_of_generations=20)
 roots_ga = f1.get_real_roots(ga_opts, use_cuda=False)
 print(f"Approximate real roots: {roots_ga[:2]}")
 # > Approximate real roots: [-1.000..., 2.500...]
 
-# 7. Find ALL roots (Real + Complex)
+# 9. Find ALL roots (Real + Complex)
 #    Use get_roots() to search the complex plane.
 f_complex = Function(2, [1, 0, 1]) # x^2 + 1
 roots_all = f_complex.get_roots(ga_opts)

pyproject.toml

@@ -5,14 +5,14 @@ build-backend = "setuptools.build_meta"
 [project]
 # --- Core Metadata ---
 name = "polysolve"
-version = "0.7.1"
+version = "0.8.1"
 authors = [
   { name="Jonathan Rampersad", email="[email protected]" },
 ]
 description = "A Python library for representing, manipulating, and solving exponential functions using analytical methods and genetic algorithms, with optional CUDA acceleration."
 readme = "README.md"
-requires-python = ">=3.8"
-license = { file="LICENSE" }
+requires-python = ">=3.10"
+license = { file = "LICENSE" }
 keywords = ["math", "polynomial", "genetic algorithm", "cuda", "equation solver"]
 
 # --- Classifiers ---
@@ -23,23 +23,23 @@ classifiers = [
     "License :: OSI Approved :: MIT License",
     "Operating System :: OS Independent",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
     "Topic :: Scientific/Engineering :: Mathematics",
 ]
 
 # --- Dependencies ---
 dependencies = [
-    "numpy>=1.21",
-    "numba"
+    "numpy>=2.0",
+    "numba>=0.60"
 ]
 
 # --- Optional Dependencies (Extras) ---
 [project.optional-dependencies]
-cuda12 = ["cupy-cuda12x"]
+cuda12 = ["cupy-cuda12x[ctk]"]
 dev = ["pytest"]
 
 [project.urls]

src/polysolve/__init__.py

@@ -1,6 +1,7 @@
 import math
 import cmath
 import numpy as np
+from numpy.polynomial import Polynomial
 import numba
 from dataclasses import dataclass
 from typing import List, Optional, Union
@@ -383,6 +384,9 @@ class Function:
             largest_exponent (int): The largest exponent (n) in the function.
         """
         self._largest_exponent = largest_exponent
+        self.poly_obj = None
+        self._initialized = False
+
         if coefficients is not None:
             self.set_coeffs(coefficients)
             # Verify user provided exponent matches if they provided both
@@ -390,7 +394,6 @@ class Function:
                 raise ValueError("Provided largest_exponent does not match coefficient list length.")
         elif largest_exponent is not None:
             self.coefficients = None
-            self._initialized = False
         else:
             raise ValueError("Must provide either coefficients or largest_exponent.")
 
@@ -425,6 +428,9 @@ class Function:
             target_dtype = np.float64
 
         self.coefficients = np.array(coefficients, dtype=target_dtype)
+        reversed_coeffs = coefficients[::-1]
+        self.poly_obj = Polynomial(reversed_coeffs)
+
         self._largest_exponent = len(coefficients) - 1
         self._initialized = True
 
@@ -454,7 +460,7 @@ class Function:
             float: The resulting y-value.
         """
         self._check_initialized()
-        return np.polyval(self.coefficients, x_val)
+        return self.poly_obj(x_val)
 
     def differential(self) -> 'Function':
         """
@@ -489,10 +495,11 @@ class Function:
             diff_func.set_coeffs([0])
             return diff_func
         
-        derivative_coefficients = np.polyder(self.coefficients)
+        deriv_poly = self.poly_obj.deriv(1)
+        new_coeffs = deriv_poly.coef[::-1].tolist()
         
         diff_func = Function(self._largest_exponent - 1)
-        diff_func.set_coeffs(derivative_coefficients.tolist())
+        diff_func.set_coeffs(new_coeffs)
         return diff_func
     
 
@@ -512,18 +519,30 @@ class Function:
             raise ValueError("Derivative order 'n' must be a positive integer.")
 
         if n > self.largest_exponent:
-            function = Function(0)
-            function.set_coeffs([0])
-            return function
+            diff_func = Function(0)
+            diff_func.set_coeffs([0])
+            return diff_func
 
-        if n == 1:
-            return self.derivative()
+        deriv_poly = self.poly_obj.deriv(n)
+        new_coeffs = deriv_poly.coef[::-1].tolist()
         
-        function = self
-        for _ in range(n):
-            function = function.derivative()
+        new_degree = self._largest_exponent - n
+        diff_func = Function(new_degree)
+        diff_func.set_coeffs(new_coeffs)
+        return diff_func
 
-        return function
+
+    def integral(self, constant_of_integration: float = 0.0) -> 'Function':
+        """Calculates the indefinite integral (antiderivative)."""
+        self._check_initialized()
+    
+        # Integration naturally increases the degree by 1
+        integ_poly = self.poly_obj.integ(m=1, k=[constant_of_integration])
+        new_coeffs = integ_poly.coef[::-1].tolist()
+    
+        integral_func = Function(self._largest_exponent + 1)
+        integral_func.set_coeffs(new_coeffs)
+        return integral_func
 
 
     def get_real_roots(self, options: Optional[GA_Options] = None, use_cuda: bool = False) -> np.ndarray:
@@ -1139,49 +1158,10 @@ class Function:
 
     def __str__(self) -> str:
         """Returns a human-readable string representation of the function."""
-        self._check_initialized()
-        parts = []
-        for i, c in enumerate(self.coefficients):
-            if c == 0:
-                continue
+        if not self._initialized:
+            return "Uninitialized Function"
 
-            power = self._largest_exponent - i
-            
-            # Coefficient part
-            coeff_val = c
-            if c == int(c):
-                coeff_val = int(c)
-
-            if coeff_val == 1 and power != 0:
-                coeff = ""
-            elif coeff_val == -1 and power != 0:
-                coeff = "-"
-            else:
-                coeff = str(coeff_val)
-
-            # Variable part
-            if power == 0:
-                var = ""
-            elif power == 1:
-                var = "x"
-            else:
-                var = f"x^{power}"
-
-            # Add sign for non-leading terms
-            sign = ""
-            if i > 0:
-                sign = " + " if c > 0 else " - "
-                coeff = str(abs(coeff_val))
-                if abs(c) == 1 and power != 0:
-                    coeff = "" # Don't show 1 for non-constant terms
-
-            parts.append(f"{sign}{coeff}{var}")
-        
-        # Join parts and clean up
-        result = "".join(parts)
-        if result.startswith(" + "):
-            result = result[3:]
-        return result if result else "0"
+        return str(self.poly_obj)
 
     def __repr__(self) -> str:
         return f"Function(str='{self}')"
@@ -1191,11 +1171,12 @@ class Function:
         self._check_initialized()
         other._check_initialized()
 
-        new_coefficients = np.polyadd(self.coefficients, other.coefficients)
-        new_coefficients = self._strip_leading_zeros(new_coefficients)
+        # The Polynomial objects handle degree alignment and zero-trimming natively
+        new_poly = self.poly_obj + other.poly_obj
+        new_coeffs = new_poly.coef[::-1].tolist()
         
-        result_func = Function(len(new_coefficients) - 1)
-        result_func.set_coeffs(new_coefficients.tolist())
+        result_func = Function(len(new_coeffs) - 1)
+        result_func.set_coeffs(new_coeffs)
         return result_func
     
     def _strip_leading_zeros(self, coeffs: np.ndarray) -> np.ndarray:
@@ -1209,11 +1190,11 @@ class Function:
         self._check_initialized()
         other._check_initialized()
 
-        new_coefficients = np.polysub(self.coefficients, other.coefficients)
-        new_coefficients = self._strip_leading_zeros(new_coefficients)
+        new_poly = self.poly_obj - other.poly_obj
+        new_coeffs = new_poly.coef[::-1].tolist()
         
-        result_func = Function(len(new_coefficients) - 1)
-        result_func.set_coeffs(new_coefficients.tolist())
+        result_func = Function(len(new_coeffs) - 1)
+        result_func.set_coeffs(new_coeffs)
         return result_func
     
     def _multiply_by_scalar(self, scalar: Union[int, float, complex]) -> 'Function':
@@ -1237,15 +1218,11 @@ class Function:
         self._check_initialized()
         other._check_initialized()
 
-        # np.polymul performs convolution of coefficients to multiply polynomials
-        new_coefficients = np.polymul(self.coefficients, other.coefficients)
-        new_coefficients = self._strip_leading_zeros(new_coefficients)
+        new_poly = self.poly_obj * other.poly_obj
+        new_coeffs = new_poly.coef[::-1].tolist()
     
-        # The degree of the resulting polynomial is derived from the new coefficients
-        new_degree = len(new_coefficients) - 1
-    
-        result_func = Function(new_degree)
-        result_func.set_coeffs(new_coefficients.tolist())
+        result_func = Function(len(new_coeffs) - 1)
+        result_func.set_coeffs(new_coeffs)
         return result_func
         
     def __mul__(self, other: Union['Function', int, float, complex]) -> 'Function':
@@ -1269,16 +1246,18 @@ class Function:
     
         if isinstance(other, (int, float, complex)):
             if other == 0:
-                self.coefficients = np.array([0], dtype=self.coefficients.dtype)
-                self._largest_exponent = 0
+                self.set_coeffs([0])
             else:
-                self.coefficients *= other
+                # Update the underlying poly_obj, then sync the legacy array
+                self.poly_obj *= other
+                self.coefficients = self.poly_obj.coef[::-1]
+                self._largest_exponent = len(self.coefficients) - 1
             
         elif isinstance(other, self.__class__):
             other._check_initialized()
-            self.coefficients = np.polymul(self.coefficients, other.coefficients)
+            self.poly_obj *= other.poly_obj
+            self.coefficients = self.poly_obj.coef[::-1]
             self._largest_exponent = len(self.coefficients) - 1
-        
         else:
             return NotImplemented
         
@@ -1289,20 +1268,20 @@ class Function:
         Checks if two Function objects are equal by comparing
         their coefficients.
         """
-        # Check if the 'other' object is even a Function
         if not isinstance(other, Function):
             return NotImplemented
         
-        # Ensure both are initialized before trying to access .coefficients
         if not self._initialized or not other._initialized:
             return False
 
-        c1 = self._strip_leading_zeros(self.coefficients)
-        c2 = self._strip_leading_zeros(other.coefficients)
+        # Extract the auto-trimmed arrays from the modern API
+        c1 = self.poly_obj.coef
+        c2 = other.poly_obj.coef
         
         if c1.shape != c2.shape:
             return False
             
+        # Retain np.allclose to protect against floating-point inaccuracies
         return np.allclose(c1, c2)
 
 
@@ -1320,10 +1299,9 @@ class Function:
         a, b, c = self.coefficients
 
         discriminant = (b**2) - (4*a*c)
-
         sqrt_discriminant = cmath.sqrt(discriminant)
 
-        if b >= 0:
+        if b.real >= 0:
             sign_b = 1
         else:
             sign_b = -1

tests/test_polysolve.py

@@ -18,6 +18,13 @@ def quadratic_func() -> Function:
     f.set_coeffs([2, -3, -5])
     return f
 
+@pytest.fixture
+def int_func() -> Function:
+    """Provides a function for easy integral calculation."""
+    f = Function(largest_exponent=2)
+    f.set_coeffs([3, -4, -5])
+    return f
+
 @pytest.fixture
 def linear_func() -> Function:
     """Provides a standard linear function: x + 10."""
@@ -70,6 +77,20 @@ def test_nth_derivative(quadratic_func):
     # The derivative of 2x^2 - 3x - 5 is 4x - 3
     assert np.array_equal(derivative.coefficients, [4])
 
+def test_integral(int_func):
+    """Tests the calculation of the function's integral."""
+    integral = int_func.integral()
+    assert integral.largest_exponent == 3
+    # The integral of 2x^2 -3x - 5 is x^3 - 2x^2 - 5x
+    assert np.array_equal(integral.coefficients, [1, -2, -5, 0])
+
+def test_integral_custom_constant(int_func):
+    """Tests the calculation of the function's integral with custom constant of integration."""
+    integral = int_func.integral(constant_of_integration=42)
+    assert integral.largest_exponent == 3
+    # The integral of 2x^2 -3x - 5 is x^3 - 2x^2 - 5x
+    assert np.array_equal(integral.coefficients, [1, -2, -5, 42])
+
 def test_quadratic_solve(quadratic_func):
     """Tests the analytical quadratic solver for exact roots."""
     roots = quadratic_func.quadratic_solve()
@@ -139,7 +160,7 @@ def test_get_real_roots_numpy(quadratic_func):
     Tests that the NumPy-based genetic algorithm approximates the roots correctly.
     """
     # Using more generations for higher accuracy in testing
-    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
+    ga_opts = GA_Options(num_of_generations=100, data_size=500000, selection_percentile=0.66, root_precision=2)
     
     roots = quadratic_func.get_real_roots(ga_opts, use_cuda=False)
     
@@ -158,7 +179,7 @@ def test_get_real_roots_cuda(quadratic_func):
     It will be skipped automatically if CuPy is not available.
     """
     
-    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
+    ga_opts = GA_Options(num_of_generations=100, data_size=500000, selection_percentile=0.66, root_precision=2)
     
     roots = quadratic_func.get_real_roots(ga_opts, use_cuda=True)
     
@@ -172,7 +193,7 @@ def test_get_roots_numpy(complex_func):
     Tests that the NumPy-based genetic algorithm approximates the roots correctly.
     """
     # Using more generations for higher accuracy in testing
-    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
+    ga_opts = GA_Options(num_of_generations=100, data_size=500000, selection_percentile=0.66, root_precision=2)
     
     roots = complex_func.get_roots(ga_opts, use_cuda=False)
     
@@ -191,7 +212,7 @@ def test_get_roots_cuda(complex_func):
     It will be skipped automatically if CuPy is not available.
     """
     
-    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
+    ga_opts = GA_Options(num_of_generations=100, data_size=500000, selection_percentile=0.66, root_precision=2)
     
     roots = complex_func.get_roots(ga_opts, use_cuda=True)