diff --git a/src/combination/fubini.h b/src/combination/fubini.h
index 769f217..6d33218 100644
--- a/src/combination/fubini.h
+++ b/src/combination/fubini.h
@@ -83,45 +83,26 @@ IntegratorFubini<IntegratorFirst,IntegratorRest,N> integrator_fubini(IntegratorF
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
-    std::array<std::size_t,1> reorder{std::get<0>(dims)};
-    return integrator_reorder(integrator_fubini<1>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
+    std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}};
+    return integrator_reorder(integrator_fubini<1>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),{0,std::get<0>(dims)});
 } 
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
-    std::array<std::size_t,2> reorder{std::get<0>(dims),std::get<1>(dims)};
+    std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}, {1, std::get<1>(dims)}};
     return integrator_reorder(integrator_fubini<2>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
 } 
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
-    std::array<std::size_t,3> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims)};
+    std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}, {1, std::get<1>(dims)}, {2, std::get<2>(dims)}};
     return integrator_reorder(integrator_fubini<3>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
 } 
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
-    std::array<std::size_t,4> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims)};
+    std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}, {1, std::get<1>(dims)}, {2, std::get<2>(dims)}, {3, std::get<3>(dims)}};
     return integrator_reorder(integrator_fubini<4>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
 } 
 
-template<typename IntegratorFirst, typename IntegratorRest>
-auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
-    std::array<std::size_t,5> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims),std::get<4>(dims)};
-    return integrator_reorder(integrator_fubini<5>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
-} 
-
-template<typename IntegratorFirst, typename IntegratorRest>
-auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
-    std::array<std::size_t,6> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims),std::get<4>(dims),std::get<5>(dims)};
-    return integrator_reorder(integrator_fubini<6>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
-} 
-
-template<typename IntegratorFirst, typename IntegratorRest>
-auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
-    std::array<std::size_t,7> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims),std::get<4>(dims),std::get<5>(dims),std::get<6>(dims)};
-    return integrator_reorder(integrator_fubini<7>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
-} 
-
-
 };
\ No newline at end of file
diff --git a/src/combination/reorder.h b/src/combination/reorder.h
index 5df21c8..e5d21a3 100644
--- a/src/combination/reorder.h
+++ b/src/combination/reorder.h
@@ -11,37 +11,55 @@
 
 namespace viltrum {
 
-template<typename Integrator, std::size_t N>
+template<typename Integrator>
 class IntegratorReorder {
     Integrator integrator;
-    std::array<std::size_t,N> first_dimensions;
+    std::vector<std::tuple<std::size_t,std::size_t>> to_from;
+
+    template<typename V>
+    std::vector<typename V::value_type> reorder(const V& v, const typename V::value_type& def = typename V::value_type()) const {
+
+        std::vector<bool> tos(v.size(),false), froms(v.size(),false);
+        for (auto [to, from] : to_from) {
+            if (from<v.size()) { 
+                if (to>=tos.size()) tos.resize(to+1,false);
+                froms[from]=true; tos[to]=true;
+            } 
+        }
+        std::size_t ntos = 0, nfroms = 0;
+        std::for_each(tos.begin(),tos.end(),[&ntos] (bool b) { if (b) ++ntos; });
+        std::for_each(froms.begin(),froms.end(),[&nfroms] (bool b) { if (b) ++nfroms; });
+
+        std::vector<typename V::value_type> sol(std::max(tos.size(),ntos + v.size()-nfroms),def);
+        for (auto [to, from] : to_from) {
+            if ((from<v.size()) && (to<sol.size())) {
+                sol[to]=v[from];    
+            } 
+        }
+        std::size_t from, to = 0;
+        for (from = 0; from<v.size(); ++from) {
+            if (!froms[from]) { //We don't have this value yet
+                while ((to<sol.size()) && tos[to]) ++to; //Find where to put it
+                if (to<sol.size()) sol[to++] = v[from]; //Actually put it  
+            }   
+        } 
+        return sol;   
+    }
+
 public:
-    IntegratorReorder(const Integrator& i, const std::array<std::size_t,N>& fd) :
-        integrator(i), first_dimensions(fd) {}
+    IntegratorReorder(const Integrator& i, const std::vector<std::tuple<std::size_t,std::size_t>>& fd) :
+        integrator(i), to_from(fd) {}
 
 	template<typename Bins, std::size_t DIMBINS, typename F, typename Float, std::size_t DIM, typename Logger>
 	void integrate(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
 		const F& f, const Range<Float,DIM>& range, Logger& logger) const {
 
-        std::array<Float,DIM> reorder_min = range.min(); 
-        std::array<Float,DIM> reorder_max = range.max();
-        for (std::size_t i = 0; i<N; ++i) {
-            if (first_dimensions[i]<DIM) {
-                std::swap(reorder_min[i],reorder_min[first_dimensions[i]]);
-                std::swap(reorder_max[i],reorder_max[first_dimensions[i]]);
-            }
-        }
-        Range<Float,DIM> range_reordered(reorder_min,reorder_max);
-
-
+        auto range_reordered = viltrum::range<DIM>(reorder(range.min(),Float(0)),reorder(range.max(),Float(1)));
         const auto f_reordered = [&] (const std::array<Float,DIM>& x) {
-            std::array<Float,DIM> x_reordered = x;
-            for (std::size_t i = 0; i<N; ++i) {
-                if (first_dimensions[i]<DIM) {
-                    std::swap(x_reordered[i],x_reordered[first_dimensions[i]]);
-                }
-            }
-            return f(x_reordered);
+            auto re = reorder(x);
+            std::array<Float,DIM> x_re;
+            for (std::size_t i = 0; (i<DIM) && (i<re.size()); ++i) x_re[i] = re[i];  
+            return (f(x_re));
         };
 
         integrator.integrate(bins,bin_resolution, f_reordered, range_reordered, logger);
@@ -51,38 +69,31 @@ class IntegratorReorder {
 	void integrate(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
 		const F& f, const RangeInfinite<Float>& range, Logger& logger) const {
 
-        std::size_t size_reorder = 0; //For checking max size of the reorder vector
-        for (std::size_t i : first_dimensions) if (size_reorder<=i) size_reorder = (i+1);
-        
-        std::vector<Float> reorder_min = range.min(); 
-        if (reorder_min.size()<size_reorder) reorder_min.resize(size_reorder,0);
-        std::vector<Float> reorder_max = range.max();
-        if (reorder_max.size()<size_reorder) reorder_max.resize(size_reorder,1);
-
-        for (std::size_t i = 0; i<N; ++i) {
-            std::swap(reorder_min[i],reorder_min[first_dimensions[i]]);
-            std::swap(reorder_max[i],reorder_max[first_dimensions[i]]);
-        } 
-        RangeInfinite<Float> range_reordered(reorder_min,reorder_max);       
-        
+        RangeInfinite<Float> range_reordered(reorder(range.min(),Float(0)),reorder(range.max(),Float(1)));       
+         std::size_t size_reorder = 0;
+        for (auto [to, from] : to_from) {
+            if (size_reorder<=from) {
+                size_reorder = from + 1; 
+            } 
+        }       
         const auto f_reordered = [&] (const auto& x) {
-            std::vector<Float> x_reordered(size_reorder);
+            std::vector<Float> xr(size_reorder);
             auto it = x.begin();
-            for (Float& xr : x_reordered) { xr = (*it); ++it; }
-            
-            for (std::size_t i = 0; i<N; ++i) {
-                std::swap(x_reordered[i],x_reordered[first_dimensions[i]]);
-            }
-            return f(concat(x_reordered,sequence_from_iterators(it,x.end())));
+            for (Float& v : xr) { v = (*it); ++it; }
+/*           for (Float v : xr) std::cerr<<v<<" ";
+            std::cerr<<" -> ";
+            for (Float v : reorder(xr)) std::cerr<<v<<" ";
+            std::cerr<<std::endl; */
+            return f(concat(reorder(xr),sequence_from_iterators(it,x.end())));
         };
 
         integrator.integrate(bins,bin_resolution, f_reordered, range_reordered, logger);
     }
 };
 
-template<typename Integrator, std::size_t N>
-IntegratorReorder<Integrator,N> integrator_reorder(Integrator&& i, const std::array<std::size_t,N>& fd) {
-    return IntegratorReorder<std::decay_t<Integrator>,N>(std::forward<Integrator>(i),fd);
+template<typename Integrator>
+IntegratorReorder<Integrator> integrator_reorder(Integrator&& i, const std::vector<std::tuple<std::size_t,std::size_t>>& fd) {
+    return IntegratorReorder<std::decay_t<Integrator>>(std::forward<Integrator>(i),fd);
 }
 
 };
\ No newline at end of file
diff --git a/src/range.h b/src/range.h
index 1ddfd16..376a717 100644
--- a/src/range.h
+++ b/src/range.h
@@ -161,6 +161,18 @@ Range<T,N> range_all(const T& va, const T& vb) {
     return range(a,b);
 }
 
+template<std::size_t N, typename T>
+Range<T,N> range(const std::vector<T>& a, const std::vector<T>& b) {
+    std::array<T,N> ar, br;
+    std::size_t i;
+    for (i = 0; (i<N) && (i<a.size()); ++i) ar[i] = a[i];
+    for (;i<N; ++i) ar[i]=0;   
+    for (i = 0; (i<N) && (i<b.size()); ++i) br[i] = b[i];
+    for (;i<N; ++i) br[i]=1;
+    return range(ar,br); 
+}   
+
+
 template<std::size_t N, typename T = float>
 Range<T,N> range_primary() {
     return range_all<N>(T(0),T(1));