algebra_util.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include <array>
#include <cmath>
#include <limits>
#include <optional>
 
#include "akreon/linalg/constants.hpp"
#include "vector_concepts.hpp"
 
#define MAX_STEP_EVAL 256
 
#if defined(__cpp_lib_constexpr_cmath) && __cpp_lib_constexpr_cmath >= 202306L
#define LINALG_HAS_CXX26_CONSTEXPR_CMATH 1
#else
#define LINALG_HAS_CXX26_CONSTEXPR_CMATH 0
#endif
 
#if defined(__cpp_lib_constexpr_complex) &&                                    \
    __cpp_lib_constexpr_complex >= 202306L
#define LINALG_HAS_CXX26_CONSTEXPR_COMPLEX_MATH 1
#else
#define LINALG_HAS_CXX26_CONSTEXPR_COMPLEX_MATH 0
#endif
 
namespace linalg {
 
template <typename T, std::size_t N>
constexpr auto array_to_tuple(const std::array<T, N>& x) {
  auto to_tuple = [&]<std::size_t... Is>(const auto& arr,
                                         std::index_sequence<Is...>) {
    return std::make_tuple(arr[Is]...);
  };
  return to_tuple(x, std::make_index_sequence<N>{});
}
 
template <Field T>
constexpr bool is_real(const std::complex<T>& z, real_type_t<T> tol = 1e-10) {
  return z.imag() <= tol;
}
 
template <Real T>
constexpr std::optional<T> constexpr_div(T num, T den) {
  if (!(den != T(0))) return std::nullopt;
  return num / den;
}
 
template <Real T>
constexpr T constexpr_div_safe(T num, T den, T fallback = T(0)) {
  if (const auto r = constexpr_div(num, den); r.has_value()) return *r;
  return fallback;
}
 
template <Real T>
constexpr T constexpr_pow(T x, std::size_t power) {
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  return std::pow(x, power);
#else
  T res(1);
  for (std::size_t i = 0; i < power; ++i) res *= x;
  return res;
#endif
}
 
template <Real T>
constexpr std::optional<T> constexpr_sqrt(T x, T g = 1e-10, T tol = 1e-10) {
  if (!(x >= T(0))) return std::nullopt;
  if (x == T(0)) return T(0);
 
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  (void)g;
  (void)tol;
  return std::sqrt(x);
#else
  if (!(g > T(0))) g = T(1);
 
  auto abs = [](const auto& val) { return val < T(0) ? -val : val; };
 
  for (std::size_t k = 0; k < MAX_STEP_EVAL; ++k) {
    if (!(g > T(0))) g = T(1);
    if (abs(g * g - x) <= tol) return g;
    g = T(0.5) * (g + constexpr_div_safe(x, g, g));
  }
  return g;
#endif
}
 
template <Real T>
constexpr T constexpr_sqrt_safe(T x, T floor = T(0), T g = 1e-10,
                                T tol = 1e-10) {
  if (const auto r = constexpr_sqrt(x, g, tol); r.has_value()) return *r;
 
  const auto floor_safe = (floor >= T(0)) ? floor : T(0);
  if (const auto rf = constexpr_sqrt(floor_safe, g, tol); rf.has_value())
    return *rf;
 
  return T(0);
}
 
template <Real T>
constexpr T constexpr_abs(T x) {
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  return std::abs(x);
#else
  return x < T(0) ? T(-1) * x : x;
#endif
}
 
template <Complex T>
constexpr auto constexpr_abs(const T& x) {
#if LINALG_HAS_CXX26_CONSTEXPR_COMPLEX_MATH
  return std::abs(x);
#else
  return constexpr_sqrt_safe(x.real() * x.real() + x.imag() * x.imag());
#endif
}
 
template <Real T>
constexpr std::optional<T> constexpr_nrt(T x, std::size_t root, T g = 1e-10,
                                         T tol = 1e-10) {
  if (root == 0) return std::nullopt;
  if (x == T(0)) return T(0);
 
  T sign = T(1);
  if (x < T(0)) {
    if (root % 2 == 0) return std::nullopt;
    sign = T(-1);
    x = -x;
  }
 
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  (void)g;
  (void)tol;
  if (root == 1) return sign * x;
  if (root == 2) return sign * std::sqrt(x);
  if (root == 3) return sign * std::cbrt(x);
  return sign * std::pow(x, T(1) / T(root));
#else
  if (!(g > T(0))) g = T(1);
 
  auto pow_n = [&](T base, std::size_t n) {
    T pow(1);
    for (std::size_t i = 0; i < n; ++i) pow *= base;
    return pow;
  };
 
  for (std::size_t k = 0; k < MAX_STEP_EVAL; ++k) {
    const auto gkm1 = pow_n(g, root - 1);
    if (!(gkm1 > T(0))) return std::nullopt;
 
    g = (T(root - 1) * g + constexpr_div_safe(x, gkm1, g)) / T(root);
    if (constexpr_abs(pow_n(g, root) - x) <= tol) return sign * g;
  }
  return sign * g;
#endif
}
 
template <Real T>
constexpr T constexpr_nrt_safe(T x, std::size_t root, T floor = T(0),
                               T g = 1e-10, T tol = 1e-10) {
  if (const auto r = constexpr_nrt(x, root, g, tol); r.has_value()) return *r;
 
  const auto floor_safe = (floor >= T(0)) ? floor : T(0);
  if (const auto rf = constexpr_nrt(floor_safe, root, g, tol); rf.has_value())
    return *rf;
 
  return T(0);
}
 
template <Real T>
constexpr auto constexpr_sin(T x, int terms = 16) {
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  (void)terms;
  return std::sin(x);
#else
  T sum = x, term = x;
  for (int k = 1; k < terms; ++k) {
    term *= (-(x * x)) / (T(2 * k) * T(2 * k + 1));
    sum += term;
  }
  return sum;
#endif
}
 
template <Real T>
constexpr auto constexpr_cos(T x, int terms = 16) {
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  (void)terms;
  return std::cos(x);
#else
  T sum = T(1), term = T(1);
  for (int k = 1; k < terms; ++k) {
    term *= (-(x * x)) / (T(2 * k - 1) * T(2 * k));
    sum += term;
  }
  return sum;
#endif
}
 
template <Real T>
constexpr std::optional<T> constexpr_atanh(T x, int terms = 16) {
  if (constexpr_abs(x) >= T(1)) return std::nullopt;
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  (void)terms;
  return std::atanh(x);
#else
  const T x2 = x * x;
  T term(x), sum(x);
  for (int k = 1; k < terms; ++k) {
    term *= x2;
    sum += (term / T(2 * k + 1));
  }
  return sum;
#endif
}
 
template <Real T>
constexpr T constexpr_atanh_safe(T x, T fallback = T(0), int terms = 16) {
  if (const auto r = constexpr_atanh(x, terms); r.has_value()) return *r;
  return fallback;
}
 
template <Real T>
constexpr auto constexpr_expi(T theta, int terms = 16) {
#if LINALG_HAS_CXX26_CONSTEXPR_COMPLEX_MATH
  (void)terms;
  return std::exp(std::complex<T>(T(0), theta));
#else
  return std::complex<T>(constexpr_cos(theta, terms),
                         constexpr_sin(theta, terms));
#endif
}
 
template <Real T>
constexpr auto constexpr_exp(T x, int terms = 16) {
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  (void)terms;
  return std::exp(x);
#else
  T sum(1), term(1);
  for (int k = 1; k < terms; ++k) {
    term *= x / k;
    sum += term;
  }
  return sum;
#endif
}
 
template <Complex T>
constexpr auto constexpr_exp_lambda(const T& lambda, int terms = 16) {
#if LINALG_HAS_CXX26_CONSTEXPR_COMPLEX_MATH
  (void)terms;
  return std::exp(lambda);
#else
  return constexpr_exp(lambda.real()) * constexpr_expi(lambda.imag(), terms);
#endif
}
 
template <Real T>
constexpr std::optional<T> constexpr_ln(T x_in, int terms = 16) {
  if (!(x_in > T(0))) return std::nullopt;
 
#if LINALG_HAS_CXX26_CONSTEXPR_CMATH
  (void)terms;
  return std::log(x_in);
#else
  T x = x_in;
  int e = 0;
 
  for (int i = 0; i < 64 && x > T(2); ++i) {
    x *= T(0.5);
    ++e;
  }
  for (int i = 0; i < 64 && x < T(0.5); ++i) {
    x *= T(2);
    --e;
  }
 
  const T y = constexpr_div_safe(x - T(1), x + T(1), T(0));
  const auto a = constexpr_atanh(y, terms);
  if (!a) return std::nullopt;
 
  return T(2) * (*a) + T(e) * ln2_v<T>;
#endif
}
 
template <Real T>
constexpr T constexpr_ln_safe(T x, T floor = std::numeric_limits<T>::min(),
                              int terms = 16) {
  const T floor_safe = (floor > T(0)) ? floor : std::numeric_limits<T>::min();
  const T x_safe = (x > floor_safe) ? x : floor_safe;
  const auto l = constexpr_ln(x_safe, terms);
  if (l.has_value()) return *l;
 
  const auto l_floor = constexpr_ln(floor_safe, terms);
  return l_floor.has_value() ? *l_floor : T(0);
}
 
template <std::size_t Start, std::size_t End, typename Func>
  requires(End - Start >= 0)
constexpr void constexpr_for(Func&& f) {
  if constexpr (Start < End) {
    f(std::integral_constant<std::size_t, Start>{});
    constexpr_for<Start + 1, End>(std::forward<Func>(f));
  }
}
 
template <Complex T>
constexpr bool is_same(const T& a, const T& b, real_type_t<T> tol = 1e-10) {
  return constexpr_abs(a - b) < tol;
}
 
template <VectorContainer Vec>
  requires(Real<vector_element_t<Vec>>)
constexpr auto argsort_desc(const Vec& x) {
  static constexpr std::size_t N = vector_dim_v<Vec>;
 
  std::array<std::size_t, N> p;
  for (std::size_t i = 0; i < N; ++i) p[i] = i;
  for (std::size_t i = 0; i + 1 < N; ++i) {
    std::size_t largest = i;
    for (std::size_t j = i + 1; j < N; ++j)
      if (x(p[j]) > x(p[largest])) largest = j;
    if (largest != i) std::swap(p[i], p[largest]);
  }
  return p;
}
 
} // namespace linalg