diff --git a/src/qha/fitting.py b/src/qha/fitting.py
index 163dc00..a20b7f5 100644
--- a/src/qha/fitting.py
+++ b/src/qha/fitting.py
@@ -17,7 +17,10 @@
 __all__ = ["polynomial_least_square_fitting", "apply_finite_strain_fitting"]
 
 
-def polynomial_least_square_fitting(xs, ys, new_xs, order: Optional[int] = 3):
+@jit(float64[:](float64[:], float64[:], float64[:], int64), nopython=True, cache=True)
+def polynomial_least_square_fitting(
+    xs: Vector, ys: Vector, new_xs: Vector, order: Optional[int] = 3
+):
     """
     The algorithm is referenced from the
     `Wolfram MathWorld <http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html>`_.
@@ -33,9 +36,9 @@ def polynomial_least_square_fitting(xs, ys, new_xs, order: Optional[int] = 3):
     xx = np.vander(
         xs, order, increasing=True
     )  # This will make a Vandermonde matrix that will be used in EoS fitting.
-    a, _, _, _ = np.linalg.lstsq(xx, ys, rcond=None)
+    a, _, _, _ = np.linalg.lstsq(xx, ys)  # See https://stackoverflow.com/a/64224087
     new_y = np.vander(new_xs, order, increasing=True) @ a
-    return a, new_y
+    return new_y
 
 
 @jit(
@@ -69,8 +72,7 @@ def apply_finite_strain_fitting(
     )  # Initialize the F(T,V) array
 
     for i in range(temperature_amount):
-        _, f_i = polynomial_least_square_fitting(
+        f_v_t[i] = polynomial_least_square_fitting(
             strains_sparse, free_energies[i], strains_dense, order
         )
-        f_v_t[i] = f_i
     return f_v_t
diff --git a/src/qha/grid_interpolation.py b/src/qha/grid_interpolation.py
index bab89a3..2bdf35c 100644
--- a/src/qha/grid_interpolation.py
+++ b/src/qha/grid_interpolation.py
@@ -203,7 +203,7 @@ def approach_to_best_ratio(
         vr.interpolate_volumes()
         strains, finer_volumes = vr.strains, vr.out_volumes
         eulerian_strain = calculate_eulerian_strain(volumes[0], volumes)
-        _, f_v_tmax = polynomial_least_square_fitting(
+        f_v_tmax = polynomial_least_square_fitting(
             eulerian_strain, free_energies, strains, self.option
         )
         p_v_tmax = -np.gradient(f_v_tmax) / np.gradient(finer_volumes)
diff --git a/src/qha/tools.py b/src/qha/tools.py
index a35c61c..16c9051 100644
--- a/src/qha/tools.py
+++ b/src/qha/tools.py
@@ -122,7 +122,7 @@ def calibrate_energy_on_reference(
         strains_after_calibration = calculate_eulerian_strain(
             volumes_before_calibration[i, 0], volumes_for_reference
         )
-        _, energies_after_calibration[i, :] = polynomial_least_square_fitting(
+        energies_after_calibration[i, :] = polynomial_least_square_fitting(
             strains_before_calibration,
             energies_before_calibration[i],
             strains_after_calibration,