v0.7.0 - Complex Number Support (#25)

Reviewed-on: #25
Co-authored-by: Jonathan Rampersad <rampersad.jonathan@gmail.com>
Co-committed-by: Jonathan Rampersad <rampersad.jonathan@gmail.com>

README.md

@@ -12,6 +12,7 @@ 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.
@@ -75,12 +76,19 @@ roots_analytic = f1.quadratic_solve()
 print(f"Analytic roots: {sorted(roots_analytic)}")
 # > Analytic roots: [-1.0, 2.5]
 
-# 6. Find roots with the genetic algorithm (Numba CPU)
-#    This is the default, JIT-compiled CPU solver.
+# 6. 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 roots from GA: {roots_ga[:2]}")
-# > Approximate roots from GA: [-1.000..., 2.500...]
+print(f"Approximate real roots: {roots_ga[:2]}")
+# > Approximate real roots: [-1.000..., 2.500...]
+
+# 7. 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]
 
 # 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)
@@ -114,7 +122,10 @@ 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
+    blend_alpha=0.75,
+
+    # Enable complex root finding (default is True)
+    find_complex=True
 )
 
 # Pass the custom options to the solver

src/polysolve/__init__.py

@@ -343,7 +343,12 @@ class GA_Options:
                 stacklevel=2
             )
         if self.min_range != 0.0 or self.max_range != 0.0:
-            warnings.warn("min_range and max_range are no longer used, instead cauchy's bound is used to find these values")
+            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
+            )
 
 def _get_cauchy_bound(coeffs: np.ndarray) -> float:
     """
@@ -688,7 +693,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 + noise
+            dst_solutions[idx_cross_end:idx_mut_end] = mutation_candidates * (1.0 + 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)
@@ -800,7 +805,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 + noise_real) + 1j * (mut_candidates.imag + noise_imag)
+            dst_solutions[idx_cross_end:idx_mut_end] = (mut_candidates.real * (1.0 + noise_real)) + 1j * (mut_candidates.imag * (1.0 + noise_imag))
 
             # 4. New Randoms: Add new blood to prevent getting stuck
             rand_real = np.random.uniform(min_r, max_r, random_size)
@@ -930,7 +935,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 + noise
+            d_dst_solutions[idx_cross_end:idx_mut_end] = d_mutation_candidates * (1.0 + noise)
 
             # 4. New Randoms
             d_dst_solutions[idx_mut_end:] = cupy.random.uniform(
@@ -1095,7 +1100,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 + noise_real) + 1j * (mut_candidates.imag + noise_imag)
+            d_dst_solutions[idx_cross_end:idx_mut_end] = (mut_candidates.real * (1.0 + noise_real)) + 1j * (mut_candidates.imag * (1.0 + noise_imag))
 
             # 4. Random Injection: Fresh genetic material
             rand_real = cupy.random.uniform(min_r, max_r, random_size, dtype=cupy.float64)
@@ -1127,7 +1132,9 @@ class Function:
 
         d_unique = cupy.unique(rounded_real + 1j * rounded_imag)
         
-        return cupy.asnumpy(d_unique)
+        # Sort the unique roots and copy back to CPU
+        final_solutions_gpu = cupy.sort(d_unique)
+        return final_solutions_gpu.get()
 
 
     def __str__(self) -> str:
@@ -1299,7 +1306,7 @@ class Function:
         return np.allclose(c1, c2)
 
 
-    def quadratic_solve(self) -> Optional[List[complex]]:
+    def quadratic_solve(self) -> Optional[List[Union[complex, float]]]:
         """
         Calculates the roots (real or complex) of a quadratic function.
 
@@ -1328,9 +1335,14 @@ class Function:
         else:
             # Standard case: Use Vieta's formula
             root2 = (c / a) / root1
+        
+        roots = np.array([root1, root2])
+        roots.sort()
 
-        # Return roots in a consistent order
-        return [root1, root2]
+        if np.all(np.abs(roots.imag) < 1e-15):
+            return roots.real.astype(np.float64)
+    
+        return roots
 
 # Example Usage
 if __name__ == '__main__':