diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root
new file mode 120000
index 0000000..945c9b4
--- /dev/null
+++ b/_codeql_detected_source_root
@@ -0,0 +1 @@
+.
\ No newline at end of file
diff --git a/include/astray/core/geodesic.hpp b/include/astray/core/geodesic.hpp
index e750519..d5f73b2 100644
--- a/include/astray/core/geodesic.hpp
+++ b/include/astray/core/geodesic.hpp
@@ -40,11 +40,15 @@ class geodesic
         {
           value_type dydt;
           dydt.head(4) = y.tail(4);
-          auto christoffel_symbols = metric.christoffel_symbols(y.head(4));
+          const auto christoffel_symbols = metric.christoffel_symbols(y.head(4));
+          const auto& velocity = y.tail(4);
           for (auto i = 0; i < 4; ++i)
+          {
+            const auto vel_i = velocity[i];
             for (auto j = 0; j < 4; ++j)
               for (auto k = 0; k < 4; ++k)
-                dydt.tail(4)[k] -= christoffel_symbols(i, j, k) * y.tail(4)[i] * y.tail(4)[j];
+                dydt.tail(4)[k] -= christoffel_symbols(i, j, k) * vel_i * velocity[j];
+          }
           return dydt;
         }
       }, 
diff --git a/include/astray/core/ray_tracer.hpp b/include/astray/core/ray_tracer.hpp
index 18198d3..f98011a 100644
--- a/include/astray/core/ray_tracer.hpp
+++ b/include/astray/core/ray_tracer.hpp
@@ -207,7 +207,7 @@ class ray_tracer
     resized_data_type_  = mpi::data_type(subarray_data_type_, 0    , partitioner_.block_size()[0] * sizeof(pixel_type));
 #endif
     
-    std::visit([value] (auto& visited) 
+    std::visit([&value] (auto& visited) 
     { 
       visited.aspect_ratio = static_cast<scalar_type>(value[0]) / static_cast<scalar_type>(value[1]); 
     }, observer_.get_projection());
diff --git a/include/astray/math/coordinate_system.hpp b/include/astray/math/coordinate_system.hpp
index 611a575..d573b33 100644
--- a/include/astray/math/coordinate_system.hpp
+++ b/include/astray/math/coordinate_system.hpp
@@ -4,6 +4,7 @@
 
 #include <astray/math/angle.hpp>
 #include <astray/math/constants.hpp>
+#include <astray/math/ipow.hpp>
 #include <astray/math/linear_algebra.hpp>
 
 namespace ast
@@ -47,7 +48,9 @@ __device__ __host__ constexpr void convert    (type& value)
   {
     if      constexpr (target == coordinate_system_type::cylindrical)
     {
-      const scalar rho = std::sqrt (static_cast<scalar>(std::pow(value[1], 2)) + static_cast<scalar>(std::pow(value[2], 2)));
+      const scalar x_sq = ipow<2>(value[1]);
+      const scalar y_sq = ipow<2>(value[2]);
+      const scalar rho = std::sqrt (x_sq + y_sq);
       const scalar phi = std::atan2(value[2], value[1]);
 
       value[1] = rho;
@@ -57,10 +60,10 @@ __device__ __host__ constexpr void convert    (type& value)
     }
     else if constexpr (target == coordinate_system_type::spherical)
     {
-      const scalar r     = std::sqrt (
-        static_cast<scalar>(std::pow(value[1], 2)) + 
-        static_cast<scalar>(std::pow(value[2], 2)) + 
-        static_cast<scalar>(std::pow(value[3], 2)));
+      const scalar x_sq = ipow<2>(value[1]);
+      const scalar y_sq = ipow<2>(value[2]);
+      const scalar z_sq = ipow<2>(value[3]);
+      const scalar r     = std::sqrt (x_sq + y_sq + z_sq);
       const scalar theta = std::acos (value[3] / r);
       const scalar phi   = std::atan2(value[2], value[1]);
 
@@ -87,11 +90,11 @@ __device__ __host__ constexpr void convert    (type& value)
     {
       wrap_angles<source>(value);
 
-      const scalar r     = std::sqrt(
-        static_cast<scalar>(std::pow(value[1], 2)) +
-        static_cast<scalar>(std::pow(value[3], 2)));
-      const scalar theta = std::atan2(value[1], value[3]);
-      const scalar phi   = value[2];
+      const scalar rho_sq = ipow<2>(value[1]);
+      const scalar z_sq   = ipow<2>(value[3]);
+      const scalar r      = std::sqrt(rho_sq + z_sq);
+      const scalar theta  = std::atan2(value[1], value[3]);
+      const scalar phi    = value[2];
 
       value[1] = r;
       value[2] = theta;
@@ -139,9 +142,9 @@ __device__ __host__ constexpr void convert    (type& value, const scalar free_pa
     {
       wrap_angles<source>(value);
 
-      const scalar sq_r2_p_a2 = std::sqrt(
-        static_cast<scalar>(std::pow(value[1]      , 2)) + 
-        static_cast<scalar>(std::pow(free_parameter, 2)));
+      const scalar r_sq   = ipow<2>(value[1]);
+      const scalar a_sq   = ipow<2>(free_parameter);
+      const scalar sq_r2_p_a2 = std::sqrt(r_sq + a_sq);
 
       const scalar x = sq_r2_p_a2 * std::sin(value[2]) * std::cos(value[3]);
       const scalar y = sq_r2_p_a2 * std::sin(value[2]) * std::sin(value[3]);
@@ -166,13 +169,14 @@ __device__ __host__ constexpr void convert    (type& value, const scalar free_pa
   {
     if      constexpr (target == coordinate_system_type::boyer_lindquist)
     {
-      const scalar w = 
-        static_cast<scalar>(std::pow(value[1]      , 2)) + 
-        static_cast<scalar>(std::pow(value[2]      , 2)) + 
-        static_cast<scalar>(std::pow(value[3]      , 2)) - 
-        static_cast<scalar>(std::pow(free_parameter, 2));
+      const scalar x_sq = ipow<2>(value[1]);
+      const scalar y_sq = ipow<2>(value[2]);
+      const scalar z_sq = ipow<2>(value[3]);
+      const scalar a_sq = ipow<2>(free_parameter);
+      const scalar w = x_sq + y_sq + z_sq - a_sq;
 
-      const scalar r     = std::sqrt (static_cast<scalar>(0.5) * (w + std::sqrt(static_cast<scalar>(std::pow(w, 2)) + static_cast<scalar>(4) * static_cast<scalar>(std::pow(free_parameter, 2)) * static_cast<scalar>(std::pow(value[3], 2)))));
+      const scalar w_sq = ipow<2>(w);
+      const scalar r    = std::sqrt (static_cast<scalar>(0.5) * (w + std::sqrt(w_sq + static_cast<scalar>(4) * a_sq * z_sq)));
       const scalar theta = std::acos (value[3] / r);
       const scalar phi   = std::atan2(value[2], value[1]);
 
@@ -184,10 +188,12 @@ __device__ __host__ constexpr void convert    (type& value, const scalar free_pa
     }
     else if constexpr (target == coordinate_system_type::prolate_spheroidal)
     {
-      const scalar x_squared        = static_cast<scalar>(std::pow(value[1], 2));
-      const scalar y_squared        = static_cast<scalar>(std::pow(value[2], 2));
-      const scalar first_component  = std::sqrt(x_squared + y_squared + static_cast<scalar>(std::pow(value[3] + free_parameter, 2)));
-      const scalar second_component = std::sqrt(x_squared + y_squared + static_cast<scalar>(std::pow(value[3] - free_parameter, 2)));
+      const scalar x_squared        = ipow<2>(value[1]);
+      const scalar y_squared        = ipow<2>(value[2]);
+      const scalar z_plus_a         = value[3] + free_parameter;
+      const scalar z_minus_a        = value[3] - free_parameter;
+      const scalar first_component  = std::sqrt(x_squared + y_squared + ipow<2>(z_plus_a));
+      const scalar second_component = std::sqrt(x_squared + y_squared + ipow<2>(z_minus_a));
       const scalar _2a              = 2 * free_parameter;
 
       const scalar sigma = (first_component + second_component) / _2a;
@@ -220,7 +226,9 @@ __device__ __host__ constexpr void convert    (type& value, const scalar free_pa
     {
       wrap_angles<source>(value);
 
-      const scalar common = free_parameter * std::sqrt((static_cast<scalar>(std::pow(value[1], 2)) - static_cast<scalar>(1)) * (static_cast<scalar>(1) - static_cast<scalar>(std::pow(value[2], 2))));
+      const scalar s_sq   = ipow<2>(value[1]);
+      const scalar t_sq   = ipow<2>(value[2]);
+      const scalar common = free_parameter * std::sqrt((s_sq - static_cast<scalar>(1)) * (static_cast<scalar>(1) - t_sq));
       const scalar x      = common * std::cos(value[3]);
       const scalar y      = common * std::sin(value[3]);
       const scalar z      = free_parameter * value[1] * value[2];
@@ -376,7 +384,9 @@ __device__ __host__ constexpr void convert_ray(type& value, const scalar free_pa
       const scalar& r = value.position[1];
       const scalar& t = value.position[2];
       const scalar& p = value.position[3];
-      const scalar  k = std::sqrt(static_cast<scalar>(std::pow(free_parameter, 2)) + static_cast<scalar>(std::pow(r, 2)));
+      const scalar  a_sq = ipow<2>(free_parameter);
+      const scalar  r_sq = ipow<2>(r);
+      const scalar  k = std::sqrt(a_sq + r_sq);
 
       matrix44<scalar> transform;
       transform << 
@@ -409,7 +419,9 @@ __device__ __host__ constexpr void convert_ray(type& value, const scalar free_pa
       const scalar& r = value.position[1];
       const scalar& t = value.position[2];
       const scalar& p = value.position[3];
-      const scalar  k = std::sqrt(static_cast<scalar>(std::pow(free_parameter, 2)) + static_cast<scalar>(std::pow(r, 2)));
+      const scalar  a_sq = ipow<2>(free_parameter);
+      const scalar  r_sq = ipow<2>(r);
+      const scalar  k = std::sqrt(a_sq + r_sq);
 
       matrix44<scalar> transform;
       transform << 
@@ -429,13 +441,15 @@ __device__ __host__ constexpr void convert_ray(type& value, const scalar free_pa
       const scalar& t = value.position[2];
       const scalar& p = value.position[3];
       const scalar& a = free_parameter;
-      const scalar  k = std::sqrt(-(static_cast<scalar>(std::pow(s, 2)) - static_cast<scalar>(1)) * (static_cast<scalar>(std::pow(t, 2)) - static_cast<scalar>(1)));
+      const scalar  s_sq = ipow<2>(s);
+      const scalar  t_sq = ipow<2>(t);
+      const scalar  k = std::sqrt(-(s_sq - static_cast<scalar>(1)) * (t_sq - static_cast<scalar>(1)));
 
       matrix44<scalar> transform;
       transform << 
         1,                                                                    0,                                                                    0,                    0,
-        0, -a * s * (static_cast<scalar>(std::pow(t, 2)) - 1) * std::cos(p) / k, -a * t * (static_cast<scalar>(std::pow(s, 2)) - 1) * std::cos(p) / k, -a * std::sin(p) * k,
-        0, -a * s * (static_cast<scalar>(std::pow(t, 2)) - 1) * std::sin(p) / k, -a * t * (static_cast<scalar>(std::pow(s, 2)) - 1) * std::sin(p) / k,  a * std::cos(p) * k,
+        0, -a * s * (t_sq - 1) * std::cos(p) / k, -a * t * (s_sq - 1) * std::cos(p) / k, -a * std::sin(p) * k,
+        0, -a * s * (t_sq - 1) * std::sin(p) / k, -a * t * (s_sq - 1) * std::sin(p) / k,  a * std::cos(p) * k,
         0,  a * t                                                              ,  a * s                                                              ,                    0;
       transform = transform.inverse().eval(); // TODO: Ideally inversion should be analytic without runtime overhead.
 
@@ -463,13 +477,15 @@ __device__ __host__ constexpr void convert_ray(type& value, const scalar free_pa
       const scalar& t = value.position[2];
       const scalar& p = value.position[3];
       const scalar& a = free_parameter;
-      const scalar  k = std::sqrt(-(static_cast<scalar>(std::pow(s, 2)) - 1) * (static_cast<scalar>(std::pow(t, 2)) - 1));
+      const scalar  s_sq = ipow<2>(s);
+      const scalar  t_sq = ipow<2>(t);
+      const scalar  k = std::sqrt(-(s_sq - static_cast<scalar>(1)) * (t_sq - static_cast<scalar>(1)));
 
       matrix44<scalar> transform;
       transform << 
         1,                                                                    0,                                                                    0,                    0,
-        0, -a * s * (static_cast<scalar>(std::pow(t, 2)) - 1) * std::cos(p) / k, -a * t * (static_cast<scalar>(std::pow(s, 2)) - 1) * std::cos(p) / k, -a * std::sin(p) * k,
-        0, -a * s * (static_cast<scalar>(std::pow(t, 2)) - 1) * std::sin(p) / k, -a * t * (static_cast<scalar>(std::pow(s, 2)) - 1) * std::sin(p) / k,  a * std::cos(p) * k,
+        0, -a * s * (t_sq - 1) * std::cos(p) / k, -a * t * (s_sq - 1) * std::cos(p) / k, -a * std::sin(p) * k,
+        0, -a * s * (t_sq - 1) * std::sin(p) / k, -a * t * (s_sq - 1) * std::sin(p) / k,  a * std::cos(p) * k,
         0,  a * t                                                              ,  a * s                                                              ,                    0;
 
       convert<source, target>(value.position, free_parameter);
diff --git a/include/astray/math/ipow.hpp b/include/astray/math/ipow.hpp
new file mode 100644
index 0000000..699951a
--- /dev/null
+++ b/include/astray/math/ipow.hpp
@@ -0,0 +1,36 @@
+#pragma once
+
+namespace ast
+{
+// Compile-time integer power function for improved code clarity and intent.
+// Computes base^N where N is a compile-time constant integer.
+// More efficient than std::pow for small integer powers while maintaining intent.
+template <int N, typename T>
+__device__ __host__ constexpr T ipow(const T& base)
+{
+  static_assert(N >= 0, "ipow only supports non-negative integer powers");
+  
+  if constexpr (N == 0)
+    return T(1);
+  else if constexpr (N == 1)
+    return base;
+  else if constexpr (N == 2)
+    return base * base;
+  else if constexpr (N == 3)
+    return base * base * base;
+  else if constexpr (N == 4)
+  {
+    const T base_sq = base * base;
+    return base_sq * base_sq;
+  }
+  else if constexpr (N % 2 == 0)
+  {
+    const T half_power = ipow<N / 2>(base);
+    return half_power * half_power;
+  }
+  else
+  {
+    return base * ipow<N - 1>(base);
+  }
+}
+}
diff --git a/include/astray/math/ode/error/controller/integral_controller.hpp b/include/astray/math/ode/error/controller/integral_controller.hpp
index 2416311..28ab485 100644
--- a/include/astray/math/ode/error/controller/integral_controller.hpp
+++ b/include/astray/math/ode/error/controller/integral_controller.hpp
@@ -4,6 +4,7 @@
 #include <cmath>
 #include <complex>
 
+#include <astray/math/ipow.hpp>
 #include <astray/math/ode/algebra/quantity_operations.hpp>
 #include <astray/math/ode/error/error_evaluation.hpp>
 #include <astray/math/ode/tableau/tableau_traits.hpp>
@@ -22,7 +23,8 @@ struct integral_controller
     type squared_sum(0);
     operations::for_each([&] (const auto& p, const auto& r, const auto& e)
     {
-      squared_sum += static_cast<type>(std::pow(std::abs(e) / (absolute_tolerance + relative_tolerance * std::max(std::abs(p), std::abs(r))), 2));
+      const auto normalized_error = std::abs(e) / (absolute_tolerance + relative_tolerance * std::max(std::abs(p), std::abs(r)));
+      squared_sum += ipow<2>(normalized_error);
     }, problem.value, result.value, result.error);
 
     type error   = std::sqrt(squared_sum / operations::size(problem.value)); // std::real(squared_sum) unavailable in CUDA until C++20 support.
diff --git a/include/astray/math/ode/error/controller/proportional_integral_controller.hpp b/include/astray/math/ode/error/controller/proportional_integral_controller.hpp
index c3098ac..2906afa 100644
--- a/include/astray/math/ode/error/controller/proportional_integral_controller.hpp
+++ b/include/astray/math/ode/error/controller/proportional_integral_controller.hpp
@@ -4,6 +4,7 @@
 #include <cmath>
 #include <complex>
 
+#include <astray/math/ipow.hpp>
 #include <astray/math/ode/algebra/quantity_operations.hpp>
 #include <astray/math/ode/error/error_evaluation.hpp>
 #include <astray/math/ode/tableau/tableau_traits.hpp>
@@ -22,7 +23,8 @@ struct proportional_integral_controller
     type squared_sum(0);
     operations::for_each([&] (const auto& p, const auto& r, const auto& e)
     {
-      squared_sum += static_cast<type>(std::pow(std::abs(e) / (absolute_tolerance + relative_tolerance * std::max(std::abs(p), std::abs(r))), 2));
+      const auto normalized_error = std::abs(e) / (absolute_tolerance + relative_tolerance * std::max(std::abs(p), std::abs(r)));
+      squared_sum += ipow<2>(normalized_error);
     }, problem.value, result.value, result.error);
 
     type error   = std::sqrt(squared_sum / operations::size(problem.value)); // std::real(squared_sum) unavailable in CUDA until C++20 support.
diff --git a/include/astray/math/ode/error/controller/proportional_integral_derivative_controller.hpp b/include/astray/math/ode/error/controller/proportional_integral_derivative_controller.hpp
index 4b5b5a5..2a39b89 100644
--- a/include/astray/math/ode/error/controller/proportional_integral_derivative_controller.hpp
+++ b/include/astray/math/ode/error/controller/proportional_integral_derivative_controller.hpp
@@ -5,6 +5,7 @@
 #include <cmath>
 #include <complex>
 
+#include <astray/math/ipow.hpp>
 #include <astray/math/ode/algebra/quantity_operations.hpp>
 #include <astray/math/ode/error/error_evaluation.hpp>
 #include <astray/math/ode/tableau/tableau_traits.hpp>
@@ -24,7 +25,8 @@ struct proportional_integral_derivative_controller
     type squared_sum(0);
     operations::for_each([&] (const auto& p, const auto& r, const auto& e)
     {
-      squared_sum += static_cast<type>(std::pow(std::abs(e) / (absolute_tolerance + relative_tolerance * std::max(std::abs(p), std::abs(r))), 2));
+      const auto normalized_error = std::abs(e) / (absolute_tolerance + relative_tolerance * std::max(std::abs(p), std::abs(r)));
+      squared_sum += ipow<2>(normalized_error);
     }, problem.value, result.value, result.error);
 
     error[0]     = type(1) / std::sqrt(squared_sum / operations::size(problem.value)); // std::real(squared_sum) unavailable in CUDA until C++20 support.
diff --git a/include/astray/math/transform.hpp b/include/astray/math/transform.hpp
index d7368d7..bd65b49 100644
--- a/include/astray/math/transform.hpp
+++ b/include/astray/math/transform.hpp
@@ -13,15 +13,15 @@ struct transform
   using quaternion_type = quaternion<scalar_type>;
   using angle_axis_type = angle_axis<scalar_type>;
 
-  vector_type right              () const
+  __device__ __host__ inline vector_type right              () const
   {
     return rotation * vector_type::UnitX();
   }
-  vector_type up                 () const
+  __device__ __host__ inline vector_type up                 () const
   {
     return rotation * vector_type::UnitY();
   }
-  vector_type forward            () const
+  __device__ __host__ inline vector_type forward            () const
   {
     return rotation * vector_type::UnitZ();
   }
diff --git a/include/astray/metrics/boyer_lindquist/kerr.hpp b/include/astray/metrics/boyer_lindquist/kerr.hpp
index 4726d86..1d474cb 100644
--- a/include/astray/metrics/boyer_lindquist/kerr.hpp
+++ b/include/astray/metrics/boyer_lindquist/kerr.hpp
@@ -4,6 +4,7 @@
 
 #include <astray/core/metric.hpp>
 #include <astray/math/constants.hpp>
+#include <astray/math/ipow.hpp>
 
 namespace ast::metrics
 {
@@ -22,7 +23,9 @@ class kerr : public metric<coordinate_system_type::boyer_lindquist, scalar_type,
   }
   __device__ termination_reason       check_termination          (const vector_type& position, const vector_type& direction) const override
   {
-    const auto event_horizon = mass + std::sqrt(static_cast<scalar_type>(std::pow(mass, 2)) - static_cast<scalar_type>(std::pow(angular_momentum, 2)));
+    const auto mass_sq = ipow<2>(mass);
+    const auto ang_mom_sq = ipow<2>(angular_momentum);
+    const auto event_horizon = mass + std::sqrt(mass_sq - ang_mom_sq);
     if (position[1] < static_cast<scalar_type>(0) || position[1] <= (static_cast<scalar_type>(1) + consts::epsilon) * event_horizon)
       return termination_reason::spacetime_breakdown;
     return termination_reason::none;
@@ -30,64 +33,64 @@ class kerr : public metric<coordinate_system_type::boyer_lindquist, scalar_type,
 
   __device__ christoffel_symbols_type christoffel_symbols        (const vector_type& position) const override
   {
-    const auto t1   = static_cast<scalar_type>(std::pow(position[1], 2));
+    const auto r_sq = ipow<2>(position[1]);
     const auto t2   = mass * position[1];
-    const auto t4   = static_cast<scalar_type>(std::pow(angular_momentum, 2));
-    const auto t5   = t1 - static_cast<scalar_type>(2) * t2 + t4;
+    const auto a_sq = ipow<2>(angular_momentum);
+    const auto t5   = r_sq - static_cast<scalar_type>(2) * t2 + a_sq;
     const auto t6   = std::cos(position[2]);
-    const auto t7   = static_cast<scalar_type>(std::pow(t6, 2));
-    const auto t8   = t4 * t7;
-    const auto t9   = t1 + t8;
-    const auto t10  = static_cast<scalar_type>(std::pow(t9, 2));
-    const auto t12  = static_cast<scalar_type>(1) / t10 / t9;
-    const auto t14  = -t1 + t8;
+    const auto t6_sq = ipow<2>(t6);
+    const auto t8   = a_sq * t6_sq;
+    const auto t9   = r_sq + t8;
+    const auto t9_sq = ipow<2>(t9);
+    const auto t12  = static_cast<scalar_type>(1) / t9_sq / t9;
+    const auto t14  = -r_sq + t8;
     const auto t20  = std::sin(position[2]);
-    const auto t21  = t4 * t6 * t20;
-    const auto t24  = t4 + t1;
+    const auto t21  = a_sq * t6 * t20;
+    const auto t24  = a_sq + r_sq;
     const auto t26  = static_cast<scalar_type>(1) / t9;
-    const auto t28  = static_cast<scalar_type>(std::pow(t4, 2));
-    const auto t29  = t28 * t7;
-    const auto t30  = t1 * t4;
-    const auto t31  = t30 * t7;
-    const auto t32  = t4 * mass;
-    const auto t35  = static_cast<scalar_type>(2) * t32 * position[1] * t7;
-    const auto t36  = static_cast<scalar_type>(std::pow(t1, 2));
-    const auto t37  = t1 * position[1];
+    const auto a_sq_sq = ipow<2>(a_sq);
+    const auto t29  = a_sq_sq * t6_sq;
+    const auto t30  = r_sq * a_sq;
+    const auto t31  = t30 * t6_sq;
+    const auto t32  = a_sq * mass;
+    const auto t35  = static_cast<scalar_type>(2) * t32 * position[1] * t6_sq;
+    const auto r_sq_sq = ipow<2>(r_sq);
+    const auto t37  = r_sq * position[1];
     const auto t38  = mass * t37;
-    const auto t41  = static_cast<scalar_type>(1) / (t29 + t31 + t30 - t35 + t36 - static_cast<scalar_type>(2) * t38);
+    const auto t41  = static_cast<scalar_type>(1) / (t29 + t31 + t30 - t35 + r_sq_sq - static_cast<scalar_type>(2) * t38);
     const auto t42  = t14 * t26 * t41;
     const auto t43  = t24 * mass * t42;
     const auto t44  = mass * angular_momentum;
     const auto t45  = t44 * t42;
-    const auto t46  = static_cast<scalar_type>(1) / t10;
+    const auto t46  = static_cast<scalar_type>(1) / t9_sq;
     const auto t49  = static_cast<scalar_type>(2) * t2 * t46 * t21;
     const auto t50  = t44 * position[1];
     const auto t51  = static_cast<scalar_type>(1) / t20;
     const auto t55  = static_cast<scalar_type>(2) * t50 * t6 * t51 * t46;
-    const auto t58  = -static_cast<scalar_type>(1) + t7;
+    const auto t58  = -static_cast<scalar_type>(1) + t6_sq;
     const auto t61  = t5 * mass * angular_momentum * t58 * t14 * t12;
     const auto t62  = t20 * t6;
     const auto t66  = static_cast<scalar_type>(2) * t50 * t62 * t24 * t12;
     const auto t68  = static_cast<scalar_type>(1) / t5 * t26;
-    const auto t69  = mass * t1;
-    const auto t77  = t26 * t4 * t62;
+    const auto t69  = mass * r_sq;
+    const auto t77  = t26 * a_sq * t62;
     const auto t78  = t26 * position[1];
-    const auto t85  = (t29 - t31 - t30 - static_cast<scalar_type>(3) * t36) * mass * angular_momentum * t58 * t26 * t41;
-    const auto t87  = static_cast<scalar_type>(std::pow(t7, 2));
-    const auto t88  = position[1] * t28 * t87;
-    const auto t89  = mass * t28;
-    const auto t90  = t89 * t7;
-    const auto t91  = t89 * t87;
+    const auto t85  = (t29 - t31 - t30 - static_cast<scalar_type>(3) * r_sq_sq) * mass * angular_momentum * t58 * t26 * t41;
+    const auto t6_sq_sq = ipow<2>(t6_sq);
+    const auto t88  = position[1] * a_sq_sq * t6_sq_sq;
+    const auto t89  = mass * a_sq_sq;
+    const auto t90  = t89 * t6_sq;
+    const auto t91  = t89 * t6_sq_sq;
     const auto t92  = t69 * t8;
-    const auto t95  = static_cast<scalar_type>(2) * t37 * t4 * t7;
-    const auto t96  = t32 * t1;
-    const auto t97  = t36 * position[1];
-    const auto t102 = (t88 + t90 - t91 - t92 + t95 - t96 + t97 - static_cast<scalar_type>(2) * mass * t36) * t26 * t41;
-    const auto t112 = static_cast<scalar_type>(2) * t20 * t58 * mass * t4 * angular_momentum * position[1] * t6 * t46;
-    const auto t113 = t87 * t28;
-    const auto t120 = t6 * (t113 + static_cast<scalar_type>(2) * t32 * position[1] - t35 + static_cast<scalar_type>(2) * t31 + t36) * t51 * t46;
-    const auto t126 = t36 * t4;
-    const auto t129 = t1 * t28;
+    const auto t95  = static_cast<scalar_type>(2) * t37 * a_sq * t6_sq;
+    const auto t96  = t32 * r_sq;
+    const auto t97  = r_sq_sq * position[1];
+    const auto t102 = (t88 + t90 - t91 - t92 + t95 - t96 + t97 - static_cast<scalar_type>(2) * mass * r_sq_sq) * t26 * t41;
+    const auto t112 = static_cast<scalar_type>(2) * t20 * t58 * mass * a_sq * angular_momentum * position[1] * t6 * t46;
+    const auto t113 = t6_sq_sq * a_sq_sq;
+    const auto t120 = t6 * (t113 + static_cast<scalar_type>(2) * t32 * position[1] - t35 + static_cast<scalar_type>(2) * t31 + r_sq_sq) * t51 * t46;
+    const auto t126 = r_sq_sq * a_sq;
+    const auto t129 = r_sq * a_sq_sq;
 
     christoffel_symbols_type symbols;
     symbols.setZero();
@@ -101,7 +104,7 @@ class kerr : public metric<coordinate_system_type::boyer_lindquist, scalar_type,
     symbols(0, 3, 2) =  t66;
     symbols(1, 0, 0) = -t43;
     symbols(1, 0, 3) = -t45;
-    symbols(1, 1, 1) = -t68 * (t69 - position[1] * t4 + t8 * position[1] - t8 * mass);
+    symbols(1, 1, 1) = -t68 * (t69 - position[1] * a_sq + t8 * position[1] - t8 * mass);
     symbols(1, 1, 2) =  t68 * t21;
     symbols(1, 2, 1) = -t77;
     symbols(1, 2, 2) =  t78;
@@ -122,9 +125,9 @@ class kerr : public metric<coordinate_system_type::boyer_lindquist, scalar_type,
     symbols(3, 2, 0) = -t112;
     symbols(3, 2, 3) =  t120;
     symbols(3, 3, 1) =  t5 * t58 * (t97 + t95 + t88 - t96 + t92 + t90 - t91) * t12;
-    symbols(3, 3, 2) = -t62 * (t36 * t1 + static_cast<scalar_type>(2) * t126 * t7 + t129 * t87 + t126
-      + static_cast<scalar_type>(2) * t129 * t7   + t28 * t4 * t87 
-      + static_cast<scalar_type>(4) * t38  * t4   - static_cast<scalar_type>(4) * t38 * t8 
+    symbols(3, 3, 2) = -t62 * (r_sq_sq * r_sq + static_cast<scalar_type>(2) * t126 * t6_sq + t129 * t6_sq_sq + t126
+      + static_cast<scalar_type>(2) * t129 * t6_sq   + a_sq_sq * a_sq * t6_sq_sq 
+      + static_cast<scalar_type>(4) * t38  * a_sq   - static_cast<scalar_type>(4) * t38 * t8 
       - static_cast<scalar_type>(2) * t2   * t113 + static_cast<scalar_type>(2) * t89 * position[1]) * t12;
     return symbols;
   }
diff --git a/include/astray/metrics/cartesian/alcubierre.hpp b/include/astray/metrics/cartesian/alcubierre.hpp
index 02741ab..252a1d4 100644
--- a/include/astray/metrics/cartesian/alcubierre.hpp
+++ b/include/astray/metrics/cartesian/alcubierre.hpp
@@ -4,6 +4,7 @@
 
 #include <astray/core/metric.hpp>
 #include <astray/math/constants.hpp>
+#include <astray/math/ipow.hpp>
 
 namespace ast::metrics
 {
@@ -19,51 +20,47 @@ class alcubierre : public metric<coordinate_system_type::cartesian, scalar_type,
   __device__ christoffel_symbols_type christoffel_symbols(const vector_type& position) const override
   {
     const auto xmvt          = position[1] - velocity * position[0];
-    const auto rs            = static_cast<scalar_type>(std::sqrt(
-      std::pow(xmvt       , 2) +
-      std::pow(position[2], 2) + 
-      std::pow(position[3], 2)));
+    const auto xmvt_sq       = ipow<2>(xmvt);
+    const auto y_sq          = ipow<2>(position[2]);
+    const auto z_sq          = ipow<2>(position[3]);
+    const auto rs            = std::sqrt(xmvt_sq + y_sq + z_sq);
 
-    const auto w1            = std::tanh(thickness * (rs + radius));
-    const auto w2            = std::tanh(thickness * (rs - radius));
-    const auto w3            = std::tanh(thickness * radius);
-    const auto f             = static_cast<scalar_type>(0.5) * (w1 / w3 - w2 / w3);
-    
     const auto tanh_positive = std::tanh(thickness * (rs + radius));
     const auto tanh_negative = std::tanh(thickness * (rs - radius));
-    const auto factor        = static_cast<scalar_type>(
-      (std::pow(tanh_negative, 2) - std::pow(tanh_positive, 2)) * 
-      (static_cast<scalar_type>(0.5) * thickness / (rs * std::tanh(thickness * radius))));
+    const auto w3            = std::tanh(thickness * radius);
+    const auto f             = static_cast<scalar_type>(0.5) * (tanh_positive / w3 - tanh_negative / w3);
+    
+    const auto tanh_pos_sq   = ipow<2>(tanh_positive);
+    const auto tanh_neg_sq   = ipow<2>(tanh_negative);
+    const auto factor        = (tanh_neg_sq - tanh_pos_sq) * 
+      (static_cast<scalar_type>(0.5) * thickness / (rs * w3));
     vector_type df {
       -velocity * xmvt * factor,
                   xmvt * factor,
       position[2]      * factor,
       position[3]      * factor};
 
-    const auto t1  = std::pow(velocity, 2);
-    const auto t2  = t1 * velocity;
-    const auto t3  = f;
-    const auto t4  = std::pow(t3, 2);
+    const auto v_sq = ipow<2>(velocity);
+    const auto t2  = v_sq * velocity;
+    const auto f_sq = ipow<2>(f);
     const auto t6  = df[1];
-    const auto t7  = std::pow(consts::speed_of_light, 2);
-    const auto t8  = static_cast<scalar_type>(1) / t7;
-    const auto t10 = t2 * t4 * t6 * t8;
+    const auto t8  = static_cast<scalar_type>(1) / consts::speed_of_light_squared;
+    const auto t10 = t2 * f_sq * t6 * t8;
     const auto t11 = df[0];
-    const auto t14 = t3 * t6;
-    const auto t22 = t1 * t3;
+    const auto t14 = f * t6;
+    const auto t22 = v_sq * f;
     const auto t23 = df[2];
     const auto t25 = df[3];
-    const auto t27 = t8 * t1;
+    const auto t27 = t8 * v_sq;
     const auto t28 = t27 * t14;
-    const auto t29 = velocity * t23;
-    const auto t30 = t29                  / static_cast<scalar_type>(2);
+    const auto t30 = velocity * t23       / static_cast<scalar_type>(2);
     const auto t31 = velocity * t25;
     const auto t32 = t31                  / static_cast<scalar_type>(2);
-    const auto t35 = t27 * t3 * t23       / static_cast<scalar_type>(2);
-    const auto t38 = (t1 * t4 + t7) * t8;
-    const auto t40 = t29 * t38            / static_cast<scalar_type>(2);
-    const auto t43 = t27 * t3 * t25       / static_cast<scalar_type>(2);
-    const auto t45 = t31 * t38            / static_cast<scalar_type>(2);
+    const auto t35 = t27 * f * t23       / static_cast<scalar_type>(2);
+    const auto t38 = (v_sq * f_sq + consts::speed_of_light_squared) * t8;
+    const auto t40 = t30 * t38;
+    const auto t43 = t27 * f * t25       / static_cast<scalar_type>(2);
+    const auto t45 = t32 * t38;
     const auto t46 = velocity * t8;
     const auto t49 = t46 * t23            / static_cast<scalar_type>(2);
     const auto t51 = t46 * t25            / static_cast<scalar_type>(2);
@@ -71,7 +68,7 @@ class alcubierre : public metric<coordinate_system_type::cartesian, scalar_type,
     christoffel_symbols_type symbols;
     symbols.setZero();
     symbols(0, 0, 0) =  t10;
-    symbols(0, 0, 1) =  velocity * (-t7 * t11 - t7 * velocity * t14 + t2 * t4 * t3 * t6) * t8;
+    symbols(0, 0, 1) =  velocity * (-consts::speed_of_light_squared * t11 - consts::speed_of_light_squared * velocity * t14 + t2 * f_sq * f * t6) * t8;
     symbols(0, 0, 2) = -t22 * t23;
     symbols(0, 0, 3) = -t22 * t25;
     symbols(0, 1, 0) = -t28;
diff --git a/include/astray/metrics/spherical/kottler.hpp b/include/astray/metrics/spherical/kottler.hpp
index e52ce71..784659f 100644
--- a/include/astray/metrics/spherical/kottler.hpp
+++ b/include/astray/metrics/spherical/kottler.hpp
@@ -4,6 +4,7 @@
 
 #include <astray/core/metric.hpp>
 #include <astray/math/constants.hpp>
+#include <astray/math/ipow.hpp>
 
 namespace ast::metrics
 {
@@ -19,8 +20,12 @@ class kottler : public metric<coordinate_system_type::spherical, scalar_type, ve
   __device__ termination_reason       check_termination  (const vector_type& position, const vector_type& direction) const override
   {
     const auto rs = consts::schwarzschild_radius(mass);
-    if (position[1] < static_cast<scalar_type>(0) || std::abs(static_cast<scalar_type>(1) - rs / position[1] 
-      - consts::cosmological_constant / static_cast<scalar_type>(3) * std::pow(position[1], 2)) <= consts::epsilon)
+    const auto r = position[1];
+    if (r < static_cast<scalar_type>(0))
+      return termination_reason::spacetime_breakdown;
+    const auto r_sq = ipow<2>(r);
+    if (std::abs(static_cast<scalar_type>(1) - rs / r 
+      - consts::cosmological_constant / static_cast<scalar_type>(3) * r_sq) <= consts::epsilon)
       return termination_reason::spacetime_breakdown;
     return termination_reason::none;
   }
@@ -28,19 +33,20 @@ class kottler : public metric<coordinate_system_type::spherical, scalar_type, ve
   __device__ christoffel_symbols_type christoffel_symbols(const vector_type& position) const override
   {
     const auto rs  = consts::schwarzschild_radius(mass);
+    const auto r   = position[1];
+    const auto r_sq = ipow<2>(r);
 
     const auto t2  = static_cast<scalar_type>(3) * rs;
-    const auto t3  = static_cast<scalar_type>(std::pow(position[1], 2));
-    const auto t4  = t3 * position[1];
+    const auto t4  = r_sq * r;
     const auto t5  = consts::cosmological_constant * t4;
-    const auto t6  = -static_cast<scalar_type>(3) * position[1] + t2 + t5;
+    const auto t6  = -static_cast<scalar_type>(3) * r + t2 + t5;
     const auto t11 = -t2 + static_cast<scalar_type>(2) * t5;
-    const auto t15 = static_cast<scalar_type>(1) / position[1];
+    const auto t15 = static_cast<scalar_type>(1) / r;
     const auto t19 = t15 / t6 * t11 / static_cast<scalar_type>(2);
     const auto t22 = std::sin(position[2]);
     const auto t24 = std::cos(position[2]);
     const auto t25 = static_cast<scalar_type>(1) / t22 * t24;
-    const auto t26 = static_cast<scalar_type>(std::pow(t22, 2));
+    const auto t22_sq = ipow<2>(t22);
 
     christoffel_symbols_type symbols;
     symbols.setZero();
@@ -51,11 +57,11 @@ class kottler : public metric<coordinate_system_type::spherical, scalar_type, ve
     symbols(1, 2, 2) =  t15;
     symbols(1, 3, 3) =  t15;
     symbols(2, 1, 2) =  t15;
-    symbols(2, 2, 1) = -position[1] + rs + t5 / static_cast<scalar_type>(3);
+    symbols(2, 2, 1) = -r + rs + t5 / static_cast<scalar_type>(3);
     symbols(2, 3, 3) =  t25;
     symbols(3, 1, 3) =  t15;
     symbols(3, 2, 3) =  t25;
-    symbols(3, 3, 1) =  t6  * t26 / static_cast<scalar_type>(3);
+    symbols(3, 3, 1) =  t6  * t22_sq / static_cast<scalar_type>(3);
     symbols(3, 3, 2) = -t22 * t24;
     return symbols;
   }
diff --git a/include/astray/metrics/spherical/morris_thorne.hpp b/include/astray/metrics/spherical/morris_thorne.hpp
index 2e524ea..961c109 100644
--- a/include/astray/metrics/spherical/morris_thorne.hpp
+++ b/include/astray/metrics/spherical/morris_thorne.hpp
@@ -3,6 +3,7 @@
 #include <cmath>
 
 #include <astray/core/metric.hpp>
+#include <astray/math/ipow.hpp>
 
 namespace ast::metrics
 {
@@ -15,24 +16,25 @@ class morris_thorne : public metric<coordinate_system_type::spherical, scalar_ty
 public:
   __device__ christoffel_symbols_type christoffel_symbols(const vector_type& position) const override
   {
-    const auto t1  = static_cast<scalar_type>(std::pow(position[1]  , 2));
-    const auto t2  = static_cast<scalar_type>(std::pow(throat_radius, 2));
-    const auto t5  = static_cast<scalar_type>(1) / (t1 + t2) * position[1];
+    const auto r   = position[1];
+    const auto r_sq = ipow<2>(r);
+    const auto b_sq = ipow<2>(throat_radius);
+    const auto t5  = static_cast<scalar_type>(1) / (r_sq + b_sq) * r;
     const auto t6  = std::sin(position[2]);
     const auto t8  = std::cos(position[2]);
     const auto t9  = static_cast<scalar_type>(1) / t6 * t8;
-    const auto t10 = static_cast<scalar_type>(std::pow(t6, 2));
+    const auto t6_sq = ipow<2>(t6);
 
     christoffel_symbols_type symbols;
     symbols.setZero();
     symbols(1, 2, 2) = t5;
     symbols(1, 3, 3) = t5;
     symbols(2, 1, 2) = t5;
-    symbols(2, 2, 1) = -position[1];
+    symbols(2, 2, 1) = -r;
     symbols(2, 3, 3) = t9;
     symbols(3, 1, 3) = t5;
     symbols(3, 2, 3) = t9;
-    symbols(3, 3, 1) = -position[1] * t10;
+    symbols(3, 3, 1) = -r * t6_sq;
     symbols(3, 3, 2) = -t6 * t8;
     return symbols;
   }
diff --git a/include/astray/metrics/spherical/reissner_nordstroem.hpp b/include/astray/metrics/spherical/reissner_nordstroem.hpp
index 3afe8a5..5c8d81d 100644
--- a/include/astray/metrics/spherical/reissner_nordstroem.hpp
+++ b/include/astray/metrics/spherical/reissner_nordstroem.hpp
@@ -4,6 +4,7 @@
 
 #include <astray/core/metric.hpp>
 #include <astray/math/constants.hpp>
+#include <astray/math/ipow.hpp>
 
 namespace ast::metrics
 {
@@ -36,25 +37,24 @@ class reissner_nordstroem : public metric<coordinate_system_type::spherical, sca
 
   __device__ christoffel_symbols_type christoffel_symbols(const vector_type& position) const override
   {
-    const auto t1  = std::pow(position[1], 2);
-    const auto t2  = consts::schwarzschild_radius(mass) * position[1];
-    const auto t3  = std::pow(charge, 2);
+    const auto r   = position[1];
+    const auto r_sq = ipow<2>(r);
+    const auto t2  = consts::schwarzschild_radius(mass) * r;
     const auto t4  = consts::characteristic_length_scale(charge);
-    const auto t5  = t1 - t2 + t4;
-    const auto t6  = std::pow(t1, 2);
-    const auto t10 = consts::speed_of_light_squared;
+    const auto t5  = r_sq - t2 + t4;
+    const auto r_sq_sq = ipow<2>(r_sq);
     const auto t12 = t2 - static_cast<scalar_type>(2) * t4;
-    const auto t16 = static_cast<scalar_type>(1) / position[1];
+    const auto t16 = static_cast<scalar_type>(1) / r;
     const auto t20 = t16 / t5 * t12 / static_cast<scalar_type>(2);
     const auto t21 = t5 * t16;
     const auto t22 = std::sin(position[2]);
     const auto t24 = std::cos(position[2]);
     const auto t25 = static_cast<scalar_type>(1) / t22 * t24;
-    const auto t26 = std::pow(t22, 2);
+    const auto t22_sq = ipow<2>(t22);
 
     christoffel_symbols_type symbols;
     symbols.setZero();
-    symbols(0, 0, 1) =  t5 / t6 / position[1] * t10 * t12 / static_cast<scalar_type>(2);
+    symbols(0, 0, 1) =  t5 / r_sq_sq / r * consts::speed_of_light_squared * t12 / static_cast<scalar_type>(2);
     symbols(0, 1, 0) =  t20;
     symbols(1, 0, 0) =  t20;
     symbols(1, 1, 1) = -t20;
@@ -65,7 +65,7 @@ class reissner_nordstroem : public metric<coordinate_system_type::spherical, sca
     symbols(2, 3, 3) =  t25;
     symbols(3, 1, 3) =  t16;
     symbols(3, 2, 3) =  t25;
-    symbols(3, 3, 1) = -t21 * t26;
+    symbols(3, 3, 1) = -t21 * t22_sq;
     symbols(3, 3, 2) = -t22 * t24;
     return symbols;
   }
diff --git a/include/astray/metrics/spherical/schwarzschild.hpp b/include/astray/metrics/spherical/schwarzschild.hpp
index a554f00..af12021 100644
--- a/include/astray/metrics/spherical/schwarzschild.hpp
+++ b/include/astray/metrics/spherical/schwarzschild.hpp
@@ -4,6 +4,7 @@
 
 #include <astray/core/metric.hpp>
 #include <astray/math/constants.hpp>
+#include <astray/math/ipow.hpp>
 
 namespace ast::metrics
 {
@@ -19,8 +20,12 @@ class schwarzschild : public metric<coordinate_system_type::spherical, scalar_ty
   __device__ termination_reason       check_termination  (const vector_type& position, const vector_type& direction) const override
   {
     const auto rs = consts::schwarzschild_radius(mass);
-    if (position[1] < static_cast<scalar_type>(0) || 
-        static_cast<scalar_type>(std::pow(position[1], 2)) <= (static_cast<scalar_type>(1) + consts::epsilon) * static_cast<scalar_type>(std::pow(rs, 2)))
+    const auto r = position[1];
+    if (r < static_cast<scalar_type>(0))
+      return termination_reason::spacetime_breakdown;
+    const auto r_sq = ipow<2>(r);
+    const auto rs_sq = ipow<2>(rs);
+    if (r_sq <= (static_cast<scalar_type>(1) + consts::epsilon) * rs_sq)
       return termination_reason::spacetime_breakdown;
     return termination_reason::none;
   }
@@ -31,18 +36,18 @@ class schwarzschild : public metric<coordinate_system_type::spherical, scalar_ty
     const auto r     = position[1];
     const auto theta = position[2];
     const auto t1    = r - rs;
-    const auto t2    = static_cast<scalar_type>(std::pow(r, 2));
+    const auto r_sq  = ipow<2>(r);
     const auto t6    = consts::speed_of_light_squared;
     const auto t10   = static_cast<scalar_type>(1) / r;
     const auto t14   = t10 / t1 * rs * static_cast<scalar_type>(0.5);
     const auto t15   = std::sin(theta);
     const auto t17   = std::cos(theta);
     const auto t18   = static_cast<scalar_type>(1) / t15 * t17;
-    const auto t19   = static_cast<scalar_type>(std::pow(t15, 2));
+    const auto t15_sq = ipow<2>(t15);
 
     christoffel_symbols_type symbols;
     symbols.setZero();
-    symbols(0, 0, 1) = t1 / t2 / r * t6 * rs / static_cast<scalar_type>(2);
+    symbols(0, 0, 1) = t1 / r_sq / r * t6 * rs / static_cast<scalar_type>(2);
     symbols(0, 1, 0) =  t14;
     symbols(1, 0, 0) =  t14;
     symbols(1, 1, 1) = -t14;
@@ -53,7 +58,7 @@ class schwarzschild : public metric<coordinate_system_type::spherical, scalar_ty
     symbols(2, 3, 3) =  t18;
     symbols(3, 1, 3) =  t10;
     symbols(3, 2, 3) =  t18;
-    symbols(3, 3, 1) = -t1  * t19;
+    symbols(3, 3, 1) = -t1  * t15_sq;
     symbols(3, 3, 2) = -t15 * t17;
     return symbols;
   }