diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h
index 9cabcd9ab..24cfd7284 100644
--- a/include/convex_bodies/hpolytope.h
+++ b/include/convex_bodies/hpolytope.h
@@ -1010,9 +1010,10 @@ class HPolytope {
 
     // updates the velocity vector v and the position vector p after a reflection
     // the real value of p is given by p + moved_dist * v
+    // MT must be sparse, in RowMajor format
     template <typename update_parameters>
     auto compute_reflection(Point &v, Point &p, update_parameters const& params) const
-         -> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> { // MT must be in RowMajor format
+         -> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> {
             NT* v_data = v.pointerToData();
             NT* p_data = p.pointerToData();
             for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
@@ -1021,15 +1022,39 @@ class HPolytope {
             }
     }
 
-    template <typename update_parameters>
-    NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &vEv, update_parameters const &params) const {
+    // function to compute reflection for GaBW random walk
+    // compatible when the polytope is both dense or sparse
+    template <typename DenseSparseMT, typename update_parameters>
+    NT compute_reflection(Point &v, Point &p, NT &vEv, DenseSparseMT const &AE, VT const &AEA, update_parameters const &params) const {
+
+            NT new_vEv;
+            if constexpr (!std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
+                Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
+                VT x = v.getCoefficients();
+                new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
+                v += a;
+            } else {
+                
+                if constexpr(!std::is_same_v<DenseSparseMT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
+                    VT x = v.getCoefficients();
+                    new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
+                } else {
+                    new_vEv = vEv + 4.0 * params.inner_vi_ak * params.inner_vi_ak * AEA(params.facet_prev);
+                    NT* v_data = v.pointerToData();
+                    for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(AE, params.facet_prev); it; ++it) {
+                        new_vEv -= 4.0 * params.inner_vi_ak * it.value() * *(v_data + it.col());
+                    }
+                }
 
-            Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
-            VT x = v.getCoefficients();
-            NT new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
-            v += a;
+                NT* v_data = v.pointerToData();
+                NT* p_data = p.pointerToData();
+                for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
+                    *(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
+                    *(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
+                }
+            }
             NT coeff = std::sqrt(vEv / new_vEv);
-            v = v * coeff;
+            vEv = new_vEv;
             return coeff;
     }
 
diff --git a/include/random_walks/accelerated_billiard_walk_utils.hpp b/include/random_walks/accelerated_billiard_walk_utils.hpp
new file mode 100644
index 000000000..24cf32d76
--- /dev/null
+++ b/include/random_walks/accelerated_billiard_walk_utils.hpp
@@ -0,0 +1,133 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+// Copyright (c) 2024 Luca Perju
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef ACCELERATED_BILLIARD_WALK_UTILS_HPP
+#define ACCELERATED_BILLIARD_WALK_UTILS_HPP
+
+#include <Eigen/Eigen>
+#include <vector>
+
+const double eps = 1e-10;
+
+// data structure which maintains the values of (b - Ar)/Av, and can extract the minimum positive value and the facet associated with it
+// vec[i].first contains the value of (b(i) - Ar(i))/Av(i) + moved_dist, where moved_dist is the total distance that the point has travelled so far
+// The heap will only contain the values from vec which are greater than moved_dist (so that they are positive)
+template<typename NT>
+class BoundaryOracleHeap {
+public:
+    int n, heap_size;
+    std::vector<std::pair<NT, int>> heap;
+    std::vector<std::pair<NT, int>> vec;
+
+private:
+    int siftDown(int index) {
+        while((index << 1) + 1 < heap_size) {
+            int child = (index << 1) + 1;
+            if(child + 1 < heap_size && heap[child + 1].first < heap[child].first - eps) {
+                child += 1;
+            }
+            if(heap[child].first < heap[index].first - eps)
+            {
+                std::swap(heap[child], heap[index]);
+                std::swap(vec[heap[child].second].second, vec[heap[index].second].second);
+                index = child;
+            } else {
+                return index;
+            }
+        }
+        return index;
+    }
+
+    int siftUp(int index) {
+        while(index > 0 && heap[(index - 1) >> 1].first - eps > heap[index].first) {
+            std::swap(heap[(index - 1) >> 1], heap[index]);
+            std::swap(vec[heap[(index - 1) >> 1].second].second, vec[heap[index].second].second);
+            index = (index - 1) >> 1;
+        }
+        return index;
+    }
+
+    // takes the index of a facet, and (in case it is in the heap) removes it from the heap.
+    void remove (int index) {
+        index = vec[index].second;
+        if(index == -1) {
+            return;
+        }
+        std::swap(heap[heap_size - 1], heap[index]);
+        std::swap(vec[heap[heap_size - 1].second].second, vec[heap[index].second].second);
+        vec[heap[heap_size - 1].second].second = -1;
+        heap_size -= 1;
+        index = siftDown(index);
+        siftUp(index);
+    }
+
+    // inserts a new value into the heap, with its associated facet
+    void insert (const std::pair<NT, int> val) {
+        vec[val.second].second = heap_size;
+        vec[val.second].first = val.first;
+        heap[heap_size++] = val;
+        siftUp(heap_size - 1);
+    }
+
+public:
+    BoundaryOracleHeap() {}
+
+    BoundaryOracleHeap(int n) : n(n), heap_size(0) {
+        heap.resize(n);
+        vec.resize(n);
+    }
+
+    // rebuilds the heap with the existing values from vec
+    // O(n)
+    void rebuild (const NT &moved_dist) {
+        heap_size = 0;
+        for(int i = 0; i < n; ++i) {
+            vec[i].second = -1;
+            if(vec[i].first - eps > moved_dist) {
+                vec[i].second = heap_size;
+                heap[heap_size++] = {vec[i].first, i};
+            }
+        }
+        for(int i = heap_size - 1; i >= 0; --i) {
+            siftDown(i);
+        }
+    }
+
+    // returns (b(i) - Ar(i))/Av(i) + moved_dist
+    // O(1)
+    NT get_val (const int &index) {
+        return vec[index].first;
+    }
+
+    // returns the nearest facet
+    // O(1)
+    std::pair<NT, int> get_min () {
+        return heap[0];
+    }
+
+    // changes the stored value for a given facet, and updates the heap accordingly
+    // O(logn)
+    void change_val(const int& index, const NT& new_val, const NT& moved_dist) {
+        if(new_val < moved_dist - eps) {
+            vec[index].first = new_val;
+            remove(index);
+        } else {
+            if(vec[index].second == -1) {
+                insert({new_val, index});
+            } else {
+                heap[vec[index].second].first = new_val;
+                vec[index].first = new_val;
+                siftDown(vec[index].second);
+                siftUp(vec[index].second);
+            }
+        }
+    }
+};
+
+
+#endif
diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp
index 532d7672e..a8a7d334e 100644
--- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp
+++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp
@@ -21,6 +21,7 @@
 #include "generators/boost_random_number_generator.hpp"
 #include "sampling/ellipsoid.hpp"
 #include "random_walks/uniform_billiard_walk.hpp"
+#include "random_walks/accelerated_billiard_walk_utils.hpp"
 
 #include "random_walks/compute_diameter.hpp"
 
@@ -47,12 +48,13 @@ struct GaussianAcceleratedBilliardWalk
     struct update_parameters
     {
         update_parameters()
-                :   facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0)
+                :   facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0), moved_dist(0.0)
         {}
         int facet_prev;
         bool hit_ball;
         double inner_vi_ak;
         double ball_inner_norm;
+        double moved_dist;
     };
 
     parameters param;
@@ -67,9 +69,14 @@ struct GaussianAcceleratedBilliardWalk
     struct Walk
     {
         typedef typename Polytope::PointType Point;
-        typedef typename Polytope::DenseMT DenseMT;
+        typedef typename Polytope::MT MT;
+        typedef typename Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> DenseMT;
         typedef typename Polytope::VT VT;
         typedef typename Point::FT NT;
+        // AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise
+        using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >; 
+        // AE is sparse rowMajor if (MT is sparse rowMajor and E is sparse), and Dense otherwise
+        using AE_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>> && std::is_base_of<Eigen::SparseMatrixBase<E_type>, E_type>::value, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>, DenseMT >; 
 
         void computeLcov(const E_type E)
         {
@@ -110,7 +117,11 @@ struct GaussianAcceleratedBilliardWalk
             _L = compute_diameter<GenericPolytope>::template compute<NT>(P);
             computeLcov(E);
             _E = E;
-            _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
+            if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::value) {
+                _AA = (P.get_mat() * P.get_mat().transpose());
+            } else {
+                _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
+            }
             _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
             initialize(P, p, rng);
         }
@@ -131,7 +142,11 @@ struct GaussianAcceleratedBilliardWalk
                                 ::template compute<NT>(P);
             computeLcov(E);
             _E = E;
-            _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
+            if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::value) {
+                _AA = (P.get_mat() * P.get_mat().transpose());
+            } else {
+                _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
+            }
             _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
             initialize(P, p, rng);
         }
@@ -162,13 +177,7 @@ struct GaussianAcceleratedBilliardWalk
                 Point p0 = _p;
 
                 it = 0;
-                std::pair<NT, int> pbpair;
-                if(!was_reset) {
-                    pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters);
-                } else {
-                    pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);
-                    was_reset = false;
-                }                                                    
+                std::pair<NT, int> pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);
 
                 if (T <= pbpair.first) {
                     _p += (T * _v);
@@ -177,18 +186,40 @@ struct GaussianAcceleratedBilliardWalk
                 }
 
                 _lambda_prev = dl * pbpair.first;
-                _p += (_lambda_prev * _v);
+                if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
+                    typename Point::Coeff b;
+                    NT* b_data;
+                    b = P.get_vec();
+                    b_data = b.data();
+                    _update_parameters.moved_dist = _lambda_prev;
+                    NT* Ar_data = _lambdas.data();
+                    NT* Av_data = _Av.data();
+                    for(int i = 0; i < P.num_of_hyperplanes(); ++i) {
+                        _distances_set.vec[i].first = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i));
+                    }
+                    // rebuild the heap with the new values of (b - Ar) / Av
+                    _distances_set.rebuild(_update_parameters.moved_dist);
+                } else {
+                    _p += (_lambda_prev * _v);
+                }
                 T -= _lambda_prev;
-                coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
+                
+                T = T * coef;
+                coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
+                T = T / coef;
+
                 it++;
 
                 while (it < _rho)
                 {
-                    std::pair<NT, int> pbpair
-                        = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters);
-                    _Av *= coef;
-                    _update_parameters.inner_vi_ak *= coef;
-                    pbpair.first /= coef;
+                    std::pair<NT, int> pbpair;
+                    if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
+                        pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev,
+                                                           _distances_set, _AA, _update_parameters);
+                    } else {
+                        pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev,
+                                                           _AA, _update_parameters);
+                    }
 
                     if (T <= pbpair.first) {
                         _p += (T * _v);
@@ -196,14 +227,23 @@ struct GaussianAcceleratedBilliardWalk
                         break;
                     }
                     _lambda_prev = dl * pbpair.first;
-                    _p += (_lambda_prev * _v);
+                    if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
+                        _update_parameters.moved_dist += _lambda_prev;
+                    } else {
+                        _p += (_lambda_prev * _v);
+                    }
                     T -= _lambda_prev;
-                    coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
+                    
+                    T = T * coef;
+                    coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
+                    T = T / coef;
+
                     it++;
                 }
+                _p += _update_parameters.moved_dist * _v;
+                _update_parameters.moved_dist = 0.0;
                 if (it == _rho){
                     _p = p0;
-                    was_reset = true;
                 }
             }
 
@@ -219,21 +259,24 @@ struct GaussianAcceleratedBilliardWalk
         {
             unsigned int n = P.dimension();
             const NT dl = 0.995;
-            was_reset = false;
             _lambdas.setZero(P.num_of_hyperplanes());
             _Av.setZero(P.num_of_hyperplanes());
             _p = p;
-            _AE.noalias() = (DenseMT)(P.get_mat() * _E);
-            _AEA = _AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum();
-            /*
-            _AEA.resize(P.num_of_hyperplanes());
-            for(int i = 0; i < P.num_of_hyperplanes(); ++i)
             {
-                _AEA(i) = _AE.row(i).dot(P.get_mat().row(i));
-            }*/
+                DenseMT temp_AE;
+                temp_AE.noalias() = (DenseMT)(P.get_mat() * _E);
+                _AEA = temp_AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum();
+                if constexpr (std::is_same_v<AE_type, DenseMT>)
+                    _AE = temp_AE;
+                else
+                    _AE = temp_AE.sparseView();
+            }
+            _distances_set = BoundaryOracleHeap<NT>(P.num_of_hyperplanes());
+
 
             _v = GetDirection<Point>::apply(n, rng, false);
             _v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());
+            
 
             NT T = -std::log(rng.sample_urdist()) * _L;
             Point p0 = _p;
@@ -251,15 +294,15 @@ struct GaussianAcceleratedBilliardWalk
             _lambda_prev = dl * pbpair.first;
             _p += (_lambda_prev * _v);
             T -= _lambda_prev;
-            coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
+            
+            T = T * coef;
+            coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
+            T = T / coef;
 
             while (it <= _rho)
             {
                 std::pair<NT, int> pbpair
                         = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters);
-                _Av *= coef;
-                _update_parameters.inner_vi_ak *= coef;
-                pbpair.first /= coef;
 
                 if (T <= pbpair.first) {
                     _p += (T * _v);
@@ -273,7 +316,11 @@ struct GaussianAcceleratedBilliardWalk
                 _lambda_prev = dl * pbpair.first;
                 _p += (_lambda_prev * _v);
                 T -= _lambda_prev;
-                coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
+
+                T = T * coef;
+                coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
+                T = T / coef;
+
                 it++;
             }
         }
@@ -282,16 +329,16 @@ struct GaussianAcceleratedBilliardWalk
         Point _p;
         Point _v;
         NT _lambda_prev;
-        DenseMT _AA;
+        AA_type _AA;
         E_type _L_cov;   // LL' = inv(E)
-        DenseMT _AE;
+        AE_type _AE;
         E_type _E;
         VT _AEA;
         unsigned int _rho;
         update_parameters _update_parameters;
         typename Point::Coeff _lambdas;
         typename Point::Coeff _Av;
-        bool was_reset;
+        BoundaryOracleHeap<NT> _distances_set;
     };
 
 };
diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp
index f91abee4b..f68a21a2c 100644
--- a/include/random_walks/uniform_accelerated_billiard_walk.hpp
+++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp
@@ -13,125 +13,7 @@
 
 #include "sampling/sphere.hpp"
 #include <Eigen/Eigen>
-#include <set>
-#include <vector>
-
-const double eps = 1e-10;
-
-// data structure which maintains the values of (b - Ar)/Av, and can extract the minimum positive value and the facet associated with it
-// vec[i].first contains the value of (b(i) - Ar(i))/Av(i) + moved_dist, where moved_dist is the total distance that the point has travelled so far
-// The heap will only contain the values from vec which are greater than moved_dist (so that they are positive)
-template<typename NT>
-class BoundaryOracleHeap {
-public:
-    int n, heap_size;
-    std::vector<std::pair<NT, int>> heap;
-    std::vector<std::pair<NT, int>> vec;
-
-private:
-    int siftDown(int index) {
-        while((index << 1) + 1 < heap_size) {
-            int child = (index << 1) + 1;
-            if(child + 1 < heap_size && heap[child + 1].first < heap[child].first - eps) {
-                child += 1;
-            }
-            if(heap[child].first < heap[index].first - eps)
-            {
-                std::swap(heap[child], heap[index]);
-                std::swap(vec[heap[child].second].second, vec[heap[index].second].second);
-                index = child;
-            } else {
-                return index;
-            }
-        }
-        return index;
-    }
-
-    int siftUp(int index) {
-        while(index > 0 && heap[(index - 1) >> 1].first - eps > heap[index].first) {
-            std::swap(heap[(index - 1) >> 1], heap[index]);
-            std::swap(vec[heap[(index - 1) >> 1].second].second, vec[heap[index].second].second);
-            index = (index - 1) >> 1;
-        }
-        return index;
-    }
-
-    // takes the index of a facet, and (in case it is in the heap) removes it from the heap.
-    void remove (int index) {
-        index = vec[index].second;
-        if(index == -1) {
-            return;
-        }
-        std::swap(heap[heap_size - 1], heap[index]);
-        std::swap(vec[heap[heap_size - 1].second].second, vec[heap[index].second].second);
-        vec[heap[heap_size - 1].second].second = -1;
-        heap_size -= 1;
-        index = siftDown(index);
-        siftUp(index);
-    }
-
-    // inserts a new value into the heap, with its associated facet
-    void insert (const std::pair<NT, int> val) {
-        vec[val.second].second = heap_size;
-        vec[val.second].first = val.first;
-        heap[heap_size++] = val;
-        siftUp(heap_size - 1);
-    }
-
-public:
-    BoundaryOracleHeap() {}
-
-    BoundaryOracleHeap(int n) : n(n), heap_size(0) {
-        heap.resize(n);
-        vec.resize(n);
-    }
-
-    // rebuilds the heap with the existing values from vec
-    // O(n)
-    void rebuild (const NT &moved_dist) {
-        heap_size = 0;
-        for(int i = 0; i < n; ++i) {
-            vec[i].second = -1;
-            if(vec[i].first - eps > moved_dist) {
-                vec[i].second = heap_size;
-                heap[heap_size++] = {vec[i].first, i};
-            }
-        }
-        for(int i = heap_size - 1; i >= 0; --i) {
-            siftDown(i);
-        }
-    }
-
-    // returns (b(i) - Ar(i))/Av(i) + moved_dist
-    // O(1)
-    NT get_val (const int &index) {
-        return vec[index].first;
-    }
-
-    // returns the nearest facet
-    // O(1)
-    std::pair<NT, int> get_min () {
-        return heap[0];
-    }
-
-    // changes the stored value for a given facet, and updates the heap accordingly
-    // O(logn)
-    void change_val(const int& index, const NT& new_val, const NT& moved_dist) {
-        if(new_val < moved_dist - eps) {
-            vec[index].first = new_val;
-            remove(index);
-        } else {
-            if(vec[index].second == -1) {
-                insert({new_val, index});
-            } else {
-                heap[vec[index].second].first = new_val;
-                vec[index].first = new_val;
-                siftDown(vec[index].second);
-                siftUp(vec[index].second);
-            }
-        }
-    }
-};
+#include "random_walks/accelerated_billiard_walk_utils.hpp"
 
 
 // Billiard walk which accelarates each step for uniform distribution
@@ -275,7 +157,7 @@ struct AcceleratedBilliardWalk
                 it++;
 
                 while (it < _rho)
-                {   
+                {
                     std::pair<NT, int> pbpair;
                     if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
                         pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev,
diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp
index a15f40102..ba71d1398 100644
--- a/test/sampling_test.cpp
+++ b/test/sampling_test.cpp
@@ -320,8 +320,11 @@ void call_test_gabw(){
     Point StartingPoint(d);
     std::list<Point> randPoints;
 
+    std::chrono::time_point<std::chrono::high_resolution_clock> start, stop;
+    start = std::chrono::high_resolution_clock::now();
+
     std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Skinny-H-cube10" << std::endl;
-    P = generate_skinny_cube<Hpolytope>(10);
+    P = generate_skinny_cube<Hpolytope>(d);
 
 
     Point p = P.ComputeInnerBall().first;
@@ -340,6 +343,11 @@ void call_test_gabw(){
     RandomPointGenerator::apply(P, p, E, numpoints, 1, randPoints,
                                 push_back_policy, rng);
 
+    stop = std::chrono::high_resolution_clock::now();
+
+    std::chrono::duration<double> total_time = stop - start;
+    std::cout << "Done in " << total_time.count() << '\n';
+
     MT samples(d, numpoints);
     unsigned int jj = 0;
     for (typename std::list<Point>::iterator rpit = randPoints.begin(); rpit != randPoints.end(); rpit++, jj++)
@@ -352,11 +360,11 @@ void call_test_gabw(){
     RNGType Srng(d);
 
     typedef Eigen::SparseMatrix<NT> SparseMT;
-    typedef HPolytope<Point, SparseMT> SparseHpoly;
+    typedef HPolytope<Point, Eigen::SparseMatrix<NT, Eigen::RowMajor>> SparseHpoly;
     std::list<Point> Points;
 
     SparseHpoly SP;
-    SP = generate_skinny_cube<SparseHpoly>(10);
+    SP = generate_skinny_cube<SparseHpoly>(d);
 
 
     std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Sparse Skinny-H-cube10" << std::endl;
@@ -371,15 +379,21 @@ void call_test_gabw(){
             > sparsewalk;
     typedef MultivariateGaussianRandomPointGenerator <sparsewalk> SparseRandomPointGenerator;
 
+    start = std::chrono::high_resolution_clock::now();
 
     ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_ELLIPSOID>
-    (SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0));
+    ((SparseMT)SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0));
 
     const SparseMT SE = get<0>(ellipsoid).sparseView();
 
     SparseRandomPointGenerator::apply(SP, p, SE, numpoints, 1, Points,
                                 push_back_policy, Srng);
 
+    stop = std::chrono::high_resolution_clock::now();
+
+    total_time = stop - start;
+    std::cout << "Done in " << total_time.count() << '\n';
+
     jj = 0;
     MT Ssamples(d, numpoints);
     for (typename std::list<Point>::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++)