diff --git a/.gitea/workflows/test.yml b/.gitea/workflows/test.yml index f0ba760..4541f68 100644 --- a/.gitea/workflows/test.yml +++ b/.gitea/workflows/test.yml @@ -6,7 +6,7 @@ jobs: runs-on: nvidia-cuda strategy: matrix: - python-version: ["3.8", "3.10", "3.12", "3.14"] + python-version: ["3.10", "3.12", "3.14"] steps: - name: Checkout repository diff --git a/README.md b/README.md index 909e139..056a21f 100644 --- a/README.md +++ b/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) diff --git a/pyproject.toml b/pyproject.toml index 4985738..b2caa7c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,14 +5,14 @@ build-backend = "setuptools.build_meta" [project] # --- Core Metadata --- name = "polysolve" -version = "0.7.1" +version = "0.8.0" authors = [ { name="Jonathan Rampersad", email="jonathan@jono-rams.work" }, ] 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,8 +23,6 @@ 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", @@ -35,8 +33,8 @@ classifiers = [ # --- Dependencies --- dependencies = [ - "numpy>=1.21", - "numba" + "numpy>=2.0", + "numba>=0.60" ] # --- Optional Dependencies (Extras) --- @@ -47,5 +45,5 @@ dev = ["pytest"] [project.urls] Homepage = "https://polysolve.jono-rams.work" Documentation = "https://polysolve.jono-rams.work/docs" -Repository = "https://github.com/jono-rams/PolySolve" +Repository = "https://github.com/jono-rams/PolySolve" "Bug Tracker" = "https://github.com/jono-rams/PolySolve/issues" diff --git a/src/polysolve/__init__.py b/src/polysolve/__init__.py index 917a8b5..29c89ca 100644 --- a/src/polysolve/__init__.py +++ b/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 diff --git a/tests/test_polysolve.py b/tests/test_polysolve.py index fb54642..f09a475 100644 --- a/tests/test_polysolve.py +++ b/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()