diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp
index e36673c..d784ba5 100644
--- a/include/integratorxx/quadratures/radial/lmg.hpp
+++ b/include/integratorxx/quadratures/radial/lmg.hpp
@@ -1,5 +1,6 @@
 #include <integratorxx/util/lambert_w.hpp>
 #include <integratorxx/util/gamma.hpp>
+#include <integratorxx/util/pow.hpp>
 #include <integratorxx/quadratures/radial/radial_transform.hpp>
 
 namespace IntegratorXX {
@@ -7,11 +8,9 @@ namespace lmg {
 
 // Eq 19 of LMG paper (DOI 10.1007/s002140100263)
 inline double r_upper_obj(int m, double alpha, double r) {
-  const double am_term = (m + 1.0) / 2.0;
-  //const double g_term  = std::tgamma((m + 3.0) / 2.0);
   const double g_term = half_integer_tgamma<double>(m + 3);
   const double x = alpha * r * r;
-  return g_term * std::pow(x, am_term) * std::exp(-x);
+  return g_term * half_integer_pow(x, m+1) * std::exp(-x);
 }
 
 // Solve Eq 19 of LMG paper (DOI 10.1007/s002140100263)
@@ -21,7 +20,6 @@ inline double r_upper(int m, double alpha, double prec) {
   // X = -L * LAMBERT_W( - (P/G)^(1/L) / L )
   // R = SQRT(X / ALPHA)
   const double am_term = (m + 1.0) / 2.0;
-  //const double g_term  = std::tgamma((m + 3.0) / 2.0);
   const double g_term = half_integer_tgamma<double>(m + 3);
   const double arg = std::pow(prec/g_term, 1.0 / am_term) / am_term;
   const double wval = lambert_wm1(-arg); // W_(-1) is the larger value here
diff --git a/include/integratorxx/util/pow.hpp b/include/integratorxx/util/pow.hpp
new file mode 100644
index 0000000..c6ae500
--- /dev/null
+++ b/include/integratorxx/util/pow.hpp
@@ -0,0 +1,23 @@
+#pragma once
+#include <cmath>
+#include <integratorxx/util/type_traits.hpp>
+
+namespace IntegratorXX {
+
+template <typename FloatType, typename IntegralType>
+FloatType integer_pow(FloatType x, IntegralType n) {
+  if(n == 0) return 1;
+  if(n == 1) return x;
+  FloatType v = x;
+  for(IntegralType i = 1; i < n; ++i) v *= x;
+  return v;
+}
+
+template <typename FloatType, typename IntegralType>
+FloatType half_integer_pow(FloatType x, IntegralType n) {
+  assert(n >= 0);
+  if(n%2 == 0) return integer_pow(x, n/2);
+  else return std::sqrt(x) * integer_pow(x, n/2);
+}
+
+}
diff --git a/test/util.cxx b/test/util.cxx
index 7165247..b051fea 100644
--- a/test/util.cxx
+++ b/test/util.cxx
@@ -1,6 +1,7 @@
 #include "catch2/catch_all.hpp"
 #include <integratorxx/util/factorial.hpp>
 #include <integratorxx/util/gamma.hpp>
+#include <integratorxx/util/pow.hpp>
 #include <integratorxx/util/pow_2.hpp>
 #include <iostream>
 
@@ -34,6 +35,17 @@ TEST_CASE("Factorial", "[util]") {
   }
 }
 
+TEST_CASE("Pow", "[util]") {
+
+  using namespace IntegratorXX;
+  using namespace Catch::Matchers;
+  for(int i = 0; i < 10; ++i) {
+    REQUIRE_THAT(integer_pow(2.0, i), WithinAbs(std::pow(2.0,i), 1e-16));
+    REQUIRE_THAT(half_integer_pow(2.0, i), WithinAbs(std::pow(2.0,i/2.0), 1e-16));
+  }
+
+}
+
 TEST_CASE("Pow 2", "[util]") {
   using namespace IntegratorXX;