feat!: add math custom layer (#238)

* add math custom layer

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* fix: correct std::sqrtf usage and restore missing array include

Co-authored-by: gabrielfrasantos <21131318+gabrielfrasantos@users.noreply.github.com>

* fix: resolve macOS constexpr compile failure in SavitzkyGolayFilter

Co-authored-by: gabrielfrasantos <21131318+gabrielfrasantos@users.noreply.github.com>

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>
Co-authored-by: gabrielfrasantos <21131318+gabrielfrasantos@users.noreply.github.com>

Author: 3 people

.github/instructions/numerical-cpp.instructions.md

@@ -27,6 +27,18 @@ Write generic `template<typename T[, std::size_t N]>` with
 Do not implement `Q15`/`Q31` — the generic `T` keeps that a cheap future add.
 Use `std::numbers::pi_v<float>` — never hardcode `3.14159265f`.
 
+## Math Functions
+
+Use `math::Sin`, `math::Cos`, `math::Abs`, `math::Sqrt`, etc. from
+`numerical/math/Math.hpp` in all production code.
+Never call `std::sin`, `std::cos`, `std::abs`, `std::sqrt`, or any other `<cmath>` function directly in production headers or source files.
+If a required `<cmath>` function is not yet wrapped in `Math.hpp`, add it there first — then use `math::FunctionName` at the call site.
+
+Each function is guarded by `#ifndef MATH_<NAME>_OVERRIDE` (e.g. `MATH_SIN_OVERRIDE`).
+To replace a function with a platform-specific implementation, define the corresponding macro
+(e.g. `add_compile_definitions(MATH_SIN_OVERRIDE)` in CMake) and provide your own
+`template<typename T> constexpr T math::Sin(T x)` in a header that is included before any call site.
+
 ## Embedded Optimizations
 
 Every algorithm header MUST include:

.github/instructions/testing.instructions.md

@@ -44,6 +44,7 @@ TEST_F(TestDare, solves_simple_system)
 - `TEST_F` macros go **outside** the anonymous namespace
 - Include `<gtest/gtest.h>` (not `<gmock/gmock.h>`) unless gmock matchers are needed
 - Use `testing::StrictMock<MockType>` for strict mock expectations
+- **Math in tests**: use `std::sin`, `std::cos`, `std::abs`, etc. directly — never `math::Sin`, `math::Cos`, `math::Abs`, etc. in test code
 - **ONLY `StrictMock`**: Never use `testing::NiceMock<>` or bare mock instantiation — `NiceMock` silences unexpected-call warnings, masking test gaps; `StrictMock` enforces all interactions explicitly
 - Test `float` only (single type) — no multi-type tests
 - **No redundant tests** — implement exactly the spec's enumerated cases; no overlapping/extra cases

AGENTS.md

@@ -54,6 +54,16 @@ resource-constrained embedded systems. Real-time, deterministic, no heap.
 `neural_network`, `optimization`, `regularization`, `solvers`, and new:
 `robust_control`, `nonlinear_control`.
 
+## Math functions
+
+- **Production code** calls `math::Sin`, `math::Cos`, `math::Abs`, etc. from
+  `numerical/math/Math.hpp`. Never call `std::` cmath functions directly in production.
+- If a needed `<cmath>` function is missing from `Math.hpp`, add a `constexpr` wrapper there first,
+  then use `math::FunctionName` at the call site.
+- Each function is individually overridable at compile time via `#ifndef MATH_<NAME>_OVERRIDE`
+  (e.g. `MATH_SIN_OVERRIDE`). Define the macro in CMake and supply your own template definition.
+- **Unit tests** call `std::sin`, `std::cos`, `std::abs`, etc. directly. Never use `math::` in tests.
+
 ## Testing
 
 - GoogleTest. **`TEST_F` on `float`** — no `TYPED_TEST`, no multi-type. **Never plain `TEST()`**.

README.md

@@ -89,6 +89,77 @@ xhost +local:docker
 > **Note:** The Dev Container sets `DISPLAY=host.docker.internal:0.0` to forward GUI windows over TCP. On Linux, if you prefer Unix socket forwarding, you can override `DISPLAY` to `:0` inside the container and add a bind mount for `/tmp/.X11-unix`.
 
 
+## Math Function Overrides
+
+All math functions used by the library (`Sin`, `Cos`, `Abs`, `Sqrt`, `Pow`, …) are centralised in
+`numerical/math/Math.hpp` under the `math::` namespace.
+Each function ships with a default `constexpr` implementation that forwards to the corresponding
+`std::` counterpart, but every one of them can be replaced at compile time with a
+platform-specific implementation — CORDIC, hardware intrinsics, look-up tables, etc.
+
+### How to override a function
+
+**1. Define the override macro in CMake** (suppresses the default definition):
+
+```cmake
+target_compile_definitions(your_target PRIVATE MATH_SIN_OVERRIDE MATH_COS_OVERRIDE)
+```
+
+**2. Provide your own template definition** in a header included before any call site:
+
+```cpp
+// platform/PlatformMath.hpp
+#pragma once
+namespace math
+{
+    template<typename T>
+    constexpr T Sin(T x) { return cordic_sin(static_cast<float>(x)); }
+
+    template<typename T>
+    constexpr T Cos(T x) { return cordic_cos(static_cast<float>(x)); }
+}
+```
+
+**3. Include your header** before `numerical/math/Math.hpp` — or use a CMake forced-include so it
+applies to every translation unit automatically:
+
+```cmake
+target_compile_options(your_target PRIVATE -include platform/PlatformMath.hpp)
+```
+
+### Available override macros
+
+| Macro                    | Function         |
+|--------------------------|------------------|
+| `MATH_ABS_OVERRIDE`      | `math::Abs`      |
+| `MATH_SQRT_OVERRIDE`     | `math::Sqrt`     |
+| `MATH_SIN_OVERRIDE`      | `math::Sin`      |
+| `MATH_COS_OVERRIDE`      | `math::Cos`      |
+| `MATH_TAN_OVERRIDE`      | `math::Tan`      |
+| `MATH_ASIN_OVERRIDE`     | `math::Asin`     |
+| `MATH_ACOS_OVERRIDE`     | `math::Acos`     |
+| `MATH_ATAN_OVERRIDE`     | `math::Atan`     |
+| `MATH_ATAN2_OVERRIDE`    | `math::Atan2`    |
+| `MATH_EXP_OVERRIDE`      | `math::Exp`      |
+| `MATH_LOG_OVERRIDE`      | `math::Log`      |
+| `MATH_LOG10_OVERRIDE`    | `math::Log10`    |
+| `MATH_LOG2_OVERRIDE`     | `math::Log2`     |
+| `MATH_POW_OVERRIDE`      | `math::Pow`      |
+| `MATH_SINH_OVERRIDE`     | `math::Sinh`     |
+| `MATH_COSH_OVERRIDE`     | `math::Cosh`     |
+| `MATH_TANH_OVERRIDE`     | `math::Tanh`     |
+| `MATH_HYPOT_OVERRIDE`    | `math::Hypot`    |
+| `MATH_COPYSIGN_OVERRIDE` | `math::Copysign` |
+| `MATH_FMOD_OVERRIDE`     | `math::Fmod`     |
+| `MATH_CEIL_OVERRIDE`     | `math::Ceil`     |
+| `MATH_FLOOR_OVERRIDE`    | `math::Floor`    |
+| `MATH_ROUND_OVERRIDE`    | `math::Round`    |
+| `MATH_ERFC_OVERRIDE`     | `math::Erfc`     |
+
+Overrides are fully compile-time: the default body is excluded from the translation unit, so the
+replacement is used even for inlined calls. Non-overridden functions continue to use the `std::`
+defaults unchanged.
+
 ## Roadmap
 
 Planned algorithms and components are tracked in [ROADMAP.md](ROADMAP.md) — a prioritized backlog

numerical/analysis/Decibels.hpp

@@ -5,8 +5,8 @@
 #endif
 
 #include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/Math.hpp"
 #include <algorithm>
-#include <cmath>
 #include <type_traits>
 
 namespace analysis
@@ -24,14 +24,14 @@ namespace analysis
         static_assert(std::is_floating_point_v<T>, "ToDecibels supports floating-point types only");
         if (ratio <= T{ 0 })
             return DecibelFloor<T>::value;
-        return std::max(T{ 20 } * std::log10(ratio), DecibelFloor<T>::value);
+        return std::max(T{ 20 } * math::Log10(ratio), DecibelFloor<T>::value);
     }
 
     template<typename T>
     OPTIMIZE_FOR_SPEED T FromDecibels(T db)
     {
         static_assert(std::is_floating_point_v<T>, "FromDecibels supports floating-point types only");
-        return std::pow(T{ 10 }, db / T{ 20 });
+        return math::Pow(T{ 10 }, db / T{ 20 });
     }
 
     template<typename T>

numerical/analysis/DiscreteCosineTransform.hpp

@@ -8,7 +8,7 @@
 #include "numerical/analysis/FastFourierTransform.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/ComplexNumber.hpp"
-#include <cmath>
+#include "numerical/math/Math.hpp"
 #include <numbers>
 
 namespace analysis
@@ -59,17 +59,17 @@ namespace analysis
 
         auto& fftResult = fft.Forward(reordered);
 
-        output[0] = QNumberType(math::ToFloat(fftResult[0].Real()) / std::sqrt(static_cast<float>(Length)));
+        output[0] = QNumberType(math::ToFloat(fftResult[0].Real()) / math::Sqrt(static_cast<float>(Length)));
 
         for (std::size_t k = 1; k < Length; ++k)
         {
             float angle = -static_cast<float>(k) * std::numbers::pi_v<float> / (2.0f * Length);
-            float scale = 2.0f / std::sqrt(static_cast<float>(Length));
+            float scale = 2.0f / math::Sqrt(static_cast<float>(Length));
 
             float real = math::ToFloat(fftResult[k].Real());
             float imag = math::ToFloat(fftResult[k].Imaginary());
 
-            output[k] = QNumberType((real * std::cos(angle) - imag * std::sin(angle)) * scale);
+            output[k] = QNumberType((real * math::Cos(angle) - imag * math::Sin(angle)) * scale);
         }
 
         return output;
@@ -78,7 +78,7 @@ namespace analysis
     template<typename QNumberType, std::size_t Length>
     typename DiscreteConsineTransform<QNumberType, Length>::VectorReal& DiscreteConsineTransform<QNumberType, Length>::Inverse(VectorReal& input)
     {
-        float sqrtN = std::sqrt(static_cast<float>(Length));
+        float sqrtN = math::Sqrt(static_cast<float>(Length));
 
         complexBuffer[0] = math::Complex<QNumberType>{ QNumberType(math::ToFloat(input[0]) * sqrtN), QNumberType(0.0f) };
 
@@ -88,8 +88,8 @@ namespace analysis
             float imag = -math::ToFloat(input[Length - k]) * sqrtN / 2.0f;
 
             float angle = static_cast<float>(k) * std::numbers::pi_v<float> / (2.0f * static_cast<float>(Length));
-            float cosine = std::cos(angle);
-            float sine = std::sin(angle);
+            float cosine = math::Cos(angle);
+            float sine = math::Sin(angle);
 
             complexBuffer[k] = math::Complex<QNumberType>{ QNumberType(real * cosine - imag * sine), QNumberType(real * sine + imag * cosine) };
         }

numerical/analysis/DiscreteWaveletTransform.hpp

@@ -6,8 +6,8 @@
 
 #include "infra/util/BoundedVector.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/Math.hpp"
 #include <array>
-#include <cmath>
 #include <cstddef>
 #include <type_traits>
 

numerical/analysis/GoertzelAlgorithm.hpp

@@ -6,7 +6,7 @@
 
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/ComplexNumber.hpp"
-#include <cmath>
+#include "numerical/math/Math.hpp"
 #include <cstddef>
 #include <numbers>
 #include <type_traits>
@@ -42,9 +42,9 @@ namespace analysis
 
     template<typename T>
     GoertzelAlgorithm<T>::GoertzelAlgorithm(std::size_t k, std::size_t blockLength)
-        : coeff{ T{ 2 } * std::cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
-        , cosine{ std::cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
-        , sine{ std::sin(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
+        : coeff{ T{ 2 } * math::Cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
+        , cosine{ math::Cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
+        , sine{ math::Sin(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
         , blockSize{ blockLength }
     {}
 
@@ -82,7 +82,7 @@ namespace analysis
     template<typename T>
     T GoertzelAlgorithm<T>::Magnitude() const
     {
-        return std::sqrt(s1 * s1 + s2 * s2 - coeff * s1 * s2);
+        return math::Sqrt(s1 * s1 + s2 * s2 - coeff * s1 * s2);
     }
 
     template<typename T>
@@ -96,7 +96,7 @@ namespace analysis
     template<typename T>
     T GoertzelAlgorithm<T>::Coefficient(std::size_t k, std::size_t blockLength)
     {
-        return T{ 2 } * std::cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength));
+        return T{ 2 } * math::Cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength));
     }
 }
 

numerical/analysis/HilbertTransform.hpp

@@ -8,8 +8,8 @@
 #include "numerical/analysis/FastFourierTransform.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/ComplexNumber.hpp"
+#include "numerical/math/Math.hpp"
 #include <array>
-#include <cmath>
 #include <numbers>
 #include <type_traits>
 
@@ -101,13 +101,13 @@ namespace analysis
     template<typename T, std::size_t N>
     T AnalyticSignalFft<T, N>::InstantaneousAmplitude(Complex a)
     {
-        return std::sqrt(a.Real() * a.Real() + a.Imaginary() * a.Imaginary());
+        return math::Sqrt(a.Real() * a.Real() + a.Imaginary() * a.Imaginary());
     }
 
     template<typename T, std::size_t N>
     T AnalyticSignalFft<T, N>::InstantaneousPhase(Complex a)
     {
-        return std::atan2(a.Imaginary(), a.Real());
+        return math::Atan2(a.Imaginary(), a.Real());
     }
 
     template<typename T, std::size_t N>
@@ -133,7 +133,7 @@ namespace analysis
             int k{ static_cast<int>(i) - static_cast<int>(centerTap) };
             if (k != 0 && (k % 2) != 0)
             {
-                T w{ T(0.54) - T(0.46) * std::cos(T(2) * std::numbers::pi_v<T> * static_cast<T>(i) / static_cast<T>(Taps - 1)) };
+                T w{ T(0.54) - T(0.46) * math::Cos(T(2) * std::numbers::pi_v<T> * static_cast<T>(i) / static_cast<T>(Taps - 1)) };
                 coeff[i] = (T(2) / (std::numbers::pi_v<T> * static_cast<T>(k))) * w;
             }
         }

numerical/analysis/SignalDetectors.hpp

@@ -5,7 +5,7 @@
 #endif
 
 #include "numerical/math/CompilerOptimizations.hpp"
-#include <cmath>
+#include "numerical/math/Math.hpp"
 #include <cstdint>
 #include <type_traits>
 
@@ -71,7 +71,7 @@ namespace analysis
     template<typename T>
     OPTIMIZE_FOR_SPEED T PeakHold<T>::Update(T x)
     {
-        T m{ std::abs(x) };
+        T m{ math::Abs(x) };
         peak = (m > peak * decay) ? m : peak * decay;
         return peak;
     }
@@ -90,7 +90,7 @@ namespace analysis
     template<typename T>
     OPTIMIZE_FOR_SPEED bool ZeroCrossingCounter<T>::Update(T x)
     {
-        bool crossed{ (previous < T{ 0 } ? x > T{ 0 } : x < T{ 0 }) && (std::abs(x) > hysteresis) };
+        bool crossed{ (previous < T{ 0 } ? x > T{ 0 } : x < T{ 0 }) && (math::Abs(x) > hysteresis) };
         if (crossed)
             ++count;
         previous = x;
@@ -119,7 +119,7 @@ namespace analysis
     OPTIMIZE_FOR_SPEED T RmsEnvelope<T>::Update(T x)
     {
         meanSquare += alpha * (x * x - meanSquare);
-        return std::sqrt(meanSquare);
+        return math::Sqrt(meanSquare);
     }
 
     template<typename T>