diff --git a/mex/sources/sobol/gaussian.hh b/mex/sources/sobol/gaussian.hh
index a10a80777427fe73a9c72bbcf75c88d184a6049d..39dfc759a92711c38a0ab29551f40236fb067fd2 100644
--- a/mex/sources/sobol/gaussian.hh
+++ b/mex/sources/sobol/gaussian.hh
@@ -28,6 +28,7 @@
 #include <algorithm>
 #include <array>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <numbers>
 #include <vector>
@@ -41,7 +42,7 @@ using namespace std;
 constexpr double lb = .02425;
 constexpr double ub = .97575;
 
-template<typename T>
+template<floating_point T>
 T
 icdf(const T uniform)
 /*
@@ -105,7 +106,7 @@ icdf(const T uniform)
 }
 
 void
-icdfm(int n, auto* U)
+icdfm(int n, floating_point auto* U)
 {
 #pragma omp parallel for
   for (int i = 0; i < n; i++)
@@ -114,7 +115,7 @@ icdfm(int n, auto* U)
 }
 
 void
-icdfmSigma(int d, int n, auto* U, const double* LowerCholSigma)
+icdfmSigma(int d, int n, floating_point auto* U, const double* LowerCholSigma)
 {
   double one = 1.0;
   double zero = 0.0;
@@ -127,7 +128,7 @@ icdfmSigma(int d, int n, auto* U, const double* LowerCholSigma)
 }
 
 void
-usphere(int d, int n, auto* U)
+usphere(int d, int n, floating_point auto* U)
 {
   icdfm(n * d, U);
 #pragma omp parallel for
@@ -145,7 +146,7 @@ usphere(int d, int n, auto* U)
 }
 
 void
-usphereRadius(int d, int n, double radius, auto* U)
+usphereRadius(int d, int n, double radius, floating_point auto* U)
 {
   icdfm(n * d, U);
 #pragma omp parallel for
diff --git a/mex/sources/sobol/sobol.hh b/mex/sources/sobol/sobol.hh
index 026ccfa261df279efc8ef555b8861753d195d4cd..1cbaee6c5bd34127a9f900397998a685dcedfd5c 100644
--- a/mex/sources/sobol/sobol.hh
+++ b/mex/sources/sobol/sobol.hh
@@ -1,6 +1,6 @@
 /* Interface to quasi Monte Carlo sequences (à la Sobol) routines.
  *
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -21,6 +21,7 @@
 #ifndef SOBOL_HH
 #define SOBOL_HH
 
+#include <concepts>
 #include <cstdint>
 #include <dynblas.h> // For the FORTRAN_WRAPPER macro
 
@@ -48,7 +49,7 @@ sobol_block(int dimension, int block_size, int64_t seed, double* block)
   return seed;
 }
 
-template<typename T>
+template<floating_point T>
 void
 expand_unit_hypercube(int dimension, int block_size, T* block, const T* lower_bound,
                       const T* upper_bound)