diff --git a/opfunu/cec_based/cec.py b/opfunu/cec_based/cec.py
index 082f914..47cb881 100644
--- a/opfunu/cec_based/cec.py
+++ b/opfunu/cec_based/cec.py
@@ -173,7 +173,7 @@ def check_solution(self, x, dim_max=None, dim_support=None):
         """
         # if not self.dim_changeable and (len(x) != self._ndim):
         if len(x) != self._ndim:
-            raise ValueError(f"{self.__class__.__name__} problem, the length of solution should has {self._ndim} variables!")
+            raise ValueError(f"{self.__class__.__name__} problem, the length of solution should have {self._ndim} variables!")
         if (dim_max is not None) and (len(x) > dim_max):
             raise ValueError(f"{self.__class__.__name__} problem is not supported ndim > {dim_max}!")
         if (dim_support is not None) and (len(x) not in dim_support):
diff --git a/opfunu/cec_based/cec2019.py b/opfunu/cec_based/cec2019.py
index fdf7437..08048d0 100644
--- a/opfunu/cec_based/cec2019.py
+++ b/opfunu/cec_based/cec2019.py
@@ -49,13 +49,13 @@ def __init__(self, ndim=None, bounds=None, f_shift="shift_data_1", f_bias=1.):
         self.f_shift = self.check_shift_data(f_shift)[:self.ndim]
         self.f_bias = f_bias
         self.f_global = f_bias
-        self.x_global = self.f_shift
+        self.x_global = np.zeros(self.ndim)
         self.paras = {"f_shift": self.f_shift, "f_bias": self.f_bias}
 
     def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
-        return operator.storn_chebyshev_polynomial_fitting_func(x) + self.f_bias
+        return operator.chebyshev_func(x) + self.f_bias
 
 
 class F22019(CecBenchmark):
@@ -97,20 +97,23 @@ def __init__(self, ndim=None, bounds=None, f_shift="shift_data_2", f_bias=1.):
         self.make_support_data_path("data_2019")
         self.f_shift = self.check_shift_data(f_shift)[:self.ndim]
         self.f_bias = f_bias
-        self.f_global = f_bias
-        self.x_global = self.f_shift
+        # the f_global and x_global was obtained by executing the cec2019 c code
+        self.f_global = int(np.sqrt(self.ndim)) + f_bias
+        self.x_global = np.zeros(self.ndim)
         self.paras = {"f_shift": self.f_shift, "f_bias": self.f_bias}
 
     def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
-        return operator.inverse_hilbert_matrix_func(x) + self.f_bias
+        return operator.inverse_hilbert_func(x) + self.f_bias
 
 
 class F32019(CecBenchmark):
     """
     .. [1] The 100-Digit Challenge: Problem Definitions and Evaluation Criteria for the 100-Digit
     Challenge Special Session and Competition on Single Objective Numerical Optimization
+
+    **Note: The CEC 2019 implementation and this implementation results match when x* = [0,...,0] and
     """
     name = "F3: Lennard-Jones Minimum Energy Cluster Problem"
     latex_formula = r'F_1(x) = \sum_{i=1}^D z_i^2 + bias, z=x-o,\\ x=[x_1, ..., x_D]; o=[o_1, ..., o_D]: \text{the shifted global optimum}'
@@ -146,14 +149,15 @@ def __init__(self, ndim=None, bounds=None, f_shift="shift_data_3", f_bias=1.):
         self.make_support_data_path("data_2019")
         self.f_shift = self.check_shift_data(f_shift)[:self.ndim]
         self.f_bias = f_bias
-        self.f_global = f_bias
+        # f_global calculated by verifying the cec2019 C value for f(x*) where x*==f_shift
+        self.f_global = 12.712062001703194 + self.f_bias
         self.x_global = self.f_shift
         self.paras = {"f_shift": self.f_shift, "f_bias": self.f_bias}
 
     def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
-        return operator.lennard_jones_minimum_energy_cluster_func(x) + self.f_bias
+        return operator.lennard_jones_func(x) + self.f_bias
 
 
 class F42019(CecBenchmark):
@@ -257,7 +261,7 @@ def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
         z = np.dot(self.f_matrix, x - self.f_shift)
-        return operator.weierstrass_func(z) + self.f_bias
+        return operator.weierstrass_norm_func(z) + self.f_bias
 
 
 class F72019(F42019):
@@ -335,7 +339,7 @@ def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
         z = np.dot(self.f_matrix, x - self.f_shift)
-        return operator.happy_cat_func(z) + self.f_bias
+        return operator.happy_cat_func(z, shift=-1.0) + self.f_bias
 
 
 class F102019(F42019):
diff --git a/opfunu/cec_based/cec2020.py b/opfunu/cec_based/cec2020.py
index c9edfc3..99e8f25 100644
--- a/opfunu/cec_based/cec2020.py
+++ b/opfunu/cec_based/cec2020.py
@@ -110,7 +110,7 @@ def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
         z = np.dot(self.f_matrix, 1000.*(x - self.f_shift)/100)
-        return operator.bent_cigar_func(z) + self.f_bias
+        return operator.modified_schwefel_func(z) + self.f_bias
 
 
 class F32020(CecBenchmark):
@@ -162,7 +162,7 @@ def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
         z = np.dot(self.f_matrix, 600.*(x - self.f_shift)/100)
-        return operator.bent_cigar_func(z) + self.f_bias
+        return operator.lunacek_bi_rastrigin_func(z, shift=2.5) + self.f_bias
 
 
 class F42020(CecBenchmark):
@@ -170,7 +170,7 @@ class F42020(CecBenchmark):
     .. [1] Problem Definitions and Evaluation Criteria for the CEC 2020
     Special Session and Competition on Single Objective Bound Constrained Numerical Optimization
     """
-    name = "F4: Expanded Rosenbrock’s plus Griewangk’s Function (F15 CEC-2014)"
+    name = "F4: Expanded Rosenbrock’s plus Griewank’s Function (F15 CEC-2014)"
     latex_formula = r'F_1(x) = \sum_{i=1}^D z_i^2 + bias, z=x-o,\\ x=[x_1, ..., x_D]; o=[o_1, ..., o_D]: \text{the shifted global optimum}'
     latex_formula_dimension = r'2 <= D <= 100'
     latex_formula_bounds = r'x_i \in [-100.0, 100.0], \forall i \in  [1, D]'
@@ -213,7 +213,7 @@ def __init__(self, ndim=None, bounds=None, f_shift="shift_data_4", f_matrix="M_4
     def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
-        z = np.dot(self.f_matrix, 5. * (x - self.f_shift) / 100) + 1
+        z = np.dot(self.f_matrix, 5. * (x - self.f_shift) / 100)
         return operator.expanded_griewank_rosenbrock_func(z) + self.f_bias
 
 
@@ -333,17 +333,16 @@ def __init__(self, ndim=None, bounds=None, f_shift="shift_data_7", f_matrix="M_7
         self.n3 = int(np.ceil(self.p[2] * self.ndim)) + self.n2
         self.idx1, self.idx2 = self.f_shuffle[:self.n1], self.f_shuffle[self.n1:self.n2]
         self.idx3, self.idx4 = self.f_shuffle[self.n2:self.n3], self.f_shuffle[self.n3:self.ndim]
-        self.g1 = operator.expanded_scaffer_f6_func
-        self.g2 = operator.hgbat_func
-        self.g3 = operator.rosenbrock_func
-        self.g4 = operator.modified_schwefel_func
         self.paras = {"f_shift": self.f_shift, "f_bias": self.f_bias, "f_matrix": self.f_matrix, "f_shuffle": self.f_shuffle}
 
     def evaluate(self, x, *args):
         self.n_fe += 1
         self.check_solution(x, self.dim_max, self.dim_supported)
         mz = np.dot(self.f_matrix, x - self.f_shift)
-        return self.g1(mz[self.idx1]) + self.g2(mz[self.idx2]) + self.g3(mz[self.idx3]) + self.g4(mz[self.idx4]) + self.f_bias
+        return (operator.expanded_schaffer_f6_func(mz[self.idx1]) +
+                operator.hgbat_func(mz[self.idx2], shift=-1.0) +
+                operator.rosenbrock_func(mz[self.idx3], shift=1.0) +
+                operator.modified_schwefel_func(mz[self.idx4]) + self.f_bias)
 
 
 class F72020(CecBenchmark):
@@ -399,11 +398,6 @@ def __init__(self, ndim=None, bounds=None, f_shift="shift_data_16", f_matrix="M_
         self.n4 = int(np.ceil(self.p[3] * self.ndim)) + self.n3
         self.idx1, self.idx2, self.idx3 = self.f_shuffle[:self.n1], self.f_shuffle[self.n1:self.n2], self.f_shuffle[self.n2:self.n3]
         self.idx4, self.idx5 = self.f_shuffle[self.n3:self.n4], self.f_shuffle[self.n4:self.ndim]
-        self.g1 = operator.expanded_scaffer_f6_func
-        self.g2 = operator.hgbat_func
-        self.g3 = operator.rosenbrock_func
-        self.g4 = operator.modified_schwefel_func
-        self.g5 = operator.elliptic_func
         self.paras = {"f_shift": self.f_shift, "f_bias": self.f_bias, "f_matrix": self.f_matrix, "f_shuffle": self.f_shuffle}
 
     def evaluate(self, x, *args):
@@ -412,8 +406,11 @@ def evaluate(self, x, *args):
         z = x - self.f_shift
         z1 = np.concatenate((z[self.idx1], z[self.idx2], z[self.idx3], z[self.idx4], z[self.idx5]))
         mz = np.dot(self.f_matrix, z1)
-        return self.g1(mz[:self.n1]) + self.g2(mz[self.n1:self.n2]) + self.g3(mz[self.n2:self.n3]) +\
-               self.g4(mz[self.n3:self.n4]) + self.g5(mz[self.n4:]) + self.f_bias
+        return (operator.expanded_scaffer_f6_func(mz[:self.n1]) +
+                operator.hgbat_func(mz[self.n1:self.n2], shift=-1.0) +
+                operator.rosenbrock_func(mz[self.n2:self.n3], shift=1.0) +
+                operator.modified_schwefel_func(mz[self.n3:self.n4]) +
+                operator.elliptic_func(mz[self.n4:]) + self.f_bias)
 
 
 class F82020(CecBenchmark):
diff --git a/opfunu/utils/operator.py b/opfunu/utils/operator.py
index 5863dea..5af7cc1 100644
--- a/opfunu/utils/operator.py
+++ b/opfunu/utils/operator.py
@@ -341,8 +341,8 @@ def expanded_schaffer_f6_func(x):
     f += 0.5 + (temp1_last - 0.5) / (temp2_last ** 2)
 
     return f
-	
-	
+
+
 def schaffer_f7_func(x):
     x = np.array(x).ravel()
     ndim = len(x)
@@ -353,58 +353,93 @@ def schaffer_f7_func(x):
     return (result / (ndim - 1)) ** 2
 
 
-def storn_chebyshev_polynomial_fitting_func(x, d=72.661):
+def chebyshev_func(x):
+    """
+    The following was converted from the cec2019 C code
+    Storn's Tchebychev - a 2nd ICEO function - generalized version
+    """
     x = np.array(x).ravel()
     ndim = len(x)
-    m = 32 * ndim
-    j1 = np.arange(0, ndim)
-    upper = ndim - j1
-
-    u = np.sum(x * 1.2 ** upper)
-    v = np.sum(x * (-1.2) ** upper)
-    p1 = 0 if u >= d else (u - d) ** 2
-    p2 = 0 if v >= d else (v - d) ** 2
-
-    wk = np.array([np.sum(x * (2. * k / m - 1) ** upper) for k in range(0, m + 1)])
-    conditions = [wk < 1, (1 <= wk) & (wk <= 1), wk > 1]
-    t1 = (wk + 1) ** 2
-    t2 = np.zeros(len(wk))
-    t3 = (wk - 1) ** 2
-    choices = [t1, t2, t3]
-    pk = np.select(conditions, choices, default=np.nan)
-    p3 = np.sum(pk)
-    return p1 + p2 + p3
-
-
-def inverse_hilbert_matrix_func(x):
+    sample = 32 * ndim
+
+    dx_arr = np.zeros(ndim)
+    dx_arr[:2] = [1.0, 1.2]
+    for i in range(2, ndim):
+        dx_arr[i] = 2.4 * dx_arr[i-1] - dx_arr[i-2]
+    dx = dx_arr[-1]
+
+    dy = 2.0 / sample
+
+    px, y, sum_val = 0, -1, 0
+    for i in range(sample + 1):
+        px = x[0]
+        for j in range(1, ndim):
+            px = y * px + x[j]
+        if px < -1 or px > 1:
+            sum_val += (1.0 - abs(px)) ** 2
+        y += dy
+
+    for _ in range(2):
+        px = np.sum(1.2 * x[1:]) + x[0]
+        mask = px < dx
+        sum_val += np.sum(px[mask] ** 2)
+
+    return sum_val
+
+
+def inverse_hilbert_func(x):
+    """
+    This is a direct conversion of the cec2019 C code for python optimized to use numpy
+    """
     x = np.array(x).ravel()
     ndim = len(x)
-    n = int(np.sqrt(ndim))
-    I = np.identity(n)
-    H = np.zeros((n, n))
-    Z = np.zeros((n, n))
-    for i in range(0, n):
-        for k in range(0, n):
-            Z[i, k] = x[i + n * k]
-            H[i, k] = 1. / (i + k + 1)
-    W = np.matmul(H, Z) - I
-    return np.sum(W)
-
-
-def lennard_jones_minimum_energy_cluster_func(x):
+    b = int(np.sqrt(ndim))
+
+    # Create the Hilbert matrix
+    i, j = np.indices((b, b))
+    hilbert = 1.0 / (i + j + 1)
+
+    # Reshape x and compute H*x
+    x = x.reshape((b, b))
+    y = np.dot(hilbert, x).dot(hilbert.T)
+
+    # Compute the absolute deviations
+    result = np.sum(np.abs(y - np.eye(b)))
+    return result
+
+
+def lennard_jones_func(x):
+    """
+    This version is a direct python conversion from the C-Code of CEC2019 implementation.
+    Find the atomic configuration with minimum energy (Lennard-Jones potential)
+    Valid for any dimension, D = 3 * k, k = 2, 3, 4, ..., 25.
+    k is the number of atoms in 3-D space.
+    """
     x = np.array(x).ravel()
     ndim = len(x)
-    result = 12.7120622568
-    n_upper = int(ndim / 3)
-    for i in range(1, n_upper):
-        for j in range(i + 1, n_upper + 1):
-            idx1, idx2 = 3 * i, 3 * j
-            dij = ((x[idx1 - 3] - x[idx2 - 3]) ** 2 + (x[idx1 - 2] - x[idx2 - 2]) ** 2 + (x[idx1 - 1] - x[idx2 - 1]) ** 2) ** 3
-            if dij == 0:
-                result += 0
+    # Minima values from Cambridge cluster database: http://www-wales.ch.cam.ac.uk/~jon/structures/LJ/tables.150.html
+    minima = np.array([-1., -3., -6., -9.103852, -12.712062, -16.505384, -19.821489, -24.113360,
+                       -28.422532, -32.765970, -37.967600, -44.326801, -47.845157, -52.322627, -56.815742,
+                       -61.317995, -66.530949, -72.659782, -77.1777043, -81.684571, -86.809782, -02.844472,
+                       -97.348815, -102.372663])
+
+    k = ndim // 3
+    sum_val = 0
+
+    x_matrix = x.reshape((k, 3))
+    for i in range(k-1):
+        for j in range(i + 1, k):
+            # Use slicing to get the differences between points i and j
+            diff = x_matrix[i] - x_matrix[j]
+            # Calculate the squared Euclidean distance
+            ed = np.sum(diff ** 2)
+            # Calculate ud and update sum_val accordingly
+            ud = ed ** 3
+            if ud > 1.0e-10:
+                sum_val += (1.0 / ud - 2.0) / ud
             else:
-                result += (1. / dij ** 2 - 2. / dij)
-    return result
+                sum_val += 1.0e20  # cec2019 version penalizes when ud is <=1e-10
+    return sum_val - minima[k - 2]  # Subtract known minima for k
 
 
 expanded_griewank_rosenbrock_func = grie_rosen_cec_func
diff --git a/tests/cec_based/test_cec2019.py b/tests/cec_based/test_cec2019.py
index 226df39..86f0ad7 100644
--- a/tests/cec_based/test_cec2019.py
+++ b/tests/cec_based/test_cec2019.py
@@ -146,3 +146,14 @@ def test_F102019_results():
     assert len(problem.lb) == ndim
     assert problem.bounds.shape[0] == ndim
     assert len(problem.x_global) == ndim
+
+
+def test_all_optimal_results():
+    known_failing = ['']
+    all_functions = [x for x in opfunu.get_all_cec_functions()
+                     if x.__name__[-4:] == '2019' and x.__name__ not in known_failing]
+    for function in all_functions:
+        problem = function(10)
+        x = problem.x_global
+        result = problem.evaluate(x)
+        assert abs(result - problem.f_global) <= problem.epsilon, f'{function.__name__} Failed Optimal Test'
diff --git a/tests/cec_based/test_cec2020.py b/tests/cec_based/test_cec2020.py
index aebc95a..0a87dc0 100644
--- a/tests/cec_based/test_cec2020.py
+++ b/tests/cec_based/test_cec2020.py
@@ -146,3 +146,15 @@ def test_F102020_results():
     assert len(problem.lb) == ndim
     assert problem.bounds.shape[0] == ndim
     assert len(problem.x_global) == ndim
+
+
+def test_all_optimal_results():
+    ndim = 30
+    known_failing = []
+    all_functions = [x for x in opfunu.get_all_cec_functions()
+                     if x.__name__[-4:] == '2020' and x.__name__ not in known_failing]
+    for function in all_functions:
+        problem = function(ndim=ndim)
+        x = problem.x_global
+        result = problem.evaluate(x)
+        assert abs(result - problem.f_global) <= problem.epsilon, f'{function.__name__} Failed Optimal Test'
diff --git a/tests/utils/__init__.py b/tests/utils/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/utils/test_operator.py b/tests/utils/test_operator.py
new file mode 100644
index 0000000..dc866c3
--- /dev/null
+++ b/tests/utils/test_operator.py
@@ -0,0 +1,11 @@
+import numpy as np
+import opfunu
+from opfunu.utils import operator
+
+
+def test_lennard_jones_func_zero_result():
+    """
+    The CEC2019 version when zero results in penalization of 1.0e20
+    """
+    x = np.zeros(18)
+    assert abs(operator.lennard_jones_func(x) - 1.5e21) <= 1e-8