diff --git a/src/Makefile.am b/src/Makefile.am
index 5f650a54853edc8c2b350cc76c77eb6065e46337..377c4e50132928435738dba70e0205fb1aa0ece0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -5,6 +5,7 @@ bin_PROGRAMS = $(top_builddir)/sisso_cpp
 
 __top_builddir__sisso_cpp_SOURCES = \
     mpi_interface/MPI_Interface.cpp \
+    utils/math_funcs.cpp \
     feature_creation/node/value_storage/nodes_value_containers.cpp \
     feature_creation/units/Unit.cpp \
     feature_creation/node/Node.cpp \
diff --git a/src/descriptor_identifier/SISSORegressor.cpp b/src/descriptor_identifier/SISSORegressor.cpp
index 63dff34d862a482d024bcd38de27020ff951e24b..27b42e72048fbfeb42e3629989a52d78a0bf1b3a 100644
--- a/src/descriptor_identifier/SISSORegressor.cpp
+++ b/src/descriptor_identifier/SISSORegressor.cpp
@@ -137,32 +137,11 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
 
     std::vector<double> min_errors(_n_residual, util_funcs::norm(_prop.data(), _n_samp));
 
-    int n = _feat_space->phi_selected().size();
-
-    int length = n;
-    for(int nn = n - 1; nn > n - n_dim; --nn)
-        length *= nn;
-    length /= n_dim;
-
-    std::array<int, 2> start_end = _mpi_comm->get_start_end_for_iterator(length, 0);
-    std::vector<bool> v(n);
-    std::fill(v.end() - n_dim, v.end(), true);
-    int rr = 0;
-    do{
-        ++rr;
-    }while(std::next_permutation(v.begin(), v.end()) && (rr < start_end[0]));
+    for(int rr = 0; rr < n_dim; ++rr)
+        inds[rr] = _feat_space->phi_selected().size() - 1 - rr;
+    util_funcs::iterate(inds, inds.size(), _mpi_comm->rank());
 
     do {
-        int nn = 0;
-        int ii = 0;
-        do{
-            if (v[nn]) {
-                inds[ii] = nn;
-                ++ii;
-            }
-            ++nn;
-        }while(ii < n_dim);
-
         least_squares(inds, coefs.data());
         set_error(inds, coefs.data());
         double error = util_funcs::norm(_error.get(), _n_samp);
@@ -177,8 +156,7 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
             inds_min.erase(inds_min.end() - n_dim, inds_min.end());
             min_errors.pop_back();
         }
-        ++rr;
-    } while(std::next_permutation(v.begin(), v.end()) && (rr < start_end[1]));
+    } while(util_funcs::iterate(inds, inds.size(), _mpi_comm->size()));
 
     std::vector<double> all_min_error(_mpi_comm->size() * _n_residual);
     // std::vector<std::vector<std::vector<int>>> all_inds_min;
diff --git a/src/utils/math_funcs.hpp b/src/utils/math_funcs.hpp
index 32352c5813510731587942200533d40b87f59c1a..c4a5577fc0bb6946ace1224bfd4d8d67d14019da 100644
--- a/src/utils/math_funcs.hpp
+++ b/src/utils/math_funcs.hpp
@@ -86,6 +86,8 @@ namespace util_funcs
     {
         return std::abs(*std::max_element(start, start + size, [](T a, T b){return std::abs(a) < std::abs(b);}));
     }
+
+    bool iterate(std::vector<int>& inds, int size, int incriment);
 }
 
 #endif
\ No newline at end of file