diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index d0350adbfadccb6dd5e32ae55115f59fc48f7a41..ce40c2c1de3a680f0a811f07e5d9ad6d44d511f2 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -35,16 +35,19 @@
 
 #define hydro_gamma 1.66666666666666667f
 #define hydro_gamma_minus_one 0.66666666666666667f
+#define hydro_one_over_gamma_minus_one 1.5f
 
 #elif defined(HYDRO_GAMMA_4_3)
 
 #define hydro_gamma 1.33333333333333333f
 #define hydro_gamma_minus_one 0.33333333333333333f
+#define hydro_one_over_gamma_minus_one 3.f
 
 #elif defined(HYDRO_GAMMA_2_1)
 
 #define hydro_gamma 2.f
 #define hydro_gamma_minus_one 1.f
+#define hydro_one_over_gamma_minus_one 1.f
 
 #endif
 
@@ -72,8 +75,39 @@ __attribute__((always_inline)) INLINE static float pow_gamma(float x) {
   return x * x;
 
 #else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Returns the argument to the power given by the adiabatic index minus
+ * one
+ *
+ * Computes $x^(\gamma - 1)$.
+ */
+__attribute__((always_inline)) INLINE static float pow_gamma_minus_one(
+    float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return cbrtf(x * x);
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return cbrtf(x);
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return x;
+
+#else
+
   error("The adiabatic index is not defined !");
   return 0.f;
+
 #endif
 }
 
@@ -99,8 +133,10 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
   return 1.f / x;
 
 #else
+
   error("The adiabatic index is not defined !");
   return 0.f;
+
 #endif
 }
 
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index d365c595f2bff9aad3754435377ea548af68313b..753025b2c7d454be8613ef69b6555c786219f378 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -914,7 +914,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
   Pi_ij.v *= (wi_dr.v + wj_dr.v);
 
   /* Thermal conductivity */
-  v_sig_u.v = vec_sqrt(vec_set1(2.f * (const_hydro_gamma - 1.f)) *
+  v_sig_u.v = vec_sqrt(vec_set1(2.f * hydro_gamma_minus_one) *
                        vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) /
                        (pirho.v + pjrho.v));
   tc.v = vec_set1(const_conductivity_alpha) * v_sig_u.v / (pirho.v + pjrho.v);
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index d0ca5817dda1243957ad1b28c13cee18ab9aea21..7789d4a80db6b37003fcbff92e5f8b42b1ca5ff3 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -247,7 +247,5 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
 
   const float entropy = p->entropy + p->entropy_dt * dt;
 
-  return entropy * powf(p->rho, hydro_gamma - 1.f) *
-         (1.f / (hydro_gamma - 1.f));
-  return 1;
+  return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
 }
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index b5cceacc31bfaeec0941c855a1c193a9bd36c4a6..848de3b5de8b916c35a3c6da7e414fd9a35d966b 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -17,8 +17,8 @@
  *
  ******************************************************************************/
 
-#include "approx_math.h"
 #include "adiabatic_index.h"
+#include "approx_math.h"
 
 /**
  * @brief Computes the hydro time-step of a given particle