feat: add complex root finding and dynamic CUDA shared memory optimization

Major update extending the library to solve for complex roots and optimizing GPU performance using Shared Memory.

Complex Number Support:
- Implemented `_solve_complex_cuda` and `_solve_complex_numpy` to find roots in the complex plane.
- Added specialized CUDA kernels (`_FITNESS_KERNEL_COMPLEX`, `_FITNESS_KERNEL_COMPLEX_DYNAMIC`) handling complex arithmetic (multiplication/addition) directly on the GPU.
- Updated `Function` class and `set_coeffs` to handle `np.complex128` data types.
- Updated `quadratic_solve` to return complex roots using `cmath`.

CUDA Performance & Optimization:
- Implemented Dynamic Shared Memory kernels (`extern __shared__`) to cache polynomial coefficients on the GPU block, significantly reducing global memory latency.
- Added intelligent fallback logic: The solver checks `MaxSharedMemoryPerBlock`. If the polynomial is too large for Shared Memory, it falls back to the standard Global Memory kernel to prevent crashes.
- Split complex coefficients into separate Real and Imaginary arrays for CUDA kernel efficiency.

Polynomial Logic:
- Added `_strip_leading_zeros` helper to ensure polynomial degree is correctly maintained after arithmetic operations (e.g., preventing `0x^2 + x` from being treated as degree 2).
- Updated `__init__` to allow direct coefficient injection.

GA Algorithm:
- Updated crossover logic to support 2D search space (Real + Imaginary) for complex solutions.
- Refined fitness function to explicitly handle `isinf`/`isnan` for numerical stability.

.gitea/workflows/publish-python.yml

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

.gitea/workflows/test.yml

@@ -6,7 +6,7 @@ jobs:
     runs-on: nvidia-cuda
     strategy:
       matrix:
-        python-version: ["3.10", "3.12", "3.14"]
+        python-version: ["3.8", "3.10", "3.12"]
 
     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, 3.14
+* **Python Versions Tested:** 3.8, 3.10, 3.12
 * **CI Docker Image:** You can pull the exact image used by our runners from Docker Hub:
 
 ```bash

README.md

@@ -12,7 +12,6 @@ A Python library for representing, manipulating, and solving polynomial equation
 ## Key Features
 
 * **Numerically Stable Solver**: Makes complex calculations **practical**. Leverage your GPU to power the robust genetic algorithm, solving high-degree polynomials accurately in a reasonable timeframe.
-* **Complex Number Support**: Fully supports complex coefficients and finding roots in the complex plane (e.g., $x^2 + 1 = 0 \to \pm i$).
 * **Numba Accelerated CPU Solver**: The default genetic algorithm is JIT-compiled with Numba for high-speed CPU performance, right out of the box.
 * **CUDA Accelerated**: Leverage NVIDIA GPUs for a massive performance boost when finding roots in large solution spaces.
 * **Create and Manipulate Polynomials**: Easily define polynomials of any degree using integer or float coefficients, and perform arithmetic operations like addition, subtraction, multiplication, and scaling.
@@ -51,8 +50,6 @@ 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
@@ -72,35 +69,18 @@ ddf1 = f1.nth_derivative(2)
 print(f"2nd Derivative of f1: {ddf1}")
 # > Derivative of f1: 4
 
-# 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
+# 5. 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]
 
-# 8. Find REAL roots with the genetic algorithm (Numba CPU)
-#    This is the default, JIT-compiled CPU solver.
+# 6. Find 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...]
-
-# 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)
-print(f"Approximate complex roots: {roots_all}")
-# > Approximate complex roots: [-1.00...j, 1.00...j]
+print(f"Approximate roots from GA: {roots_ga[:2]}")
+# > Approximate roots from GA: [-1.000..., 2.500...]
 
 # If you installed a CUDA extra, you can run it on the GPU:
 # roots_ga_gpu = f1.get_real_roots(ga_opts, use_cuda=True)
@@ -134,10 +114,7 @@ ga_robust_search = GA_Options(
     # Increase the crossover blend factor to 0.75.
     # This allows new solutions to be created further
     # away from their parents, increasing exploration.
-    blend_alpha=0.75,
-
-    # Enable complex root finding (default is True)
-    find_complex=True
+    blend_alpha=0.75
 )
 
 # Pass the custom options to the solver

pyproject.toml

@@ -5,13 +5,13 @@ build-backend = "setuptools.build_meta"
 [project]
 # --- Core Metadata ---
 name = "polysolve"
-version = "0.8.1"
+version = "0.7.0"
 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.10"
+requires-python = ">=3.8"
 license = { file="LICENSE" }
 keywords = ["math", "polynomial", "genetic algorithm", "cuda", "equation solver"]
 
@@ -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>=2.0",
-    "numba>=0.60"
+    "numpy>=1.21",
+    "numba"
 ]
 
 # --- Optional Dependencies (Extras) ---
 [project.optional-dependencies]
-cuda12 = ["cupy-cuda12x[ctk]"]
+cuda12 = ["cupy-cuda12x"]
 dev = ["pytest"]
 
 [project.urls]

src/polysolve/__init__.py

@@ -1,7 +1,6 @@
 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
@@ -344,12 +343,7 @@ class GA_Options:
                 stacklevel=2
             )
         if self.min_range != 0.0 or self.max_range != 0.0:
-            warnings.warn(
-                "The 'min_range' and 'max_range' parameters are deprecated and will be ignored. "
-                "Search bounds are now automatically calculated using Cauchy's bound.",
-                DeprecationWarning,
-                stacklevel=2
-            )
+            warnings.warn("min_range and max_range are no longer used, instead cauchy's bound is used to find these values")
 
 def _get_cauchy_bound(coeffs: np.ndarray) -> float:
     """
@@ -384,9 +378,6 @@ 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
@@ -394,6 +385,7 @@ 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.")
 
@@ -428,9 +420,6 @@ 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
 
@@ -460,7 +449,7 @@ class Function:
             float: The resulting y-value.
         """
         self._check_initialized()
-        return self.poly_obj(x_val)
+        return np.polyval(self.coefficients, x_val)
 
     def differential(self) -> 'Function':
         """
@@ -495,11 +484,10 @@ class Function:
             diff_func.set_coeffs([0])
             return diff_func
         
-        deriv_poly = self.poly_obj.deriv(1)
-        new_coeffs = deriv_poly.coef[::-1].tolist()
+        derivative_coefficients = np.polyder(self.coefficients)
         
         diff_func = Function(self._largest_exponent - 1)
-        diff_func.set_coeffs(new_coeffs)
+        diff_func.set_coeffs(derivative_coefficients.tolist())
         return diff_func
     
 
@@ -519,30 +507,18 @@ class Function:
             raise ValueError("Derivative order 'n' must be a positive integer.")
 
         if n > self.largest_exponent:
-            diff_func = Function(0)
-            diff_func.set_coeffs([0])
-            return diff_func
+            function = Function(0)
+            function.set_coeffs([0])
+            return function
 
-        deriv_poly = self.poly_obj.deriv(n)
-        new_coeffs = deriv_poly.coef[::-1].tolist()
+        if n == 1:
+            return self.derivative()
         
-        new_degree = self._largest_exponent - n
-        diff_func = Function(new_degree)
-        diff_func.set_coeffs(new_coeffs)
-        return diff_func
+        function = self
+        for _ in range(n):
+            function = function.derivative()
 
-
-    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
+        return function
 
 
     def get_real_roots(self, options: Optional[GA_Options] = None, use_cuda: bool = False) -> np.ndarray:
@@ -712,7 +688,7 @@ class Function:
             
             # Use mutation_strength
             noise = np.random.normal(0, options.mutation_strength, mutation_size)
-            dst_solutions[idx_cross_end:idx_mut_end] = mutation_candidates * (1.0 + noise)
+            dst_solutions[idx_cross_end:idx_mut_end] = mutation_candidates + noise
 
             # 4. New Randoms: Add new blood to prevent getting stuck
             dst_solutions[idx_mut_end:] = np.random.uniform(min_r, max_r, random_size)
@@ -824,7 +800,7 @@ class Function:
             noise_real = np.random.normal(0, options.mutation_strength, mutation_size)
             noise_imag = np.random.normal(0, options.mutation_strength, mutation_size)
             
-            dst_solutions[idx_cross_end:idx_mut_end] = (mut_candidates.real * (1.0 + noise_real)) + 1j * (mut_candidates.imag * (1.0 + noise_imag))
+            dst_solutions[idx_cross_end:idx_mut_end] = (mut_candidates.real + noise_real) + 1j * (mut_candidates.imag + noise_imag)
 
             # 4. New Randoms: Add new blood to prevent getting stuck
             rand_real = np.random.uniform(min_r, max_r, random_size)
@@ -954,7 +930,7 @@ class Function:
             
             # Use mutation_strength
             noise = cupy.random.normal(0, options.mutation_strength, mutation_size)
-            d_dst_solutions[idx_cross_end:idx_mut_end] = d_mutation_candidates * (1.0 + noise)
+            d_dst_solutions[idx_cross_end:idx_mut_end] = d_mutation_candidates + noise
 
             # 4. New Randoms
             d_dst_solutions[idx_mut_end:] = cupy.random.uniform(
@@ -1119,7 +1095,7 @@ class Function:
             noise_imag = cupy.random.normal(0, options.mutation_strength, mutation_size)
 
             #    Apply Mutation: Scale Real/Imag independently to allow "rotation" off the line
-            d_dst_solutions[idx_cross_end:idx_mut_end] = (mut_candidates.real * (1.0 + noise_real)) + 1j * (mut_candidates.imag * (1.0 + noise_imag))
+            d_dst_solutions[idx_cross_end:idx_mut_end] = (mut_candidates.real + noise_real) + 1j * (mut_candidates.imag + noise_imag)
 
             # 4. Random Injection: Fresh genetic material
             rand_real = cupy.random.uniform(min_r, max_r, random_size, dtype=cupy.float64)
@@ -1151,17 +1127,54 @@ class Function:
 
         d_unique = cupy.unique(rounded_real + 1j * rounded_imag)
         
-        # Sort the unique roots and copy back to CPU
-        final_solutions_gpu = cupy.sort(d_unique)
-        return final_solutions_gpu.get()
+        return cupy.asnumpy(d_unique)
 
 
     def __str__(self) -> str:
         """Returns a human-readable string representation of the function."""
-        if not self._initialized:
-            return "Uninitialized Function"
+        self._check_initialized()
+        parts = []
+        for i, c in enumerate(self.coefficients):
+            if c == 0:
+                continue
 
-        return str(self.poly_obj)
+            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"
 
     def __repr__(self) -> str:
         return f"Function(str='{self}')"
@@ -1171,12 +1184,11 @@ class Function:
         self._check_initialized()
         other._check_initialized()
 
-        # 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()
+        new_coefficients = np.polyadd(self.coefficients, other.coefficients)
+        new_coefficients = self._strip_leading_zeros(new_coefficients)
         
-        result_func = Function(len(new_coeffs) - 1)
-        result_func.set_coeffs(new_coeffs)
+        result_func = Function(len(new_coefficients) - 1)
+        result_func.set_coeffs(new_coefficients.tolist())
         return result_func
     
     def _strip_leading_zeros(self, coeffs: np.ndarray) -> np.ndarray:
@@ -1190,11 +1202,11 @@ class Function:
         self._check_initialized()
         other._check_initialized()
 
-        new_poly = self.poly_obj - other.poly_obj
-        new_coeffs = new_poly.coef[::-1].tolist()
+        new_coefficients = np.polysub(self.coefficients, other.coefficients)
+        new_coefficients = self._strip_leading_zeros(new_coefficients)
         
-        result_func = Function(len(new_coeffs) - 1)
-        result_func.set_coeffs(new_coeffs)
+        result_func = Function(len(new_coefficients) - 1)
+        result_func.set_coeffs(new_coefficients.tolist())
         return result_func
     
     def _multiply_by_scalar(self, scalar: Union[int, float, complex]) -> 'Function':
@@ -1218,11 +1230,15 @@ class Function:
         self._check_initialized()
         other._check_initialized()
 
-        new_poly = self.poly_obj * other.poly_obj
-        new_coeffs = new_poly.coef[::-1].tolist()
+        # 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)
     
-        result_func = Function(len(new_coeffs) - 1)
-        result_func.set_coeffs(new_coeffs)
+        # 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())
         return result_func
         
     def __mul__(self, other: Union['Function', int, float, complex]) -> 'Function':
@@ -1246,18 +1262,16 @@ class Function:
     
         if isinstance(other, (int, float, complex)):
             if other == 0:
-                self.set_coeffs([0])
+                self.coefficients = np.array([0], dtype=self.coefficients.dtype)
+                self._largest_exponent = 0
             else:
-                # 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
+                self.coefficients *= other
             
         elif isinstance(other, self.__class__):
             other._check_initialized()
-            self.poly_obj *= other.poly_obj
-            self.coefficients = self.poly_obj.coef[::-1]
+            self.coefficients = np.polymul(self.coefficients, other.coefficients)
             self._largest_exponent = len(self.coefficients) - 1
+        
         else:
             return NotImplemented
         
@@ -1268,24 +1282,24 @@ 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
 
-        # Extract the auto-trimmed arrays from the modern API
-        c1 = self.poly_obj.coef
-        c2 = other.poly_obj.coef
+        c1 = self._strip_leading_zeros(self.coefficients)
+        c2 = self._strip_leading_zeros(other.coefficients)
         
         if c1.shape != c2.shape:
             return False
             
-        # Retain np.allclose to protect against floating-point inaccuracies
         return np.allclose(c1, c2)
 
 
-    def quadratic_solve(self) -> Optional[List[Union[complex, float]]]:
+    def quadratic_solve(self) -> Optional[List[complex]]:
         """
         Calculates the roots (real or complex) of a quadratic function.
 
@@ -1299,9 +1313,10 @@ class Function:
         a, b, c = self.coefficients
 
         discriminant = (b**2) - (4*a*c)
+
         sqrt_discriminant = cmath.sqrt(discriminant)
 
-        if b.real >= 0:
+        if b >= 0:
             sign_b = 1
         else:
             sign_b = -1
@@ -1314,13 +1329,8 @@ class Function:
             # Standard case: Use Vieta's formula
             root2 = (c / a) / root1
 
-        roots = np.array([root1, root2])
-        roots.sort()
-
-        if np.all(np.abs(roots.imag) < 1e-15):
-            return roots.real.astype(np.float64)
-    
-        return roots
+        # Return roots in a consistent order
+        return [root1, root2]
 
 # Example Usage
 if __name__ == '__main__':

tests/test_polysolve.py

@@ -18,13 +18,6 @@ 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."""
@@ -77,20 +70,6 @@ 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()
@@ -160,7 +139,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=100, data_size=500000, selection_percentile=0.66, root_precision=2)
+    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
     
     roots = quadratic_func.get_real_roots(ga_opts, use_cuda=False)
     
@@ -179,7 +158,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=100, data_size=500000, selection_percentile=0.66, root_precision=2)
+    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
     
     roots = quadratic_func.get_real_roots(ga_opts, use_cuda=True)
     
@@ -193,7 +172,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=100, data_size=500000, selection_percentile=0.66, root_precision=2)
+    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
     
     roots = complex_func.get_roots(ga_opts, use_cuda=False)
     
@@ -212,7 +191,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=100, data_size=500000, selection_percentile=0.66, root_precision=2)
+    ga_opts = GA_Options(num_of_generations=50, data_size=200000, selection_percentile=0.66, root_precision=3)
     
     roots = complex_func.get_roots(ga_opts, use_cuda=True)