diff --git a/src/polysolve/__init__.py b/src/polysolve/__init__.py
index 5a2813b..b37b691 100644
--- a/src/polysolve/__init__.py
+++ b/src/polysolve/__init__.py
@@ -310,15 +310,17 @@ class Function:
             parent_pool = solutions[:data_size // 2]
 
             # 2. Crossover: Breed two parents to create a child
-            parent1_indices = np.random.randint(0, parent_pool.size, crossover_size)
-            parent2_indices = np.random.randint(0, parent_pool.size, crossover_size)
-            parents1 = parent_pool[parent1_indices]
-            parents2 = parent_pool[parent2_indices]
+            # Select from the full list (indices 0 to data_size-1)
+            parent1_indices = np.random.randint(0, data_size, crossover_size)
+            parent2_indices = np.random.randint(0, data_size, crossover_size)
+            parents1 = solutions[parent1_indices]
+            parents2 = solutions[parent2_indices]
             # Simple "average" crossover
             crossover_solutions = (parents1 + parents2) / 2.0
 
-            # 3. Mutation: Apply your original mutation logic
-            mutation_candidates = parent_pool[np.random.randint(0, parent_pool.size, mutation_size)]
+            # 3. Mutation:
+            # Select from the full list (indices 0 to data_size-1)
+            mutation_candidates = solutions[np.random.randint(0, data_size, mutation_size)]
             
             # Use mutation_strength (the new name)
             mutation_factors = np.random.uniform(
@@ -400,15 +402,17 @@ class Function:
             parent_pool_size = d_parent_pool.size
 
             # 2. Crossover
-            parent1_indices = cupy.random.randint(0, parent_pool_size, crossover_size)
-            parent2_indices = cupy.random.randint(0, parent_pool_size, crossover_size)
-            d_parents1 = d_parent_pool[parent1_indices]
-            d_parents2 = d_parent_pool[parent2_indices]
+            # Select from the full list (indices 0 to data_size-1)
+            parent1_indices = cupy.random.randint(0, data_size, crossover_size)
+            parent2_indices = cupy.random.randint(0, data_size, crossover_size)
+            d_parents1 = d_solutions[parent1_indices]
+            d_parents2 = d_solutions[parent2_indices]
             d_crossover_solutions = (d_parents1 + d_parents2) / 2.0
 
             # 3. Mutation
-            mutation_indices = cupy.random.randint(0, parent_pool_size, mutation_size)
-            d_mutation_candidates = d_parent_pool[mutation_indices]
+            # Select from the full list (indices 0 to data_size-1)
+            mutation_indices = cupy.random.randint(0, data_size, mutation_size)
+            d_mutation_candidates = d_solutions[mutation_indices]
             
             # Use mutation_strength (the new name)
             d_mutation_factors = cupy.random.uniform(