diff --git a/numerical/analysis/ConvolutionCorrelation.hpp b/numerical/analysis/ConvolutionCorrelation.hpp
index 550ac9b5..b9ad0662 100644
--- a/numerical/analysis/ConvolutionCorrelation.hpp
+++ b/numerical/analysis/ConvolutionCorrelation.hpp
@@ -127,11 +127,8 @@ namespace analysis
         for (std::size_t i = 0; i < M + K - 1; ++i)
             y[i] = yFull[i];
     }
-}
 
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
-namespace analysis
-{
     extern template void LinearConvolution<float, 3, 3>(
         const infra::BoundedVector<float>::WithMaxSize<3>&,
         const infra::BoundedVector<float>::WithMaxSize<3>&,
@@ -159,5 +156,5 @@ namespace analysis
         const infra::BoundedVector<float>::WithMaxSize<3>&,
         infra::BoundedVector<float>::WithMaxSize<5>&,
         FastFourierTransform<float>&);
-}
 #endif
+}
diff --git a/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp b/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp
index 046a3fa3..a428ecc3 100644
--- a/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp
+++ b/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp
@@ -101,6 +101,13 @@ namespace
     };
 }
 
+TYPED_TEST(TestFastFourierTransform, log2_runtime_both_branches)
+{
+    auto& fftInst = *this->fft;
+    EXPECT_EQ(fftInst.Log2(1), 0u);
+    EXPECT_EQ(fftInst.Log2(8), 3u);
+}
+
 TYPED_TEST(TestFastFourierTransform, zero_input_produces_zero_output)
 {
     this->timeDomain.clear();
diff --git a/numerical/analysis/windowing/test/TestWindowing.cpp b/numerical/analysis/windowing/test/TestWindowing.cpp
index 3a1f1af1..d6db172b 100644
--- a/numerical/analysis/windowing/test/TestWindowing.cpp
+++ b/numerical/analysis/windowing/test/TestWindowing.cpp
@@ -123,3 +123,27 @@ TYPED_TEST(WindowingTest, WindowSymmetry)
         }
     }
 }
+
+TYPED_TEST(WindowingTest, HammingWindowPower)
+{
+    windowing::HammingWindow<TypeParam> w;
+    EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.397f, this->kEpsilon);
+}
+
+TYPED_TEST(WindowingTest, HanningWindowPower)
+{
+    windowing::HanningWindow<TypeParam> w;
+    EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.375f, this->kEpsilon);
+}
+
+TYPED_TEST(WindowingTest, BlackmanWindowPower)
+{
+    windowing::BlackmanWindow<TypeParam> w;
+    EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.305f, this->kEpsilon);
+}
+
+TYPED_TEST(WindowingTest, RectangularWindowPower)
+{
+    windowing::RectangularWindow<TypeParam> w;
+    EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.999f, this->kEpsilon);
+}
diff --git a/numerical/control_analysis/ControllabilityObservability.hpp b/numerical/control_analysis/ControllabilityObservability.hpp
index 685f8e28..f0dab178 100644
--- a/numerical/control_analysis/ControllabilityObservability.hpp
+++ b/numerical/control_analysis/ControllabilityObservability.hpp
@@ -131,8 +131,11 @@ namespace control_analysis
         std::size_t rank{ 0 };
         std::size_t pivotRow{ 0 };
 
-        for (std::size_t col = 0; col < Cols && pivotRow < Rows; ++col)
+        for (std::size_t col = 0; col < Cols; ++col)
         {
+            if (pivotRow >= Rows)
+                break;
+
             std::size_t maxRow{ pivotRow };
             T maxVal{ T(0) };
             for (std::size_t r = pivotRow; r < Rows; ++r)
diff --git a/numerical/control_analysis/test/TestFrequencyResponse.cpp b/numerical/control_analysis/test/TestFrequencyResponse.cpp
index bc6a85c3..4f156f71 100644
--- a/numerical/control_analysis/test/TestFrequencyResponse.cpp
+++ b/numerical/control_analysis/test/TestFrequencyResponse.cpp
@@ -317,3 +317,18 @@ TEST_F(TestFrequencyResponseHighpass, magnitude_at_quarter_nyquist_is_minus_3db)
 
     EXPECT_NEAR(bestMag, -3.0103f, 0.5f);
 }
+
+TEST_F(TestFrequencyResponseUnity, zero_denominator_coefficients_produce_finite_output)
+{
+    std::array<float, 1> bz{ 1.0f };
+    std::array<float, 1> az{ 0.0f };
+    control_analysis::FrequencyResponse<float, 64> frZeroDenom{ bz, az, kSampleFrequency };
+
+    auto [frequencies, magnitudes, phases] = frZeroDenom.Calculate();
+
+    for (const auto& m : magnitudes)
+    {
+        EXPECT_FALSE(std::isnan(m));
+        EXPECT_FALSE(std::isinf(m));
+    }
+}
diff --git a/numerical/controllers/implementations/GainScheduledController.hpp b/numerical/controllers/implementations/GainScheduledController.hpp
index 6a2c33c2..49b396b1 100644
--- a/numerical/controllers/implementations/GainScheduledController.hpp
+++ b/numerical/controllers/implementations/GainScheduledController.hpp
@@ -35,13 +35,12 @@ namespace controllers
 
     private:
         std::array<SchedulePoint<T, GainSize>, N> table;
-        std::array<T, GainSize> active;
+        std::array<T, GainSize> active{};
     };
 
     template<typename T, std::size_t N, std::size_t GainSize>
     GainScheduledController<T, N, GainSize>::GainScheduledController(std::array<SchedulePoint<T, GainSize>, N> scheduleTable)
         : table{ scheduleTable }
-        , active{}
     {
         for (std::size_t i{ 0 }; i < N - 1; ++i)
             assert(table[i].breakpoint < table[i + 1].breakpoint);
diff --git a/numerical/controllers/implementations/LuenbergerObserver.hpp b/numerical/controllers/implementations/LuenbergerObserver.hpp
index e8c474cf..fde6ccc9 100644
--- a/numerical/controllers/implementations/LuenbergerObserver.hpp
+++ b/numerical/controllers/implementations/LuenbergerObserver.hpp
@@ -45,7 +45,7 @@ namespace controllers
 
         Plant plant;
         GainMatrix L;
-        StateVector xhat;
+        StateVector xhat{};
     };
 
     template<typename T, std::size_t StateSize, std::size_t InputSize, std::size_t OutputSize>
@@ -53,7 +53,6 @@ namespace controllers
         const Plant& plantModel, const GainMatrix& observerGain)
         : plant{ plantModel }
         , L{ observerGain }
-        , xhat{}
     {}
 
     template<typename T, std::size_t StateSize, std::size_t InputSize, std::size_t OutputSize>
diff --git a/numerical/controllers/implementations/test/TestIntegralStateFeedbackLqi.cpp b/numerical/controllers/implementations/test/TestIntegralStateFeedbackLqi.cpp
index 127bad11..e3ffd8bf 100644
--- a/numerical/controllers/implementations/test/TestIntegralStateFeedbackLqi.cpp
+++ b/numerical/controllers/implementations/test/TestIntegralStateFeedbackLqi.cpp
@@ -195,3 +195,32 @@ TEST_F(TestIntegralStateFeedbackLqi, control_uses_negative_feedback)
     EXPECT_NEAR(u.at(0, 0), expected, math::Tolerance<float>());
     EXPECT_LT(u.at(0, 0), 0.0f);
 }
+
+TEST_F(TestIntegralStateFeedbackLqi, direct_gain_constructor_matches_lqr_constructor)
+{
+    const auto kx = controller.GetGainState();
+    const auto ki = controller.GetGainIntegral();
+
+    Controller directController{ kx, ki, 0.01f };
+
+    math::Vector<float, 2> x{ { 0.5f }, { 0.1f } };
+    math::Vector<float, 1> reference{ { 1.0f } };
+    math::Vector<float, 1> measured{ { 0.3f } };
+
+    auto u1 = controller.ComputeControl(x, reference, measured);
+    auto u2 = directController.ComputeControl(x, reference, measured);
+
+    EXPECT_NEAR(u1.at(0, 0), u2.at(0, 0), math::Tolerance<float>());
+}
+
+TEST_F(TestIntegralStateFeedbackLqi, direct_gain_constructor_get_gains)
+{
+    math::Matrix<float, 1, 2> kx{ { 1.5f, 0.8f } };
+    math::Matrix<float, 1, 1> ki{ { 2.0f } };
+
+    Controller c{ kx, ki, 0.01f };
+
+    EXPECT_FLOAT_EQ(c.GetGainState().at(0, 0), 1.5f);
+    EXPECT_FLOAT_EQ(c.GetGainState().at(0, 1), 0.8f);
+    EXPECT_FLOAT_EQ(c.GetGainIntegral().at(0, 0), 2.0f);
+}
diff --git a/numerical/estimators/offline/ExpectationMaximization.hpp b/numerical/estimators/offline/ExpectationMaximization.hpp
index 1ca798ca..7e76f90d 100644
--- a/numerical/estimators/offline/ExpectationMaximization.hpp
+++ b/numerical/estimators/offline/ExpectationMaximization.hpp
@@ -104,8 +104,7 @@ namespace estimators
         while (iter < maxIterations)
         {
             const auto smootherOutput = smoother_.Smooth(
-                currentParams.F, currentParams.H,
-                currentParams.Q, currentParams.R,
+                { currentParams.F, currentParams.H, currentParams.Q, currentParams.R },
                 observations, numSteps,
                 currentParams.initialState, currentParams.initialCovariance);
 
diff --git a/numerical/estimators/offline/test/TestPolynomialFitting.cpp b/numerical/estimators/offline/test/TestPolynomialFitting.cpp
index 358e4afa..a2f027a4 100644
--- a/numerical/estimators/offline/test/TestPolynomialFitting.cpp
+++ b/numerical/estimators/offline/test/TestPolynomialFitting.cpp
@@ -180,6 +180,24 @@ TEST_F(TestPolynomialFitting, determinism_identical_fit_produces_identical_coeff
     EXPECT_FLOAT_EQ(c1.at(2, 0), c2.at(2, 0));
 }
 
+TEST_F(TestPolynomialFitting, predict_evaluates_fitted_polynomial)
+{
+    math::Matrix<float, 8, 1> x;
+    math::Matrix<float, 8, 1> y;
+
+    for (std::size_t i = 0; i < 8; ++i)
+    {
+        float xi = static_cast<float>(i) * 0.25f;
+        x.at(i, 0) = xi;
+        y.at(i, 0) = 1.0f - 0.5f * xi + 0.25f * xi * xi;
+    }
+
+    fitter.Fit(x, y);
+
+    EXPECT_NEAR(fitter.Predict(0.5f), 1.0f - 0.5f * 0.5f + 0.25f * 0.25f, math::Tolerance<float>());
+    EXPECT_NEAR(fitter.Predict(1.0f), 1.0f - 0.5f * 1.0f + 0.25f * 1.0f,  math::Tolerance<float>());
+}
+
 TEST_F(TestPolynomialFitting, predict_all_coefficients_finite_on_far_from_origin_data)
 {
     math::Matrix<float, 8, 1> x;
diff --git a/numerical/filters/active/AlphaBetaFilter.hpp b/numerical/filters/active/AlphaBetaFilter.hpp
index 752a9bd8..bedcff62 100644
--- a/numerical/filters/active/AlphaBetaFilter.hpp
+++ b/numerical/filters/active/AlphaBetaFilter.hpp
@@ -58,9 +58,7 @@ namespace filters
         : samplePeriod{ Ts }
         , gainAlpha{ alpha }
         , gainBeta{ beta }
-        , gainGamma{ T{} }
         , betaOverTs{ beta / Ts }
-        , twoGammaOverTs2{ T{} }
     {}
 
     template<typename T, std::size_t Order>
diff --git a/numerical/filters/active/KalmanSmoother.hpp b/numerical/filters/active/KalmanSmoother.hpp
index 16fa424a..3d20e1cd 100644
--- a/numerical/filters/active/KalmanSmoother.hpp
+++ b/numerical/filters/active/KalmanSmoother.hpp
@@ -40,11 +40,16 @@ namespace filters
 
         KalmanSmoother() = default;
 
+        struct SmootherParams
+        {
+            const StateMatrix& F;
+            const MeasurementMatrix& H;
+            const StateMatrix& Q;
+            const MeasurementCovariance& R;
+        };
+
         OPTIMIZE_FOR_SPEED SmootherOutput Smooth(
-            const StateMatrix& F,
-            const MeasurementMatrix& H,
-            const StateMatrix& Q,
-            const MeasurementCovariance& R,
+            const SmootherParams& params,
             const std::array<MeasurementVector, MaxSteps>& observations,
             std::size_t numSteps,
             const StateVector& initialState,
@@ -68,10 +73,7 @@ namespace filters
         std::array<KalmanGainMatrix, MaxSteps> kalmanGains_{};
 
         OPTIMIZE_FOR_SPEED float RunForwardPass(
-            const StateMatrix& F,
-            const MeasurementMatrix& H,
-            const StateMatrix& Q,
-            const MeasurementCovariance& R,
+            const SmootherParams& params,
             const std::array<MeasurementVector, MaxSteps>& observations,
             std::size_t numSteps,
             const StateVector& initialState,
@@ -97,10 +99,7 @@ namespace filters
     OPTIMIZE_FOR_SPEED
         typename KalmanSmoother<StateSize, MeasurementSize, MaxSteps>::SmootherOutput
         KalmanSmoother<StateSize, MeasurementSize, MaxSteps>::Smooth(
-            const StateMatrix& F,
-            const MeasurementMatrix& H,
-            const StateMatrix& Q,
-            const MeasurementCovariance& R,
+            const SmootherParams& params,
             const std::array<MeasurementVector, MaxSteps>& observations,
             std::size_t numSteps,
             const StateVector& initialState,
@@ -109,18 +108,15 @@ namespace filters
         really_assert(numSteps >= 2 && numSteps <= MaxSteps);
 
         SmootherOutput output;
-        output.logLikelihood = RunForwardPass(F, H, Q, R, observations, numSteps, initialState, initialCovariance);
-        RunBackwardPass(F, H, numSteps, output);
+        output.logLikelihood = RunForwardPass(params, observations, numSteps, initialState, initialCovariance);
+        RunBackwardPass(params.F, params.H, numSteps, output);
         return output;
     }
 
     template<std::size_t StateSize, std::size_t MeasurementSize, std::size_t MaxSteps>
     OPTIMIZE_FOR_SPEED float
     KalmanSmoother<StateSize, MeasurementSize, MaxSteps>::RunForwardPass(
-        const StateMatrix& F,
-        const MeasurementMatrix& H,
-        const StateMatrix& Q,
-        const MeasurementCovariance& R,
+        const SmootherParams& params,
         const std::array<MeasurementVector, MaxSteps>& observations,
         std::size_t numSteps,
         const StateVector& initialState,
@@ -133,18 +129,18 @@ namespace filters
 
         for (std::size_t t = 0; t < numSteps; ++t)
         {
-            const auto nu = observations[t] - H * predictedMeans_[t];
-            const auto S = math::CongruenceTransform(H, predictedCovariances_[t]) + R;
+            const auto nu = observations[t] - params.H * predictedMeans_[t];
+            const auto S = math::CongruenceTransform(params.H, predictedCovariances_[t]) + params.R;
 
             const auto K = solvers::SolveSystem<float, MeasurementSize, StateSize>(
-                S, H * predictedCovariances_[t])
+                S, params.H * predictedCovariances_[t])
                                .Transpose();
 
             filteredMeans_[t] = predictedMeans_[t] + K * nu;
 
-            const auto IminusKH = StateMatrix::Identity() - K * H;
+            const auto IminusKH = StateMatrix::Identity() - K * params.H;
             filteredCovariances_[t] =
-                math::CongruenceTransform(IminusKH, predictedCovariances_[t]) + math::CongruenceTransform(K, R);
+                math::CongruenceTransform(IminusKH, predictedCovariances_[t]) + math::CongruenceTransform(K, params.R);
 
             kalmanGains_[t] = K;
 
@@ -152,8 +148,8 @@ namespace filters
 
             if (t < numSteps - 1)
             {
-                predictedMeans_[t + 1] = F * filteredMeans_[t];
-                predictedCovariances_[t + 1] = math::CongruenceTransform(F, filteredCovariances_[t]) + Q;
+                predictedMeans_[t + 1] = params.F * filteredMeans_[t];
+                predictedCovariances_[t + 1] = math::CongruenceTransform(params.F, filteredCovariances_[t]) + params.Q;
             }
         }
 
@@ -244,7 +240,7 @@ namespace filters
             const StateVector& initialState,
             const StateMatrix& initialCovariance)
     {
-        return Smooth(plant.A, plant.C, Q, R, observations, numSteps, initialState, initialCovariance);
+        return Smooth({ plant.A, plant.C, Q, R }, observations, numSteps, initialState, initialCovariance);
     }
 
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
diff --git a/numerical/filters/active/test/TestKalmanSmoother.cpp b/numerical/filters/active/test/TestKalmanSmoother.cpp
index de5411b9..63d146a5 100644
--- a/numerical/filters/active/test/TestKalmanSmoother.cpp
+++ b/numerical/filters/active/test/TestKalmanSmoother.cpp
@@ -47,7 +47,7 @@ namespace
 
         Smoother::SmootherOutput RunSmoother()
         {
-            return smoother.Smooth(F, H, Q, R, observations, T, x0, P0);
+            return smoother.Smooth({F, H, Q, R}, observations, T, x0, P0);
         }
     };
 
@@ -80,7 +80,7 @@ namespace
 
         Smoother::SmootherOutput RunSmoother()
         {
-            return smoother.Smooth(F, H, Q, R, obs, T, x0, P0);
+            return smoother.Smooth({F, H, Q, R}, obs, T, x0, P0);
         }
     };
 }
@@ -210,7 +210,7 @@ TEST_F(TestKalmanSmoother, smoothed_rmse_is_less_than_filtered_rmse)
     }
 
     Smoother synthSmoother;
-    const auto output = synthSmoother.Smooth(F, H, Q, R, syntheticObs, T, x0, P0);
+    const auto output = synthSmoother.Smooth({F, H, Q, R}, syntheticObs, T, x0, P0);
 
     float filteredMse = 0.0f;
     float smoothedMse = 0.0f;
@@ -270,7 +270,7 @@ TEST_F(TestKalmanSmoother, smoothed_covariance_at_final_step_equals_filtered_cov
 TEST_F(TestKalmanSmoother, minimum_steps_boundary_produces_finite_output)
 {
     constexpr std::size_t minSteps = 2;
-    const auto output = smoother.Smooth(F, H, Q, R, observations, minSteps, x0, P0);
+    const auto output = smoother.Smooth({F, H, Q, R}, observations, minSteps, x0, P0);
 
     EXPECT_TRUE(std::isfinite(output.logLikelihood));
     for (std::size_t t = 0; t < minSteps; ++t)
@@ -281,7 +281,7 @@ TEST_F(TestKalmanSmoother, minimum_steps_boundary_produces_finite_output)
 TEST_F(TestKalmanSmoother, near_zero_process_noise_produces_no_nan)
 {
     StateMatrix Qsmall{ { 1e-6f, 0.0f }, { 0.0f, 1e-6f } };
-    const auto output = smoother.Smooth(F, H, Qsmall, R, observations, T, x0, P0);
+    const auto output = smoother.Smooth({F, H, Qsmall, R}, observations, T, x0, P0);
 
     EXPECT_TRUE(std::isfinite(output.logLikelihood));
     for (std::size_t t = 0; t < T; ++t)
@@ -321,7 +321,7 @@ TEST_F(TestKalmanSmoother, nees_within_chi_squared_band_on_simulated_system)
             trialObs[t] = MeasurementVector{ { trialTrue[t].at(0, 0) + measNoise(rng) } };
 
         Smoother trialSmoother;
-        const auto out = trialSmoother.Smooth(F, H, Q, R, trialObs, T, x0, P0);
+        const auto out = trialSmoother.Smooth({F, H, Q, R}, trialObs, T, x0, P0);
 
         for (std::size_t t = 0; t < T; ++t)
         {
diff --git a/numerical/filters/passive/MedianFilter.hpp b/numerical/filters/passive/MedianFilter.hpp
index 6372b131..dc7e6ec5 100644
--- a/numerical/filters/passive/MedianFilter.hpp
+++ b/numerical/filters/passive/MedianFilter.hpp
@@ -26,12 +26,11 @@ namespace filters::passive
     private:
         std::array<T, N> window;
         std::array<T, N> scratch;
-        std::size_t head;
+        std::size_t head{ 0 };
     };
 
     template<typename T, std::size_t N>
     MedianFilter<T, N>::MedianFilter(T initial) noexcept
-        : head{ 0 }
     {
         window.fill(initial);
         scratch.fill(T{});
diff --git a/numerical/filters/passive/SavitzkyGolayFilter.hpp b/numerical/filters/passive/SavitzkyGolayFilter.hpp
index 97954fd4..b3bdf459 100644
--- a/numerical/filters/passive/SavitzkyGolayFilter.hpp
+++ b/numerical/filters/passive/SavitzkyGolayFilter.hpp
@@ -22,7 +22,7 @@ namespace filters::passive
             math::Matrix<T, Window, Order + 1> A{};
             for (std::size_t row = 0; row < Window; ++row)
             {
-                T offset = static_cast<T>(static_cast<int>(row) - static_cast<int>(half));
+                auto offset = static_cast<T>(static_cast<int>(row) - static_cast<int>(half));
                 T power = T{ 1 };
                 for (std::size_t col = 0; col < Order + 1; ++col)
                 {
@@ -73,12 +73,11 @@ namespace filters::passive
     private:
         static inline const std::array<T, Window> coeffs{ detail::ComputeKernel<T, Window, Order, Deriv>() };
         std::array<T, Window> line;
-        std::size_t head;
+        std::size_t head{ 0 };
     };
 
     template<typename T, std::size_t Window, std::size_t Order, std::size_t Deriv>
     SavitzkyGolayFilter<T, Window, Order, Deriv>::SavitzkyGolayFilter(T initial) noexcept
-        : head{ 0 }
     {
         line.fill(initial);
     }
diff --git a/numerical/filters/passive/test/TestIir.cpp b/numerical/filters/passive/test/TestIir.cpp
index 26b1bf64..7df2d71a 100644
--- a/numerical/filters/passive/test/TestIir.cpp
+++ b/numerical/filters/passive/test/TestIir.cpp
@@ -8,12 +8,12 @@ namespace
     class TestIir : public ::testing::Test
     {
     protected:
-        static constexpr std::size_t Order = 2;
+        static constexpr std::size_t Order = 3;
 
-        math::RecursiveBuffer<float, Order> MakeCoeffs(float c0, float c1)
+        math::RecursiveBuffer<float, Order> MakeCoeffs(float c0, float c1, float c2 = 0.0f)
         {
             math::RecursiveBuffer<float, Order> buf{};
-            buf = { c0, c1 };
+            buf = { c0, c1, c2 };
             return buf;
         }
     };
diff --git a/numerical/math/ComplexNumber.hpp b/numerical/math/ComplexNumber.hpp
index 51f332d7..2638b537 100644
--- a/numerical/math/ComplexNumber.hpp
+++ b/numerical/math/ComplexNumber.hpp
@@ -123,8 +123,8 @@ namespace math
     }
 
     template<typename QNumberType>
-    std::enable_if_t<std::is_floating_point_v<QNumberType>, QNumberType>
-    Abs(const Complex<QNumberType>& c)
+    requires std::is_floating_point_v<QNumberType>
+    QNumberType Abs(const Complex<QNumberType>& c)
     {
         return math::Hypot(c.Real(), c.Imaginary());
     }
diff --git a/numerical/math/Math.cpp b/numerical/math/Math.cpp
index e92aa64e..09a2d68d 100644
--- a/numerical/math/Math.cpp
+++ b/numerical/math/Math.cpp
@@ -2,80 +2,77 @@
 
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
 
-namespace math
-{
 #ifndef MATH_ABS_OVERRIDE
-    template float Abs<float>(float);
+template float math::Abs<float>(float);
 #endif
 #ifndef MATH_SQRT_OVERRIDE
-    template float Sqrt<float>(float);
+template float math::Sqrt<float>(float);
 #endif
 #ifndef MATH_SIN_OVERRIDE
-    template float Sin<float>(float);
+template float math::Sin<float>(float);
 #endif
 #ifndef MATH_COS_OVERRIDE
-    template float Cos<float>(float);
+template float math::Cos<float>(float);
 #endif
 #ifndef MATH_TAN_OVERRIDE
-    template float Tan<float>(float);
+template float math::Tan<float>(float);
 #endif
 #ifndef MATH_ASIN_OVERRIDE
-    template float Asin<float>(float);
+template float math::Asin<float>(float);
 #endif
 #ifndef MATH_ACOS_OVERRIDE
-    template float Acos<float>(float);
+template float math::Acos<float>(float);
 #endif
 #ifndef MATH_ATAN_OVERRIDE
-    template float Atan<float>(float);
+template float math::Atan<float>(float);
 #endif
 #ifndef MATH_ATAN2_OVERRIDE
-    template float Atan2<float>(float, float);
+template float math::Atan2<float>(float, float);
 #endif
 #ifndef MATH_EXP_OVERRIDE
-    template float Exp<float>(float);
+template float math::Exp<float>(float);
 #endif
 #ifndef MATH_LOG_OVERRIDE
-    template float Log<float>(float);
+template float math::Log<float>(float);
 #endif
 #ifndef MATH_LOG10_OVERRIDE
-    template float Log10<float>(float);
+template float math::Log10<float>(float);
 #endif
 #ifndef MATH_LOG2_OVERRIDE
-    template float Log2<float>(float);
+template float math::Log2<float>(float);
 #endif
 #ifndef MATH_POW_OVERRIDE
-    template float Pow<float>(float, float);
+template float math::Pow<float>(float, float);
 #endif
 #ifndef MATH_SINH_OVERRIDE
-    template float Sinh<float>(float);
+template float math::Sinh<float>(float);
 #endif
 #ifndef MATH_COSH_OVERRIDE
-    template float Cosh<float>(float);
+template float math::Cosh<float>(float);
 #endif
 #ifndef MATH_TANH_OVERRIDE
-    template float Tanh<float>(float);
+template float math::Tanh<float>(float);
 #endif
 #ifndef MATH_HYPOT_OVERRIDE
-    template float Hypot<float>(float, float);
+template float math::Hypot<float>(float, float);
 #endif
 #ifndef MATH_COPYSIGN_OVERRIDE
-    template float Copysign<float>(float, float);
+template float math::Copysign<float>(float, float);
 #endif
 #ifndef MATH_FMOD_OVERRIDE
-    template float Fmod<float>(float, float);
+template float math::Fmod<float>(float, float);
 #endif
 #ifndef MATH_CEIL_OVERRIDE
-    template float Ceil<float>(float);
+template float math::Ceil<float>(float);
 #endif
 #ifndef MATH_FLOOR_OVERRIDE
-    template float Floor<float>(float);
+template float math::Floor<float>(float);
 #endif
 #ifndef MATH_ROUND_OVERRIDE
-    template float Round<float>(float);
+template float math::Round<float>(float);
 #endif
 #ifndef MATH_ERFC_OVERRIDE
-    template float Erfc<float>(float);
+template float math::Erfc<float>(float);
 #endif
-}
 
 #endif
diff --git a/numerical/math/Math.hpp b/numerical/math/Math.hpp
index 227ced56..009ab3bc 100644
--- a/numerical/math/Math.hpp
+++ b/numerical/math/Math.hpp
@@ -224,4 +224,79 @@ namespace math
         return std::erfc(x);
     }
 #endif
+
+#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
+#ifndef MATH_ABS_OVERRIDE
+    extern template float Abs<float>(float);
+#endif
+#ifndef MATH_SQRT_OVERRIDE
+    extern template float Sqrt<float>(float);
+#endif
+#ifndef MATH_SIN_OVERRIDE
+    extern template float Sin<float>(float);
+#endif
+#ifndef MATH_COS_OVERRIDE
+    extern template float Cos<float>(float);
+#endif
+#ifndef MATH_TAN_OVERRIDE
+    extern template float Tan<float>(float);
+#endif
+#ifndef MATH_ASIN_OVERRIDE
+    extern template float Asin<float>(float);
+#endif
+#ifndef MATH_ACOS_OVERRIDE
+    extern template float Acos<float>(float);
+#endif
+#ifndef MATH_ATAN_OVERRIDE
+    extern template float Atan<float>(float);
+#endif
+#ifndef MATH_ATAN2_OVERRIDE
+    extern template float Atan2<float>(float, float);
+#endif
+#ifndef MATH_EXP_OVERRIDE
+    extern template float Exp<float>(float);
+#endif
+#ifndef MATH_LOG_OVERRIDE
+    extern template float Log<float>(float);
+#endif
+#ifndef MATH_LOG10_OVERRIDE
+    extern template float Log10<float>(float);
+#endif
+#ifndef MATH_LOG2_OVERRIDE
+    extern template float Log2<float>(float);
+#endif
+#ifndef MATH_POW_OVERRIDE
+    extern template float Pow<float>(float, float);
+#endif
+#ifndef MATH_SINH_OVERRIDE
+    extern template float Sinh<float>(float);
+#endif
+#ifndef MATH_COSH_OVERRIDE
+    extern template float Cosh<float>(float);
+#endif
+#ifndef MATH_TANH_OVERRIDE
+    extern template float Tanh<float>(float);
+#endif
+#ifndef MATH_HYPOT_OVERRIDE
+    extern template float Hypot<float>(float, float);
+#endif
+#ifndef MATH_COPYSIGN_OVERRIDE
+    extern template float Copysign<float>(float, float);
+#endif
+#ifndef MATH_FMOD_OVERRIDE
+    extern template float Fmod<float>(float, float);
+#endif
+#ifndef MATH_CEIL_OVERRIDE
+    extern template float Ceil<float>(float);
+#endif
+#ifndef MATH_FLOOR_OVERRIDE
+    extern template float Floor<float>(float);
+#endif
+#ifndef MATH_ROUND_OVERRIDE
+    extern template float Round<float>(float);
+#endif
+#ifndef MATH_ERFC_OVERRIDE
+    extern template float Erfc<float>(float);
+#endif
+#endif
 }
diff --git a/numerical/math/MatrixExponential.hpp b/numerical/math/MatrixExponential.hpp
index 41af3e25..778ea85e 100644
--- a/numerical/math/MatrixExponential.hpp
+++ b/numerical/math/MatrixExponential.hpp
@@ -34,7 +34,7 @@ namespace math
 
         void PadeNumeratorDenominator(const SquareMatrix<T, N>& as,
             SquareMatrix<T, N>& num,
-            SquareMatrix<T, N>& den);
+            SquareMatrix<T, N>& den) const;
 
         SquareMatrix<T, N> SolvePade(const SquareMatrix<T, N>& den, const SquareMatrix<T, N>& num);
     };
@@ -42,7 +42,7 @@ namespace math
     template<typename T, std::size_t N>
     void MatrixExponential<T, N>::PadeNumeratorDenominator(const SquareMatrix<T, N>& as,
         SquareMatrix<T, N>& num,
-        SquareMatrix<T, N>& den)
+        SquareMatrix<T, N>& den) const
     {
         const auto identity = SquareMatrix<T, N>::Identity();
         const auto a2 = as * as;
diff --git a/numerical/math/Quaternion.hpp b/numerical/math/Quaternion.hpp
index ab447d67..c9dedda8 100644
--- a/numerical/math/Quaternion.hpp
+++ b/numerical/math/Quaternion.hpp
@@ -27,7 +27,15 @@ namespace math
         static Quaternion FromEulerZYX(T roll, T pitch, T yaw);
         static Quaternion Slerp(const Quaternion& a, const Quaternion& b, T t);
 
-        OPTIMIZE_FOR_SPEED Quaternion operator*(const Quaternion& rhs) const;
+        [[nodiscard]] OPTIMIZE_FOR_SPEED friend Quaternion operator*(const Quaternion& lhs, const Quaternion& rhs)
+        {
+            return Quaternion{
+                lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z,
+                lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y,
+                lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x,
+                lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w
+            };
+        }
         Quaternion operator-() const;
 
         Quaternion Conjugate() const;
@@ -137,17 +145,6 @@ namespace math
         };
     }
 
-    template<typename T>
-    OPTIMIZE_FOR_SPEED Quaternion<T> Quaternion<T>::operator*(const Quaternion& rhs) const
-    {
-        return Quaternion<T>{
-            w * rhs.w - x * rhs.x - y * rhs.y - z * rhs.z,
-            w * rhs.x + x * rhs.w + y * rhs.z - z * rhs.y,
-            w * rhs.y - x * rhs.z + y * rhs.w + z * rhs.x,
-            w * rhs.z + x * rhs.y - y * rhs.x + z * rhs.w
-        };
-    }
-
     template<typename T>
     Quaternion<T> Quaternion<T>::operator-() const
     {
diff --git a/numerical/math/test/CMakeLists.txt b/numerical/math/test/CMakeLists.txt
index 673168f0..5bbcb224 100644
--- a/numerical/math/test/CMakeLists.txt
+++ b/numerical/math/test/CMakeLists.txt
@@ -12,6 +12,7 @@ numerical_add_qemu_test(numerical.math_test)
 
 target_sources(numerical.math_test PRIVATE
     TestComplexNumber.cpp
+    TestMath.cpp
     TestConsistencyMetrics.cpp
     TestCordic.cpp
     TestGivensRotation.cpp
diff --git a/numerical/math/test/TestMath.cpp b/numerical/math/test/TestMath.cpp
new file mode 100644
index 00000000..8f8d34f8
--- /dev/null
+++ b/numerical/math/test/TestMath.cpp
@@ -0,0 +1,158 @@
+#include "numerical/math/Math.hpp"
+#include <cmath>
+#include <gtest/gtest.h>
+#include <numbers>
+
+namespace
+{
+    class TestMath : public ::testing::Test
+    {
+    protected:
+        static constexpr float kEps = 1e-5f;
+        static constexpr float kPi  = std::numbers::pi_v<float>;
+        static constexpr float kPi2 = std::numbers::pi_v<float> / 2.0f;
+        static constexpr float kPi4 = std::numbers::pi_v<float> / 4.0f;
+    };
+}
+
+TEST_F(TestMath, Abs_positive_and_negative)
+{
+    EXPECT_FLOAT_EQ(math::Abs(2.0f), 2.0f);
+    EXPECT_FLOAT_EQ(math::Abs(-3.5f), 3.5f);
+}
+
+TEST_F(TestMath, Sqrt_known_value)
+{
+    EXPECT_NEAR(math::Sqrt(4.0f), 2.0f, kEps);
+}
+
+TEST_F(TestMath, Sin_known_values)
+{
+    EXPECT_NEAR(math::Sin(0.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Sin(kPi2), 1.0f, kEps);
+}
+
+TEST_F(TestMath, Cos_known_values)
+{
+    EXPECT_NEAR(math::Cos(0.0f), 1.0f, kEps);
+    EXPECT_NEAR(math::Cos(kPi), -1.0f, kEps);
+}
+
+TEST_F(TestMath, Tan_known_value)
+{
+    EXPECT_NEAR(math::Tan(0.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Tan(kPi4), 1.0f, 1e-4f);
+}
+
+TEST_F(TestMath, Asin_known_value)
+{
+    EXPECT_NEAR(math::Asin(0.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Asin(1.0f), kPi2, kEps);
+}
+
+TEST_F(TestMath, Acos_known_value)
+{
+    EXPECT_NEAR(math::Acos(1.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Acos(0.0f), kPi2, kEps);
+}
+
+TEST_F(TestMath, Atan_known_value)
+{
+    EXPECT_NEAR(math::Atan(0.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Atan(1.0f), kPi4, kEps);
+}
+
+TEST_F(TestMath, Atan2_known_value)
+{
+    EXPECT_NEAR(math::Atan2(1.0f, 1.0f), kPi4, kEps);
+    EXPECT_NEAR(math::Atan2(0.0f, 1.0f), 0.0f, kEps);
+}
+
+TEST_F(TestMath, Exp_known_value)
+{
+    EXPECT_NEAR(math::Exp(0.0f), 1.0f, kEps);
+    EXPECT_NEAR(math::Exp(1.0f), std::exp(1.0f), kEps);
+}
+
+TEST_F(TestMath, Log_known_value)
+{
+    EXPECT_NEAR(math::Log(1.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Log(std::exp(1.0f)), 1.0f, kEps);
+}
+
+TEST_F(TestMath, Log10_known_value)
+{
+    EXPECT_NEAR(math::Log10(1.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Log10(100.0f), 2.0f, kEps);
+}
+
+TEST_F(TestMath, Log2_known_value)
+{
+    EXPECT_NEAR(math::Log2(1.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Log2(8.0f), 3.0f, kEps);
+}
+
+TEST_F(TestMath, Pow_known_value)
+{
+    EXPECT_NEAR(math::Pow(2.0f, 3.0f), 8.0f, kEps);
+    EXPECT_NEAR(math::Pow(4.0f, 0.5f), 2.0f, kEps);
+}
+
+TEST_F(TestMath, Sinh_known_value)
+{
+    EXPECT_NEAR(math::Sinh(0.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Sinh(1.0f), std::sinh(1.0f), kEps);
+}
+
+TEST_F(TestMath, Cosh_known_value)
+{
+    EXPECT_NEAR(math::Cosh(0.0f), 1.0f, kEps);
+    EXPECT_NEAR(math::Cosh(1.0f), std::cosh(1.0f), kEps);
+}
+
+TEST_F(TestMath, Tanh_known_value)
+{
+    EXPECT_NEAR(math::Tanh(0.0f), 0.0f, kEps);
+    EXPECT_NEAR(math::Tanh(1.0f), std::tanh(1.0f), kEps);
+}
+
+TEST_F(TestMath, Hypot_known_value)
+{
+    EXPECT_NEAR(math::Hypot(3.0f, 4.0f), 5.0f, kEps);
+}
+
+TEST_F(TestMath, Copysign_transfers_sign)
+{
+    EXPECT_FLOAT_EQ(math::Copysign(3.0f, -1.0f), -3.0f);
+    EXPECT_FLOAT_EQ(math::Copysign(-3.0f, 1.0f), 3.0f);
+}
+
+TEST_F(TestMath, Fmod_known_value)
+{
+    EXPECT_NEAR(math::Fmod(5.5f, 2.0f), 1.5f, kEps);
+    EXPECT_NEAR(math::Fmod(7.0f, 3.0f), 1.0f, kEps);
+}
+
+TEST_F(TestMath, Ceil_rounds_up)
+{
+    EXPECT_FLOAT_EQ(math::Ceil(1.2f), 2.0f);
+    EXPECT_FLOAT_EQ(math::Ceil(-1.2f), -1.0f);
+}
+
+TEST_F(TestMath, Floor_rounds_down)
+{
+    EXPECT_FLOAT_EQ(math::Floor(1.9f), 1.0f);
+    EXPECT_FLOAT_EQ(math::Floor(-1.2f), -2.0f);
+}
+
+TEST_F(TestMath, Round_rounds_to_nearest)
+{
+    EXPECT_FLOAT_EQ(math::Round(1.5f), 2.0f);
+    EXPECT_FLOAT_EQ(math::Round(1.4f), 1.0f);
+}
+
+TEST_F(TestMath, Erfc_known_value)
+{
+    EXPECT_NEAR(math::Erfc(0.0f), 1.0f, kEps);
+    EXPECT_NEAR(math::Erfc(1.0f), std::erfc(1.0f), kEps);
+}
diff --git a/numerical/math/test/TestQNumber.cpp b/numerical/math/test/TestQNumber.cpp
index 5ab3cf50..d65f627b 100644
--- a/numerical/math/test/TestQNumber.cpp
+++ b/numerical/math/test/TestQNumber.cpp
@@ -317,3 +317,36 @@ TYPED_TEST(QNumberTest, DivideByZeroDies)
 
     EXPECT_DEATH_IF_SUPPORTED({ TypeParam result = a / zero; }, ""); // NOLINT
 }
+
+TYPED_TEST(QNumberTest, ToFloatFreeFunction_QNumber)
+{
+    TypeParam a(0.25f);
+    EXPECT_NEAR(math::ToFloat(a), 0.25f, math::Tolerance<float>());
+}
+
+TEST(QNumberUtilTest, ToFloatFreeFunction_Float)
+{
+    EXPECT_FLOAT_EQ(math::ToFloat(0.5f), 0.5f);
+    EXPECT_FLOAT_EQ(math::ToFloat(-0.3f), -0.3f);
+}
+
+TEST(QNumberUtilTest, MinMaxLowest_Float)
+{
+    EXPECT_GT(math::Min<float>(), 0.0f);
+    EXPECT_GT(math::Max<float>(), 1.0f);
+    EXPECT_LT(math::Lowest<float>(), 0.0f);
+}
+
+TEST(QNumberUtilTest, MinMaxLowest_Q31)
+{
+    EXPECT_NEAR(math::Min<math::Q31>(), -0.9999f, 1e-3f);
+    EXPECT_NEAR(math::Max<math::Q31>(), 0.9999f, 1e-3f);
+    EXPECT_NEAR(math::Lowest<math::Q31>(), -0.9999f, 1e-3f);
+}
+
+TEST(QNumberUtilTest, MinMaxLowest_Q15)
+{
+    EXPECT_NEAR(math::Min<math::Q15>(), -0.9999f, 1e-3f);
+    EXPECT_NEAR(math::Max<math::Q15>(), 0.9999f, 1e-3f);
+    EXPECT_NEAR(math::Lowest<math::Q15>(), -0.9999f, 1e-3f);
+}
diff --git a/numerical/math/test/TestToeplitz.cpp b/numerical/math/test/TestToeplitz.cpp
index e05daf12..e0b9251a 100644
--- a/numerical/math/test/TestToeplitz.cpp
+++ b/numerical/math/test/TestToeplitz.cpp
@@ -14,7 +14,7 @@ namespace
         : public ::testing::Test
     {
     protected:
-        static constexpr size_t N = 3;
+        static constexpr size_t N = 2;
         using ToeplitzType = math::ToeplitzMatrix<T, N>;
         using VectorType = math::Vector<T, N>;
         using MatrixType = math::Matrix<T, N, N>;
@@ -24,24 +24,19 @@ namespace
             return T(std::max(std::min(f, 0.09f), -0.09f));
         }
 
-        VectorType MakeVector(float a, float b, float c)
+        VectorType MakeVector(float a, float b)
         {
             return VectorType{
                 { MakeValue(a) },
-                { MakeValue(b) },
-                { MakeValue(c) }
+                { MakeValue(b) }
             };
         }
 
-        MatrixType MakeMatrix(
-            float a00, float a01, float a02,
-            float a10, float a11, float a12,
-            float a20, float a21, float a22)
+        MatrixType MakeMatrix(float a00, float a01, float a10, float a11)
         {
             return MatrixType{
-                { MakeValue(a00), MakeValue(a01), MakeValue(a02) },
-                { MakeValue(a10), MakeValue(a11), MakeValue(a12) },
-                { MakeValue(a20), MakeValue(a21), MakeValue(a22) }
+                { MakeValue(a00), MakeValue(a01) },
+                { MakeValue(a10), MakeValue(a11) }
             };
         }
     };
@@ -63,54 +58,47 @@ TYPED_TEST(ToeplitzMatrixTest, DefaultConstructorProducesAllZeroEntries)
 
 TYPED_TEST(ToeplitzMatrixTest, SymmetricConstructorProducesSymmetricToeplitzStructure)
 {
-    auto vec = this->MakeVector(0.06f, 0.03f, 0.01f);
+    auto vec = this->MakeVector(0.06f, 0.03f);
 
     typename TestFixture::ToeplitzType t(vec);
 
     EXPECT_TRUE(t.IsSymmetric());
     EXPECT_TRUE(AreMatricesNear(t.ToFullMatrix(),
         this->MakeMatrix(
-            0.06f, 0.03f, 0.01f,
-            0.03f, 0.06f, 0.03f,
-            0.01f, 0.03f, 0.06f)));
+            0.06f, 0.03f,
+            0.03f, 0.06f)));
 }
 
 TYPED_TEST(ToeplitzMatrixTest, GeneralConstructorProducesAsymmetricToeplitzStructure)
 {
-    auto row = this->MakeVector(0.06f, 0.03f, 0.01f);
-    auto col = this->MakeVector(0.06f, -0.03f, -0.01f);
+    auto row = this->MakeVector(0.06f, 0.03f);
+    auto col = this->MakeVector(0.06f, -0.03f);
 
     typename TestFixture::ToeplitzType t(row, col);
 
     EXPECT_FALSE(t.IsSymmetric());
     EXPECT_TRUE(AreMatricesNear(t.ToFullMatrix(),
         this->MakeMatrix(
-            0.06f,  0.03f,  0.01f,
-           -0.03f,  0.06f,  0.03f,
-           -0.01f, -0.03f,  0.06f)));
+            0.06f,  0.03f,
+           -0.03f,  0.06f)));
 }
 
 TYPED_TEST(ToeplitzMatrixTest, ElementAccessMatchesConstructedRowAndColumn)
 {
-    auto row = this->MakeVector(0.06f, 0.03f, 0.01f);
-    auto col = this->MakeVector(0.06f, -0.03f, -0.01f);
+    auto row = this->MakeVector(0.06f, 0.03f);
+    auto col = this->MakeVector(0.06f, -0.03f);
     typename TestFixture::ToeplitzType t(row, col);
 
     EXPECT_NEAR(math::ToFloat(t.at(0, 0)),  0.06f, math::Tolerance<float>());
     EXPECT_NEAR(math::ToFloat(t.at(0, 1)),  0.03f, math::Tolerance<float>());
-    EXPECT_NEAR(math::ToFloat(t.at(0, 2)),  0.01f, math::Tolerance<float>());
     EXPECT_NEAR(math::ToFloat(t.at(1, 0)), -0.03f, math::Tolerance<float>());
     EXPECT_NEAR(math::ToFloat(t.at(1, 1)),  0.06f, math::Tolerance<float>());
-    EXPECT_NEAR(math::ToFloat(t.at(1, 2)),  0.03f, math::Tolerance<float>());
-    EXPECT_NEAR(math::ToFloat(t.at(2, 0)), -0.01f, math::Tolerance<float>());
-    EXPECT_NEAR(math::ToFloat(t.at(2, 1)), -0.03f, math::Tolerance<float>());
-    EXPECT_NEAR(math::ToFloat(t.at(2, 2)),  0.06f, math::Tolerance<float>());
 }
 
 TYPED_TEST(ToeplitzMatrixTest, ToFullMatrixConsistentWithAtAccessor)
 {
-    auto row = this->MakeVector(0.06f, 0.03f, 0.01f);
-    auto col = this->MakeVector(0.06f, -0.03f, -0.01f);
+    auto row = this->MakeVector(0.06f, 0.03f);
+    auto col = this->MakeVector(0.06f, -0.03f);
     typename TestFixture::ToeplitzType t(row, col);
 
     auto full = t.ToFullMatrix();
@@ -122,9 +110,9 @@ TYPED_TEST(ToeplitzMatrixTest, ToFullMatrixConsistentWithAtAccessor)
 
 TYPED_TEST(ToeplitzMatrixTest, VectorMultiplicationMatchesFullMatrixMultiply)
 {
-    auto vec = this->MakeVector(0.06f, 0.03f, 0.01f);
+    auto vec = this->MakeVector(0.06f, 0.03f);
     typename TestFixture::ToeplitzType t(vec);
-    auto x = this->MakeVector(0.01f, 0.02f, 0.03f);
+    auto x = this->MakeVector(0.01f, 0.02f);
 
     auto result = t * x;
 
@@ -144,7 +132,7 @@ TYPED_TEST(ToeplitzMatrixTest, VectorMultiplicationMatchesFullMatrixMultiply)
 
 TYPED_TEST(ToeplitzMatrixTest, ZeroVectorMultiplicationProducesZeroVector)
 {
-    auto vec = this->MakeVector(0.06f, 0.03f, 0.01f);
+    auto vec = this->MakeVector(0.06f, 0.03f);
     typename TestFixture::ToeplitzType t(vec);
     typename TestFixture::VectorType zero;
 
@@ -156,22 +144,21 @@ TYPED_TEST(ToeplitzMatrixTest, ZeroVectorMultiplicationProducesZeroVector)
 
 TYPED_TEST(ToeplitzMatrixTest, AdditionProducesCorrectToeplitzSum)
 {
-    typename TestFixture::ToeplitzType t1(this->MakeVector(0.04f, 0.02f, 0.01f));
-    typename TestFixture::ToeplitzType t2(this->MakeVector(0.02f, 0.01f, 0.005f));
+    typename TestFixture::ToeplitzType t1(this->MakeVector(0.04f, 0.02f));
+    typename TestFixture::ToeplitzType t2(this->MakeVector(0.02f, 0.01f));
 
     auto result = t1 + t2;
 
     EXPECT_TRUE(AreMatricesNear(result.ToFullMatrix(),
         this->MakeMatrix(
-            0.06f, 0.03f, 0.015f,
-            0.03f, 0.06f, 0.03f,
-            0.015f, 0.03f, 0.06f)));
+            0.06f, 0.03f,
+            0.03f, 0.06f)));
 }
 
 TYPED_TEST(ToeplitzMatrixTest, AdditionIsCommutative)
 {
-    typename TestFixture::ToeplitzType t1(this->MakeVector(0.04f, 0.02f, 0.01f));
-    typename TestFixture::ToeplitzType t2(this->MakeVector(0.02f, 0.01f, 0.005f));
+    typename TestFixture::ToeplitzType t1(this->MakeVector(0.04f, 0.02f));
+    typename TestFixture::ToeplitzType t2(this->MakeVector(0.02f, 0.01f));
 
     auto ab = t1 + t2;
     auto ba = t2 + t1;
@@ -181,24 +168,22 @@ TYPED_TEST(ToeplitzMatrixTest, AdditionIsCommutative)
 
 TYPED_TEST(ToeplitzMatrixTest, SubtractionProducesCorrectToeplitzDifference)
 {
-    typename TestFixture::ToeplitzType t1(this->MakeVector(0.06f, 0.03f, 0.015f));
-    typename TestFixture::ToeplitzType t2(this->MakeVector(0.02f, 0.01f, 0.005f));
+    typename TestFixture::ToeplitzType t1(this->MakeVector(0.06f, 0.03f));
+    typename TestFixture::ToeplitzType t2(this->MakeVector(0.02f, 0.01f));
 
     auto result = t1 - t2;
 
     EXPECT_TRUE(AreMatricesNear(result.ToFullMatrix(),
         this->MakeMatrix(
-            0.04f, 0.02f, 0.01f,
-            0.02f, 0.04f, 0.02f,
-            0.01f, 0.02f, 0.04f)));
+            0.04f, 0.02f,
+            0.02f, 0.04f)));
 }
 
 TYPED_TEST(ToeplitzMatrixTest, IsToeplitzMatrixAcceptsValidPattern)
 {
     auto matrix = this->MakeMatrix(
-        0.06f, 0.03f, 0.01f,
-        0.03f, 0.06f, 0.03f,
-        0.01f, 0.03f, 0.06f);
+        0.06f, 0.03f,
+        0.03f, 0.06f);
 
     EXPECT_TRUE(TestFixture::ToeplitzType::IsToeplitzMatrix(matrix));
 }
@@ -206,9 +191,8 @@ TYPED_TEST(ToeplitzMatrixTest, IsToeplitzMatrixAcceptsValidPattern)
 TYPED_TEST(ToeplitzMatrixTest, IsToeplitzMatrixRejectsNonToeplitzPattern)
 {
     auto matrix = this->MakeMatrix(
-        0.06f, 0.03f, 0.01f,
-        0.03f, 0.06f, 0.03f,
-        0.01f, 0.03f, 0.09f);
+        0.06f, 0.03f,
+        0.03f, 0.09f);
 
     EXPECT_FALSE(TestFixture::ToeplitzType::IsToeplitzMatrix(matrix));
 }
@@ -216,17 +200,16 @@ TYPED_TEST(ToeplitzMatrixTest, IsToeplitzMatrixRejectsNonToeplitzPattern)
 TYPED_TEST(ToeplitzMatrixTest, IsToeplitzMatrixAcceptsZeroMatrix)
 {
     auto matrix = this->MakeMatrix(
-        0.0f, 0.0f, 0.0f,
-        0.0f, 0.0f, 0.0f,
-        0.0f, 0.0f, 0.0f);
+        0.0f, 0.0f,
+        0.0f, 0.0f);
 
     EXPECT_TRUE(TestFixture::ToeplitzType::IsToeplitzMatrix(matrix));
 }
 
 TYPED_TEST(ToeplitzMatrixTest, ToFullMatrixPassesIsToeplitzCheck)
 {
-    auto row = this->MakeVector(0.06f, 0.03f, 0.01f);
-    auto col = this->MakeVector(0.06f, -0.03f, -0.01f);
+    auto row = this->MakeVector(0.06f, 0.03f);
+    auto col = this->MakeVector(0.06f, -0.03f);
     typename TestFixture::ToeplitzType t(row, col);
 
     EXPECT_TRUE(TestFixture::ToeplitzType::IsToeplitzMatrix(t.ToFullMatrix()));
@@ -235,26 +218,22 @@ TYPED_TEST(ToeplitzMatrixTest, ToFullMatrixPassesIsToeplitzCheck)
 TYPED_TEST(ToeplitzMatrixTest, ExtractToeplitzVectorsRecoverRowAndColumn)
 {
     auto matrix = this->MakeMatrix(
-        0.06f, 0.03f, 0.01f,
-        0.03f, 0.06f, 0.03f,
-        0.01f, 0.03f, 0.06f);
+        0.06f, 0.03f,
+        0.03f, 0.06f);
 
     auto [row, col] = TestFixture::ToeplitzType::ExtractToeplitzVectors(matrix);
 
     EXPECT_NEAR(math::ToFloat(row.at(0, 0)), 0.06f, math::Tolerance<float>());
     EXPECT_NEAR(math::ToFloat(row.at(1, 0)), 0.03f, math::Tolerance<float>());
-    EXPECT_NEAR(math::ToFloat(row.at(2, 0)), 0.01f, math::Tolerance<float>());
     EXPECT_NEAR(math::ToFloat(col.at(0, 0)), 0.06f, math::Tolerance<float>());
     EXPECT_NEAR(math::ToFloat(col.at(1, 0)), 0.03f, math::Tolerance<float>());
-    EXPECT_NEAR(math::ToFloat(col.at(2, 0)), 0.01f, math::Tolerance<float>());
 }
 
 TYPED_TEST(ToeplitzMatrixTest, ExtractThenConstructRoundTripMatchesOriginalMatrix)
 {
     auto matrix = this->MakeMatrix(
-        0.06f, 0.03f, 0.01f,
-        -0.03f, 0.06f, 0.03f,
-        -0.01f, -0.03f, 0.06f);
+        0.06f, 0.03f,
+       -0.03f, 0.06f);
 
     auto [row, col] = TestFixture::ToeplitzType::ExtractToeplitzVectors(matrix);
     typename TestFixture::ToeplitzType t(row, col);
@@ -265,9 +244,8 @@ TYPED_TEST(ToeplitzMatrixTest, ExtractThenConstructRoundTripMatchesOriginalMatri
 TYPED_TEST(ToeplitzMatrixTest, ExtractToeplitzVectorsSymmetricMatrixYieldsEqualRowAndColumn)
 {
     auto matrix = this->MakeMatrix(
-        0.06f, 0.03f, 0.01f,
-        0.03f, 0.06f, 0.03f,
-        0.01f, 0.03f, 0.06f);
+        0.06f, 0.03f,
+        0.03f, 0.06f);
 
     auto [row, col] = TestFixture::ToeplitzType::ExtractToeplitzVectors(matrix);
 
@@ -278,9 +256,8 @@ TYPED_TEST(ToeplitzMatrixTest, ExtractToeplitzVectorsSymmetricMatrixYieldsEqualR
 TYPED_TEST(ToeplitzMatrixTest, ExtractToeplitzVectorsZeroMatrixYieldsZeroVectors)
 {
     auto matrix = this->MakeMatrix(
-        0.0f, 0.0f, 0.0f,
-        0.0f, 0.0f, 0.0f,
-        0.0f, 0.0f, 0.0f);
+        0.0f, 0.0f,
+        0.0f, 0.0f);
 
     auto [row, col] = TestFixture::ToeplitzType::ExtractToeplitzVectors(matrix);
 
@@ -293,7 +270,7 @@ TYPED_TEST(ToeplitzMatrixTest, ExtractToeplitzVectorsZeroMatrixYieldsZeroVectors
 
 TYPED_TEST(ToeplitzMatrixTest, CreateToeplitzMatrixFactoryProducesSameResultAsConstructor)
 {
-    auto vec = this->MakeVector(0.06f, 0.03f, 0.01f);
+    auto vec = this->MakeVector(0.06f, 0.03f);
 
     auto fromFactory = math::CreateToeplitzMatrix(vec);
     typename TestFixture::ToeplitzType fromCtor(vec);
@@ -303,13 +280,13 @@ TYPED_TEST(ToeplitzMatrixTest, CreateToeplitzMatrixFactoryProducesSameResultAsCo
 
 TYPED_TEST(ToeplitzMatrixTest, TwoInstancesWithSameInputProduceIdenticalOutput)
 {
-    auto row = this->MakeVector(0.06f, 0.03f, 0.01f);
-    auto col = this->MakeVector(0.06f, -0.03f, -0.01f);
+    auto row = this->MakeVector(0.06f, 0.03f);
+    auto col = this->MakeVector(0.06f, -0.03f);
 
     typename TestFixture::ToeplitzType t1(row, col);
     typename TestFixture::ToeplitzType t2(row, col);
 
-    auto x = this->MakeVector(0.02f, 0.04f, 0.03f);
+    auto x = this->MakeVector(0.02f, 0.04f);
     auto r1 = t1 * x;
     auto r2 = t2 * x;
 
diff --git a/numerical/robust_control/ActiveDisturbanceRejection.hpp b/numerical/robust_control/ActiveDisturbanceRejection.hpp
index 5a892576..e2ec9330 100644
--- a/numerical/robust_control/ActiveDisturbanceRejection.hpp
+++ b/numerical/robust_control/ActiveDisturbanceRejection.hpp
@@ -109,7 +109,7 @@ namespace robust_control
         T woPow{ wo };
         for (std::size_t i = 0; i < n; ++i)
         {
-            const T coeff = static_cast<T>(detail::BinomialCoeff(n, i + 1));
+            const auto coeff = static_cast<T>(detail::BinomialCoeff(n, i + 1));
             gains.at(i, 0) = coeff * woPow;
             woPow *= wo;
         }
@@ -124,7 +124,7 @@ namespace robust_control
         T wcPow{ wc };
         for (std::size_t i = 0; i < Order; ++i)
         {
-            const T coeff = static_cast<T>(detail::BinomialCoeff(Order, i + 1));
+            const auto coeff = static_cast<T>(detail::BinomialCoeff(Order, i + 1));
             gains.at(Order - 1 - i, 0) = coeff * wcPow;
             wcPow *= wc;
         }
diff --git a/numerical/robust_control/test/TestActiveDisturbanceRejection.cpp b/numerical/robust_control/test/TestActiveDisturbanceRejection.cpp
index 6f493d4a..c1088e94 100644
--- a/numerical/robust_control/test/TestActiveDisturbanceRejection.cpp
+++ b/numerical/robust_control/test/TestActiveDisturbanceRejection.cpp
@@ -222,3 +222,24 @@ TEST_F(TestActiveDisturbanceRejection, reset_mid_run_matches_fresh_instance)
     EXPECT_FLOAT_EQ(plant.y, freshPlant.y);
     EXPECT_FLOAT_EQ(adrc.AppliedPrev(), fresh.AppliedPrev());
 }
+
+TEST(TestBinomialCoeff, k_zero_returns_one)
+{
+    EXPECT_EQ(robust_control::detail::BinomialCoeff(5, 0), 1u);
+}
+
+TEST(TestBinomialCoeff, k_equals_n_returns_one)
+{
+    EXPECT_EQ(robust_control::detail::BinomialCoeff(4, 4), 1u);
+}
+
+TEST(TestBinomialCoeff, k_greater_than_n_returns_zero)
+{
+    EXPECT_EQ(robust_control::detail::BinomialCoeff(3, 5), 0u);
+}
+
+TEST(TestBinomialCoeff, known_interior_values)
+{
+    EXPECT_EQ(robust_control::detail::BinomialCoeff(4, 2), 6u);
+    EXPECT_EQ(robust_control::detail::BinomialCoeff(5, 3), 10u);
+}
diff --git a/numerical/solvers/SingularValueDecomposition.hpp b/numerical/solvers/SingularValueDecomposition.hpp
index fcd07cf1..649e0b38 100644
--- a/numerical/solvers/SingularValueDecomposition.hpp
+++ b/numerical/solvers/SingularValueDecomposition.hpp
@@ -40,7 +40,7 @@ namespace solvers
         using Betas = std::array<T, Cols>;
 
         void Bidiagonalize(math::Matrix<T, Rows, Cols>& bidiag,
-            LeftVectors& leftVecs, Betas& leftBeta, RightVectors& rightVecs, Betas& rightBeta);
+            LeftVectors& leftVecs, Betas& leftBeta, RightVectors& rightVecs, Betas& rightBeta) const;
         void ExtractBidiagonal(const math::Matrix<T, Rows, Cols>& bidiag);
         void AccumulateU(const LeftVectors& leftVecs, const Betas& leftBeta);
         void AccumulateV(const RightVectors& rightVecs, const Betas& rightBeta);
@@ -61,7 +61,7 @@ namespace solvers
     template<typename T, std::size_t Rows, std::size_t Cols>
     void SingularValueDecomposition<T, Rows, Cols>::Bidiagonalize(
         math::Matrix<T, Rows, Cols>& bidiag,
-        LeftVectors& leftVecs, Betas& leftBeta, RightVectors& rightVecs, Betas& rightBeta)
+        LeftVectors& leftVecs, Betas& leftBeta, RightVectors& rightVecs, Betas& rightBeta) const
     {
         math::Vector<T, Rows> lvec{};
         math::Vector<T, Cols> rvec{};