diff --git a/src/__main__.py b/src/__main__.py
index d067201..9f31ca7 100644
--- a/src/__main__.py
+++ b/src/__main__.py
@@ -1,51 +1,5 @@
-from collections.abc import Iterator
-
-import numpy as np
-
-from src.function import Function
 from src.functions import available_functions
-from src.types import INPUT_VECTOR
-from src.utils import generate_bounded_points
-
-
-def search(
-    function: Function,
-    x0: INPUT_VECTOR,
-    iterations: int,
-    neighbors_count: int,
-    std_dev: float,
-    rng: np.random.Generator | None = None,
-) -> Iterator[INPUT_VECTOR]:
-    """Search for the minimum value of the function using the local search algorithm.
-
-    On each iteration, N neighboring points will be generated around the starting point. These points
-    will then be evaluated on the function, finding the smallest one. This smallest point will become
-    the new starting point, repeating until we run out of iterations.
-
-    Params:
-        x0: Starting point (N-dimensional).
-        iterations: Maximum number of iterations.
-        include_origin: When searching for the next minimum, should the origin point be checked too?
-        neighbors_count: The amount of neighbor points.
-        std_dev: Standard deviation for the normal distribution for neighbor generating.
-        rng: Random generator instance (None for a new rng).
-
-    Yields:
-        Minimum input vector (x) found so far, yielded from each iteration.
-    """
-    if rng is None:
-        rng = np.random.default_rng()
-    x_center = x0
-
-    for _ in range(iterations):
-        y_min = function(x_center)
-        for point in generate_bounded_points(x_center, neighbors_count, std_dev, function.definition_interval, rng):
-            y_point = function(point)
-            if y_point < y_min:
-                y_min = y_point
-                x_center = point
-
-        yield x_center
+from src.search import min_search
 
 
 def main() -> None:
@@ -56,8 +10,15 @@ def main() -> None:
     stddev = 1
 
     for func_name, function in available_functions.items():
+        print("-" * 80)
         print(func_name)
-        for x in search(function, function.definition_interval.random_point(dims), iters, neighbors, stddev):
+        for x in min_search(
+            function,
+            function.definition_interval.random_point(dims),
+            iterations=iters,
+            neighbors_count=neighbors,
+            std_dev=stddev,
+        ):
             y = function(x)
             print(f"{x} -> {y}")
 
diff --git a/src/search.py b/src/search.py
new file mode 100644
index 0000000..2a7a681
--- /dev/null
+++ b/src/search.py
@@ -0,0 +1,59 @@
+from collections.abc import Iterator
+
+import numpy as np
+
+from src.function import Function
+from src.types import INPUT_VECTOR
+from src.utils import generate_bounded_points
+
+
+def min_search(
+    function: Function,
+    x0: INPUT_VECTOR,
+    *,
+    iterations: int,
+    neighbors_count: int,
+    std_dev: float,
+    include_center: bool = True,
+    rng: np.random.Generator | None = None,
+) -> Iterator[INPUT_VECTOR]:
+    """Search for the minimum value of the function using the Local Search / Stochastic Hill Climber algorithm.
+
+    On each iteration, N neighboring points will be generated around the starting point. These points
+    will then be evaluated on the function, finding the smallest one. This smallest point will become
+    the new starting point, repeating until we run out of iterations.
+
+    Params:
+        x0: Starting point (N-dimensional).
+        iterations: Maximum number of iterations.
+        include_origin: When searching for the next minimum, should the origin point be checked too?
+        neighbors_count: The amount of neighbor points.
+        std_dev: Standard deviation for the normal distribution for neighbor generating.
+        include_center:
+            - When False, this algorithm works as Stochastic Hill Climber, always picking the minimum
+              from the newly generated points only, ignoring the current point.
+            - When True, this algorithm works as Local Search, picking the minimum from all of the
+              newly generated points and also from the center point.
+        rng: Random generator instance (None for a new rng).
+
+    Yields:
+        Minimum input vector (x) found so far, yielded from each iteration.
+    """
+    if rng is None:
+        rng = np.random.default_rng()
+    x_center = x0
+    y_min = function(x_center) if include_center else float("inf")
+
+    for _ in range(iterations):
+        for point in generate_bounded_points(x_center, neighbors_count, std_dev, function.definition_interval, rng):
+            y_point = function(point)
+            if y_point < y_min:
+                y_min = y_point
+                x_center = point
+
+        # If we're using Stochastic Hill Climber, reset the minimum value, to only
+        # pick it from the newly generated points, not the current center
+        if not include_center:
+            y_min = float("inf")
+
+        yield x_center