diff --git a/pygad/utils/crossover.py b/pygad/utils/crossover.py
index 4a85f86c..0da8bd53 100644
--- a/pygad/utils/crossover.py
+++ b/pygad/utils/crossover.py
@@ -114,17 +114,22 @@ def two_points_crossover(self, parents, offspring_size):
         else:
             offspring = numpy.empty(offspring_size, dtype=object)
 
-        # Randomly generate all the first K points at which crossover takes place between each two parents. 
-        # This saves time by calling the numpy.random.randint() function only once.
-        if (parents.shape[1] == 1): # If the chromosome has only a single gene. In this case, this gene is copied from the second parent.
-            crossover_points_1 = numpy.zeros(offspring_size[0])
-        else:
-            crossover_points_1 = numpy.random.randint(low=0, 
-                                                      high=numpy.ceil(parents.shape[1]/2 + 1), 
-                                                      size=offspring_size[0])
+        # Randomly generate all the K pairs of points at which crossover takes place between each two parents.
+        # This saves time by calling the numpy.random.randint() function only twice.
+        # The 2 points of a pair are different values in [0, num_genes], and every such pair is equally likely.
+        # If the chromosome has only a single gene, the points are 0 and 1: the gene is copied from the second parent.
+        points_a = numpy.random.randint(low=0,
+                                        high=parents.shape[1] + 1,
+                                        size=offspring_size[0])
+        points_b = numpy.random.randint(low=0,
+                                        high=parents.shape[1],
+                                        size=offspring_size[0])
+        # Skip the value of the first point so that the 2 points differ.
+        points_b[points_b >= points_a] += 1
 
         # The second point must always be greater than the first point.
-        crossover_points_2 = crossover_points_1 + int(parents.shape[1]/2) 
+        crossover_points_1 = numpy.minimum(points_a, points_b)
+        crossover_points_2 = numpy.maximum(points_a, points_b)
 
         for k in range(offspring_size[0]):
 
diff --git a/tests/test_crossover_mutation.py b/tests/test_crossover_mutation.py
index acc38942..497376f1 100644
--- a/tests/test_crossover_mutation.py
+++ b/tests/test_crossover_mutation.py
@@ -241,6 +241,32 @@ def test_random_mutation_manual_call4():
     for value in comp_sorted:
         assert value in value_space
 
+def test_two_points_crossover_manual_call():
+    # Both points are random: the genes between them form 1 segment of any length from 1 to num_genes.
+    num_genes = 10
+    num_offspring = 1000
+    result, ga_instance = output_crossover_mutation(gene_type=int,
+                                                    crossover_type="two_points")
+
+    parents = numpy.array([[0] * num_genes,
+                           [1] * num_genes])
+    offspring = ga_instance.two_points_crossover(parents=parents,
+                                                 offspring_size=(num_offspring, num_genes))
+
+    # Without crossover_probability, the first parent of offspring k is parents[k % 2].
+    # Mark the genes that come from the second parent with 1.
+    from_second_parent = numpy.array([child if k % 2 == 0 else 1 - child for k, child in enumerate(offspring)])
+
+    segment_lengths = set()
+    for child in from_second_parent:
+        segment = numpy.flatnonzero(child)
+        # The genes from the second parent are consecutive.
+        assert len(segment) > 0
+        assert segment[-1] - segment[0] + 1 == len(segment)
+        segment_lengths.add(len(segment))
+
+    assert segment_lengths == set(range(1, num_genes + 1))
+
 if __name__ == "__main__":
     #### Single-objective
     print()
@@ -285,3 +311,6 @@ def test_random_mutation_manual_call4():
 
     test_random_mutation_manual_call4()
     print()
+
+    test_two_points_crossover_manual_call()
+    print()