diff --git a/src/bgmCompare_logp_and_grad.cpp b/src/bgmCompare_logp_and_grad.cpp
index 45adf60c..32794912 100644
--- a/src/bgmCompare_logp_and_grad.cpp
+++ b/src/bgmCompare_logp_and_grad.cpp
@@ -2,6 +2,7 @@
 #include "bgmCompare_helper.h"
 #include "bgmCompare_logp_and_grad.h"
 #include <cmath>
+#include "explog_switch.h"
 
 using namespace Rcpp;
 
@@ -174,7 +175,7 @@ double log_pseudoposterior(
 
   // ---- priors ----
   auto log_beta_prior = [&](double x) {
-    return x * main_alpha - std::log1p(std::exp(x)) * (main_alpha + main_beta);
+    return x * main_alpha - std::log1p(MY_EXP(x)) * (main_alpha + main_beta);
   };
 
   // Main effects prior
@@ -510,7 +511,7 @@ arma::vec gradient(
       for (int c = 0; c < num_cats; ++c) {
         off = main_index(base + c, 0);
         double value = main_effects(base + c, 0);
-        const double p = 1.0 / (1.0 + std::exp(-value));
+        const double p = 1.0 / (1.0 + MY_EXP(-value));
         grad(off) += main_alpha - (main_alpha + main_beta) * p;
 
         if (inclusion_indicator(v, v) == 0) continue;
@@ -524,12 +525,12 @@ arma::vec gradient(
     } else {
       off = main_index(base, 0);
       double value = main_effects(base, 0);
-      double p = 1.0 / (1.0 + std::exp(-value));
+      double p = 1.0 / (1.0 + MY_EXP(-value));
       grad(off) += main_alpha - (main_alpha + main_beta) * p;
 
       off = main_index(base + 1, 0);
       value = main_effects(base + 1, 0);
-      p = 1.0 / (1.0 + std::exp(-value));
+      p = 1.0 / (1.0 + MY_EXP(-value));
       grad(off) += main_alpha - (main_alpha + main_beta) * p;
 
 
@@ -733,7 +734,7 @@ double log_pseudoposterior_main_component(
   if (h == 0) {
     // overall
     auto log_beta_prior = [&](double x) {
-      return x * main_alpha - std::log1p(std::exp(x)) * (main_alpha + main_beta);
+      return x * main_alpha - std::log1p(MY_EXP(x)) * (main_alpha + main_beta);
     };
 
     // Main effects prior
diff --git a/src/bgmCompare_sampler.cpp b/src/bgmCompare_sampler.cpp
index 264898ac..1d74b2d0 100644
--- a/src/bgmCompare_sampler.cpp
+++ b/src/bgmCompare_sampler.cpp
@@ -12,6 +12,7 @@
 #include "print_mutex.h"
 #include "rng_utils.h"
 #include "sampler_output.h"
+#include "explog_switch.h"
 #include <string>
 
 using namespace Rcpp;
@@ -125,7 +126,7 @@ void impute_missing_bgmcompare(
       for(int category = 1; category <= num_categories(variable, group); category++) {
         exponent = group_main_effects(category - 1);
         exponent += category * rest_score;
-        cumsum += std::exp(exponent);
+        cumsum += MY_EXP(exponent);
         category_response_probabilities[category] = cumsum;
       }
     } else {
@@ -137,7 +138,7 @@ void impute_missing_bgmcompare(
           (category - baseline_category[variable]) *
           (category - baseline_category[variable]);
         exponent += category * rest_score;
-        cumsum += std::exp(exponent);
+        cumsum += MY_EXP(exponent);
         category_response_probabilities[category] = cumsum;
       }
     }
@@ -985,7 +986,7 @@ void tune_proposal_sd_bgmcompare(
           SamplerResult result = rwm_sampler(current, prop_sd, log_post, rng);
           current = result.state[0];
           prop_sd = update_proposal_sd_with_robbins_monro(
-            prop_sd, std::log(result.accept_prob), rm_weight, target_accept
+            prop_sd, MY_LOG(result.accept_prob), rm_weight, target_accept
           );
         }
       }
@@ -1013,7 +1014,7 @@ void tune_proposal_sd_bgmcompare(
           SamplerResult result = rwm_sampler(current, prop_sd, log_post, rng);
           current = result.state[0];
           prop_sd = update_proposal_sd_with_robbins_monro(
-            prop_sd, std::log(result.accept_prob), rm_weight, target_accept
+            prop_sd, MY_LOG(result.accept_prob), rm_weight, target_accept
           );
         }
       }
@@ -1046,7 +1047,7 @@ void tune_proposal_sd_bgmcompare(
         SamplerResult result = rwm_sampler(current, prop_sd, log_post, rng);
         current = result.state[0];
         prop_sd = update_proposal_sd_with_robbins_monro(
-          prop_sd, std::log(result.accept_prob), rm_weight, target_accept
+          prop_sd, MY_LOG(result.accept_prob), rm_weight, target_accept
         );
       }
     }
@@ -1165,11 +1166,11 @@ void update_indicator_differences_metropolis_bgmcompare (
     // Add prior inclusion probability contribution
     double inc_prob = inclusion_probability_difference(var, var);
     if(proposed_ind == 1) {
-      log_accept += std::log(inc_prob);
-      log_accept -= std::log(1 - inc_prob);
+      log_accept += MY_LOG(inc_prob);
+      log_accept -= MY_LOG(1 - inc_prob);
     } else {
-      log_accept -= std::log(inc_prob);
-      log_accept += std::log(1 - inc_prob);
+      log_accept -= MY_LOG(inc_prob);
+      log_accept += MY_LOG(1 - inc_prob);
     }
 
     // Add parameter prior contribution
@@ -1201,7 +1202,7 @@ void update_indicator_differences_metropolis_bgmcompare (
 
     // Perform Metropolis-Hastings step
     double U = runif(rng);
-    if(std::log(U) < log_accept) {
+    if(MY_LOG(U) < log_accept) {
       inclusion_indicator(var, var) = proposed_ind;
       main_effects.rows(start, stop).cols(1, num_groups - 1) =
         proposed_main_effects.rows(start, stop).cols(1, num_groups - 1);
@@ -1243,11 +1244,11 @@ void update_indicator_differences_metropolis_bgmcompare (
     // Add prior inclusion probability contribution
     double inc_prob = inclusion_probability_difference(var1, var2);
     if(proposed_ind == 1) {
-      log_accept += std::log(inc_prob);
-      log_accept -= std::log(1 - inc_prob);
+      log_accept += MY_LOG(inc_prob);
+      log_accept -= MY_LOG(1 - inc_prob);
     } else {
-      log_accept -= std::log(inc_prob);
-      log_accept += std::log(1 - inc_prob);
+      log_accept -= MY_LOG(inc_prob);
+      log_accept += MY_LOG(1 - inc_prob);
     }
 
     // Add parameter prior contribution
@@ -1280,7 +1281,7 @@ void update_indicator_differences_metropolis_bgmcompare (
 
     // Metropolis-Hastings acceptance step
     double U = runif(rng);
-    if (std::log(U) < log_accept) {
+    if (MY_LOG(U) < log_accept) {
       // Update inclusion inclusion_indicator
       inclusion_indicator(var1, var2) = proposed_ind;
       inclusion_indicator(var2, var1) = proposed_ind;
diff --git a/src/bgm_logp_and_grad.cpp b/src/bgm_logp_and_grad.cpp
index fb308927..a5e06fd8 100644
--- a/src/bgm_logp_and_grad.cpp
+++ b/src/bgm_logp_and_grad.cpp
@@ -2,6 +2,7 @@
 #include "bgm_helper.h"
 #include "bgm_logp_and_grad.h"
 #include "common_helpers.h"
+#include "explog_switch.h"
 
 using namespace Rcpp;
 
@@ -58,7 +59,7 @@ double log_pseudoposterior_main_effects_component (
   double log_posterior = 0.0;
 
   auto log_beta_prior = [&](double main_effect_param) {
-    return main_effect_param * main_alpha - std::log1p (std::exp (main_effect_param)) * (main_alpha + main_beta);
+    return main_effect_param * main_alpha - std::log1p (MY_EXP (main_effect_param)) * (main_alpha + main_beta);
   };
 
   const int num_cats = num_categories(variable);
@@ -278,7 +279,7 @@ double log_pseudoposterior (
 
   // Calculate the contribution from the data and the prior
   auto log_beta_prior = [&](double main_effect_param) {
-    return main_effect_param * main_alpha - std::log1p (std::exp (main_effect_param)) * (main_alpha + main_beta);
+    return main_effect_param * main_alpha - std::log1p (MY_EXP (main_effect_param)) * (main_alpha + main_beta);
   };
 
   for (int variable = 0; variable < num_variables; variable++) {
@@ -531,14 +532,14 @@ arma::vec gradient_log_pseudoposterior (
     if (is_ordinal_variable(variable)) {
       const int num_cats = num_categories(variable);
       for (int cat = 0; cat < num_cats; cat++) {
-        const double p = 1.0 / (1.0 + std::exp (-main_effects(variable, cat)));
+        const double p = 1.0 / (1.0 + MY_EXP (-main_effects(variable, cat)));
         gradient(offset + cat) += main_alpha - (main_alpha + main_beta) * p;
       }
       offset += num_cats;
     } else {
       for (int i = 0; i < 2; i++) {
         const double main_effect_param = main_effects(variable, i);
-        const double p = 1.0 / (1.0 + std::exp (-main_effect_param));
+        const double p = 1.0 / (1.0 + MY_EXP (-main_effect_param));
         gradient(offset + i) += main_alpha - (main_alpha + main_beta) * p;
       }
       offset += 2;
@@ -721,8 +722,8 @@ double compute_log_likelihood_ratio_for_variable (
       for (int person = 0; person < num_persons; person++) {
         const double base = main + score * residual_scores[person] - bounds[person];
 
-        const double exp_current = std::exp(base + score * interaction[person] * current_state);
-        const double exp_proposed = std::exp(base + score * interaction[person] * proposed_state);
+        const double exp_current = MY_EXP(base + score * interaction[person] * current_state);
+        const double exp_proposed = MY_EXP(base + score * interaction[person] * proposed_state);
 
         denom_current[person] += exp_current;
         denom_proposed[person] += exp_proposed;
diff --git a/src/bgm_sampler.cpp b/src/bgm_sampler.cpp
index 562348ee..ede4a489 100644
--- a/src/bgm_sampler.cpp
+++ b/src/bgm_sampler.cpp
@@ -95,14 +95,14 @@ void impute_missing_bgm (
       for (int cat = 0; cat < num_cats; cat++) {
         const int score = cat + 1;
         const double exponent = main_effects (variable, cat) + score * residual_score;
-        cumsum += std::exp (exponent);
+        cumsum += MY_EXP (exponent);
         category_probabilities[score] = cumsum;
       }
     } else {
       // Compute probabilities for Blume-Capel variable
       const int ref = baseline_category (variable);
 
-      cumsum = std::exp (main_effects (variable, 1) * ref * ref);
+      cumsum = MY_EXP (main_effects (variable, 1) * ref * ref);
       category_probabilities[0] = cumsum;
 
       for (int cat = 0; cat < num_cats; cat++) {
@@ -112,7 +112,7 @@ void impute_missing_bgm (
           main_effects (variable, 0) * score +
           main_effects (variable, 1) * centered * centered +
           score * residual_score;
-        cumsum += std::exp (exponent);
+        cumsum += MY_EXP (exponent);
         category_probabilities[score] = cumsum;
       }
     }
@@ -792,7 +792,7 @@ void tune_proposal_sd_bgm(
       }
 
       proposal_sd = update_proposal_sd_with_robbins_monro(
-        proposal_sd, std::log(result.accept_prob), rm_weight, target_accept
+        proposal_sd, MY_LOG(result.accept_prob), rm_weight, target_accept
       );
 
       proposal_sd_pairwise_effects(variable1, variable2) = proposal_sd;
@@ -891,15 +891,15 @@ void update_indicator_edges_metropolis_bgm (
     if (proposing_addition) {
       log_accept += R::dcauchy(proposed_state, 0.0, pairwise_scale, true);
       log_accept -= R::dnorm(proposed_state, current_state, sd, true);
-      log_accept += std::log (inclusion_probability_ij) - std::log (1.0 - inclusion_probability_ij);
+      log_accept += MY_LOG (inclusion_probability_ij) - MY_LOG (1.0 - inclusion_probability_ij);
     } else {
       log_accept -= R::dcauchy(current_state, 0.0, pairwise_scale, true);
       log_accept += R::dnorm(current_state, proposed_state, sd, true);
-      log_accept -= std::log (inclusion_probability_ij) - std::log (1.0 - inclusion_probability_ij);
+      log_accept -= MY_LOG (inclusion_probability_ij) - MY_LOG (1.0 - inclusion_probability_ij);
     }
 
     // Metropolis-Hastings accept step
-    if (std::log (runif(rng)) < log_accept) {
+    if (MY_LOG (runif(rng)) < log_accept) {
       const int updated_indicator = 1 - indicator(variable1, variable2);
       indicator(variable1, variable2) = updated_indicator;
       indicator(variable2, variable1) = updated_indicator;
diff --git a/src/mcmc_adaptation.h b/src/mcmc_adaptation.h
index d8257445..f414a93d 100644
--- a/src/mcmc_adaptation.h
+++ b/src/mcmc_adaptation.h
@@ -7,7 +7,7 @@
 #include <vector>
 #include "mcmc_utils.h"
 #include "mcmc_rwm.h"
-
+#include "explog_switch.h"
 
 class DualAveraging {
 public:
@@ -21,10 +21,10 @@ class DualAveraging {
   int t;
 
   DualAveraging(double initial_step_size)
-    : log_step_size(std::log(initial_step_size)),
-      log_step_size_avg(std::log(initial_step_size)),
+    : log_step_size(MY_LOG(initial_step_size)),
+      log_step_size_avg(MY_LOG(initial_step_size)),
       hbar(0.0),
-      mu(std::log(10.0 * initial_step_size)),
+      mu(MY_LOG(10.0 * initial_step_size)),
       gamma(0.05),
       t0(10.0),
       kappa(0.75),
@@ -42,15 +42,15 @@ class DualAveraging {
   }
 
   void restart(double new_step_size) {
-    log_step_size = std::log(new_step_size);
-    log_step_size_avg = std::log(new_step_size);
-    mu = std::log(10.0 * new_step_size);
+    log_step_size = MY_LOG(new_step_size);
+    log_step_size_avg = MY_LOG(new_step_size);
+    mu = MY_LOG(10.0 * new_step_size);
     hbar = 0.0;
     t = 1;
   }
 
-  double current() const { return std::exp(log_step_size); }
-  double averaged() const { return std::exp(log_step_size_avg); }
+  double current() const { return MY_EXP(log_step_size); }
+  double averaged() const { return MY_EXP(log_step_size_avg); }
 };
 
 
@@ -281,7 +281,7 @@ class RWMAdaptationController {
     for (arma::uword i = 0; i < proposal_sd.n_rows; ++i) {
       for (arma::uword j = 0; j < proposal_sd.n_cols; ++j) {
         if (index_mask(i, j) == 1) {
-          const double accept_log = std::log(accept_prob_matrix(i, j));
+          const double accept_log = MY_LOG(accept_prob_matrix(i, j));
           proposal_sd(i, j) = update_proposal_sd_with_robbins_monro(
             proposal_sd(i, j),
             accept_log,
diff --git a/src/mcmc_hmc.cpp b/src/mcmc_hmc.cpp
index 1db70c5d..8617a2e2 100644
--- a/src/mcmc_hmc.cpp
+++ b/src/mcmc_hmc.cpp
@@ -30,9 +30,9 @@ SamplerResult hmc_sampler(
   double proposed_H = -log_post(theta) + kinetic_energy(r, inv_mass_diag);
   double log_accept_prob = current_H - proposed_H;
 
-  arma::vec state = (std::log(runif(rng)) < log_accept_prob) ? theta : init_theta;
+  arma::vec state = (MY_LOG(runif(rng)) < log_accept_prob) ? theta : init_theta;
 
-  double accept_prob = std::min(1.0, std::exp(log_accept_prob));
+  double accept_prob = std::min(1.0, MY_EXP(log_accept_prob));
 
   return {state, accept_prob};
 }
diff --git a/src/mcmc_nuts.cpp b/src/mcmc_nuts.cpp
index 3e49937f..81f74d20 100644
--- a/src/mcmc_nuts.cpp
+++ b/src/mcmc_nuts.cpp
@@ -92,7 +92,7 @@ BuildTreeResult build_tree(
     int n_new = 1 * (log_u <= logp - kin);
     int s_new = 1 * (log_u <= Delta_max + logp - kin);
     bool divergent = (s_new == 0);
-    double alpha = std::min(1.0, std::exp(logp - kin - logp0 + kin0));
+    double alpha = std::min(1.0, MY_EXP(logp - kin - logp0 + kin0));
 
     return {
       theta_new, r_new, theta_new, r_new, theta_new, n_new, s_new, alpha, 1, divergent
diff --git a/src/mcmc_rwm.cpp b/src/mcmc_rwm.cpp
index 3c847f1a..bce6d827 100644
--- a/src/mcmc_rwm.cpp
+++ b/src/mcmc_rwm.cpp
@@ -34,7 +34,7 @@ SamplerResult rwm_sampler(
 ) {
   double proposed_state = rnorm(rng, current_state, step_size);
   double log_accept = log_post(proposed_state) - log_post(current_state);
-  double accept_prob = std::min(1.0, std::exp(log_accept));
+  double accept_prob = std::min(1.0, MY_EXP(log_accept));
 
   double state = (runif(rng) < accept_prob) ? proposed_state : current_state;
 
diff --git a/src/mcmc_utils.cpp b/src/mcmc_utils.cpp
index 9fe32df1..3941c178 100644
--- a/src/mcmc_utils.cpp
+++ b/src/mcmc_utils.cpp
@@ -74,10 +74,10 @@ double heuristic_initial_step_size(
   double H0 = logp0 - kin0;
   double H1 = logp1 - kin1;
 
-  int direction = 2 * (H1 - H0 > std::log(0.5)) - 1;  // +1 or -1
+  int direction = 2 * (H1 - H0 > MY_LOG(0.5)) - 1;  // +1 or -1
 
   int attempts = 0;
-  while (direction * (H1 - H0) > -direction * std::log(2.0) && attempts < max_attempts) {
+  while (direction * (H1 - H0) > -direction * MY_LOG(2.0) && attempts < max_attempts) {
     eps = (direction == 1) ? 2.0 * eps : 0.5 * eps;
 
     // One leapfrog step
diff --git a/src/mcmc_utils.h b/src/mcmc_utils.h
index d6116cbd..9c935706 100644
--- a/src/mcmc_utils.h
+++ b/src/mcmc_utils.h
@@ -5,6 +5,7 @@
 #include <cmath>
 #include <functional>
 #include <memory>
+#include "explog_switch.h"
 struct SafeRNG;
 
 // (only if <algorithm> didn’t already provide it under C++17)
@@ -114,7 +115,7 @@ inline void update_step_size_with_dual_averaging (
     arma::vec& state,
     const double target_acceptance
 ) {
-  const double target_log_step_size = std::log (10.0 * initial_step_size);
+  const double target_log_step_size = MY_LOG (10.0 * initial_step_size);
   constexpr int stabilization_offset = 10;
 
   double& log_step_size = state[0];
@@ -160,9 +161,9 @@ inline void update_step_size_with_robbins_monro (
   const double error = acceptance_probability - target_acceptance;
   const double decay = std::pow (static_cast<double> (iteration), -decay_rate);
 
-  double log_step_size = std::log (step_size);
+  double log_step_size = MY_LOG (step_size);
   log_step_size += error * decay;
-  step_size = std::exp (log_step_size);
+  step_size = MY_EXP (log_step_size);
 }
 
 
@@ -191,7 +192,7 @@ inline double update_proposal_sd_with_robbins_monro (
   // Normalize the acceptance probability
   double observed_acceptance_probability = 1.0;
   if (observed_log_acceptance_probability < 0.0) {
-    observed_acceptance_probability = std::exp (observed_log_acceptance_probability);
+    observed_acceptance_probability = MY_EXP (observed_log_acceptance_probability);
   }
 
   // Robbins-Monro update step
diff --git a/src/sbm_edge_prior.cpp b/src/sbm_edge_prior.cpp
index 17d77823..a1f9dd22 100644
--- a/src/sbm_edge_prior.cpp
+++ b/src/sbm_edge_prior.cpp
@@ -1,6 +1,7 @@
 #include <RcppArmadillo.h>
 #include <Rcpp.h>
 #include "rng_utils.h"
+#include "explog_switch.h"
 
 using namespace Rcpp;
 
@@ -71,7 +72,7 @@ arma::vec compute_Vn_mfm_sbm(arma::uword no_variables,
       arma::vec b_linspace_2 = arma::linspace((k+1)*dirichlet_alpha,(k+1)*dirichlet_alpha+no_variables-1, no_variables); // denominator b*e*(b*e+1)*...*(b*e+p-1)
       double b = arma::accu(arma::log(b_linspace_1))-arma::accu(arma::log(b_linspace_2)) + R::dpois((k+1)-1, lambda, true); // sum(log(numerator)) - sum(log(denominator)) + log(P=(k+1|lambda))
       double m = std::max(b,r);  // scaling factor for log-sum-exp formula
-      r = std::log(std::exp(r-m) +  std::exp(b-m)) + m; // update r using log-sum-exp formula to ensure numerical stability and avoid underflow
+      r = MY_LOG(MY_EXP(r-m) +  MY_EXP(b-m)) + m; // update r using log-sum-exp formula to ensure numerical stability and avoid underflow
     }
     log_Vn(t) = r;
   }
@@ -92,14 +93,14 @@ double log_likelihood_mfm_sbm(arma::uvec cluster_assign,
     if(j != node) {
       if(j < node) {
         output += indicator(j, node) *
-          std::log(cluster_probs(cluster_assign(j), cluster_assign(node)));
+          MY_LOG(cluster_probs(cluster_assign(j), cluster_assign(node)));
         output += (1 - indicator(j, node)) *
-          std::log(1 - cluster_probs(cluster_assign(j), cluster_assign(node)));
+          MY_LOG(1 - cluster_probs(cluster_assign(j), cluster_assign(node)));
       } else {
         output += indicator(node, j) *
-          std::log(cluster_probs(cluster_assign(node), cluster_assign(j)));
+          MY_LOG(cluster_probs(cluster_assign(node), cluster_assign(j)));
         output += (1 - indicator(node, j)) *
-          std::log(1 - cluster_probs(cluster_assign(node), cluster_assign(j)));
+          MY_LOG(1 - cluster_probs(cluster_assign(node), cluster_assign(j)));
       }
     }
   }
@@ -221,7 +222,7 @@ arma::uvec block_allocations_mfm_sbm(arma::uvec cluster_assign,
                                            no_variables);
 
           prob = (static_cast<double>(dirichlet_alpha) + static_cast<double>(cluster_size_node(c))) *
-            std::exp(loglike);
+            MY_EXP(loglike);
           }
           else{ // if old group, the probability is set to 0.0
             prob = 0.0;
@@ -236,8 +237,8 @@ arma::uvec block_allocations_mfm_sbm(arma::uvec cluster_assign,
                                          beta_bernoulli_beta);
 
           prob = static_cast<double>(dirichlet_alpha) *
-            std::exp(logmarg) *
-            std::exp(log_Vn(no_clusters - 1) - log_Vn(no_clusters - 2));
+            MY_EXP(logmarg) *
+            MY_EXP(log_Vn(no_clusters - 1) - log_Vn(no_clusters - 2));
         }
 
         cluster_prob(c) = prob;
@@ -278,7 +279,7 @@ arma::uvec block_allocations_mfm_sbm(arma::uvec cluster_assign,
                                            no_variables);
 
           prob = (static_cast<double>(dirichlet_alpha) + static_cast<double>(cluster_size_node(c))) *
-            std::exp(loglike);
+            MY_EXP(loglike);
         } else {
           logmarg = log_marginal_mfm_sbm(cluster_assign_tmp,
                                          indicator,
@@ -288,8 +289,8 @@ arma::uvec block_allocations_mfm_sbm(arma::uvec cluster_assign,
                                          beta_bernoulli_beta);
 
           prob = static_cast<double>(dirichlet_alpha) *
-            std::exp(logmarg) *
-            std::exp(log_Vn(no_clusters) - log_Vn(no_clusters-1));
+            MY_EXP(logmarg) *
+            MY_EXP(log_Vn(no_clusters) - log_Vn(no_clusters-1));
         }
 
         cluster_prob(c) = prob;