Added simd::pow(T, int)

include/dsp/ode.hpp

@@ -8,7 +8,7 @@ namespace dsp {
 
 /** The callback function `f` in each of these stepping functions must have the signature
 
-	void f(float t, const float x[], float dxdt[])
+	void f(T t, const T x[], T dxdt[])
 
 A capturing lambda is ideal for this.
 For example, the following solves the system x''(t) = -x(t) using a fixed timestep of 0.01 and initial conditions x(0) = 1, x'(0) = 0.

include/simd/functions.hpp

@@ -121,6 +121,18 @@ inline float_4 pow(float a, float_4 b) {
 	return exp(b * std::log(a));
 }
 
+template <typename T>
+T pow(T a, int b) {
+	// Optimal with `-O3 -ffast-math` when b is known at compile-time
+	T p = 1;
+	for (int i = 1; i <= b; i *= 2) {
+		if (i & b)
+			p *= a;
+		a *= a;
+	}
+	return p;
+}
+
 // Nonstandard functions
 
 inline float ifelse(bool cond, float a, float b) {