vec.hpp

 
#ifndef UFO_NUMERIC_VEC_HPP
#define UFO_NUMERIC_VEC_HPP
 
// UFO
#include <ufo/numeric/math.hpp>
 
// STL
#include <algorithm>
#include <array>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <format>
#include <functional>
#include <limits>
#include <ostream>
#include <ranges>
#include <type_traits>
#include <utility>
 
namespace ufo
{
 
template <std::size_t Dim, class T>
  requires(Dim <= 4 && std::is_arithmetic_v<T>)
struct Vec {
    using value_type = T;
    using size_type  = std::size_t;
 
    std::array<T, Dim> fields{};
 
    constexpr Vec() noexcept = default;
 
    template <std::convertible_to<T>... Args>
      requires(sizeof...(Args) <= Dim)
    explicit(sizeof...(Args) == 1) constexpr Vec(Args... args) noexcept
        : fields{static_cast<T>(args)...}
    {
    }
 
    template <class U>
    constexpr explicit(!std::is_same_v<T, U>) Vec(Vec<Dim, U> const& other) noexcept
    {
        std::ranges::transform(other, begin(), [](U const& v) { return static_cast<T>(v); });
    }
 
    /**************************************************************************************
    |                                                                                     |
    |                                     Accessors                                      |
    |                                                                                     |
    **************************************************************************************/
 
    constexpr auto& operator[](this auto& self, size_type i) noexcept
    {
        return self.fields[i];
    }
 
    [[nodiscard]] constexpr auto& x(this auto& self) noexcept
      requires(Dim >= 1)
    {
        return self[0];
    }
 
    [[nodiscard]] constexpr auto& y(this auto& self) noexcept
      requires(Dim >= 2)
    {
        return self[1];
    }
 
    [[nodiscard]] constexpr auto& z(this auto& self) noexcept
      requires(Dim >= 3)
    {
        return self[2];
    }
 
    [[nodiscard]] constexpr auto& w(this auto& self) noexcept
      requires(Dim >= 4)
    {
        return self[3];
    }
 
    [[nodiscard]] constexpr auto begin(this auto& self) noexcept
    {
        return self.fields.begin();
    }
 
    [[nodiscard]] constexpr auto end(this auto& self) noexcept { return self.fields.end(); }
 
    [[nodiscard]] constexpr auto data(this auto& self) noexcept
    {
        return self.fields.data();
    }
 
    [[nodiscard]] static constexpr std::size_t size() noexcept { return Dim; }
 
    constexpr bool operator==(Vec const&) const = default;
 
    /**************************************************************************************
    |                                                                                     |
    |                              Unary arithmetic operator                              |
    |                                                                                     |
    **************************************************************************************/
 
    constexpr Vec operator+() const noexcept { return *this; }
 
    constexpr Vec operator-() const noexcept
    {
        Vec result;
        std::ranges::transform(*this, result.begin(), std::negate{});
        return result;
    }
 
    constexpr Vec operator~() const noexcept
      requires std::is_integral_v<T>
    {
        Vec result;
        std::ranges::transform(*this, result.begin(), std::bit_not{});
        return result;
    }
 
    /**************************************************************************************
    |                                                                                     |
    |                            Compound assignment operator                             |
    |                                                                                     |
    **************************************************************************************/
 
    constexpr Vec& operator+=(Vec const& rhs) noexcept
    {
        std::ranges::transform(*this, rhs, begin(), std::plus{});
        return *this;
    }
 
    constexpr Vec& operator-=(Vec const& rhs) noexcept
    {
        std::ranges::transform(*this, rhs, begin(), std::minus{});
        return *this;
    }
 
    constexpr Vec& operator*=(Vec const& rhs) noexcept
    {
        std::ranges::transform(*this, rhs, begin(), std::multiplies{});
        return *this;
    }
 
    constexpr Vec& operator/=(Vec const& rhs) noexcept
    {
        std::ranges::transform(*this, rhs, begin(), std::divides{});
        return *this;
    }
 
    constexpr Vec& operator%=(Vec const& rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, rhs, begin(), std::modulus{});
        return *this;
    }
 
    constexpr Vec& operator&=(Vec const& rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, rhs, begin(), std::bit_and{});
        return *this;
    }
 
    constexpr Vec& operator|=(Vec const& rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, rhs, begin(), std::bit_or{});
        return *this;
    }
 
    constexpr Vec& operator^=(Vec const& rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, rhs, begin(), std::bit_xor{});
        return *this;
    }
 
    constexpr Vec& operator<<=(Vec const& rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, rhs, begin(), [](T a, T b) { return a << b; });
        return *this;
    }
 
    constexpr Vec& operator>>=(Vec const& rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, rhs, begin(), [](T a, T b) { return a >> b; });
        return *this;
    }
 
    /**************************************************************************************
    |                                                                                     |
    |                             Scalar compound assignment operator                     |
    |                                                                                     |
    **************************************************************************************/
 
    constexpr Vec& operator+=(T rhs) noexcept
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v + rhs; });
        return *this;
    }
 
    constexpr Vec& operator-=(T rhs) noexcept
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v - rhs; });
        return *this;
    }
 
    constexpr Vec& operator*=(T rhs) noexcept
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v * rhs; });
        return *this;
    }
 
    constexpr Vec& operator/=(T rhs) noexcept
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v / rhs; });
        return *this;
    }
 
    constexpr Vec& operator%=(T rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v % rhs; });
        return *this;
    }
 
    constexpr Vec& operator&=(T rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v & rhs; });
        return *this;
    }
 
    constexpr Vec& operator|=(T rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v | rhs; });
        return *this;
    }
 
    constexpr Vec& operator^=(T rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v ^ rhs; });
        return *this;
    }
 
    constexpr Vec& operator<<=(T rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v <<= rhs; });
        return *this;
    }
 
    constexpr Vec& operator>>=(T rhs) noexcept
      requires std::is_integral_v<T>
    {
        std::ranges::transform(*this, begin(), [rhs](T v) { return v >>= rhs; });
        return *this;
    }
};
 
/**************************************************************************************
|                                                                                     |
|                                   Type Aliases                                      |
|                                                                                     |
**************************************************************************************/
 
using Vec1b = Vec<1, bool>;
using Vec2b = Vec<2, bool>;
using Vec3b = Vec<3, bool>;
using Vec4b = Vec<4, bool>;
using Vec1i = Vec<1, int>;
using Vec2i = Vec<2, int>;
using Vec3i = Vec<3, int>;
using Vec4i = Vec<4, int>;
using Vec1u = Vec<1, unsigned>;
using Vec2u = Vec<2, unsigned>;
using Vec3u = Vec<3, unsigned>;
using Vec4u = Vec<4, unsigned>;
using Vec1f = Vec<1, float>;
using Vec2f = Vec<2, float>;
using Vec3f = Vec<3, float>;
using Vec4f = Vec<4, float>;
using Vec1d = Vec<1, double>;
using Vec2d = Vec<2, double>;
using Vec3d = Vec<3, double>;
using Vec4d = Vec<4, double>;
 
/**************************************************************************************
|                                                                                     |
|                                  Binary operators                                   |
|                                                                                     |
**************************************************************************************/
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator+(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs += rhs;
    return lhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator-(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs -= rhs;
    return lhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator*(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs *= rhs;
    return lhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator/(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs /= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator%(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs %= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator&(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs &= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator|(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs |= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator^(Vec<Dim, T>        lhs,
                                              Vec<Dim, T> const& rhs) noexcept
{
    lhs ^= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator<<(Vec<Dim, T>        lhs,
                                               Vec<Dim, T> const& rhs) noexcept
{
    lhs <<= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator>>(Vec<Dim, T>        lhs,
                                               Vec<Dim, T> const& rhs) noexcept
{
    lhs >>= rhs;
    return lhs;
}
 
template <std::size_t Dim>
[[nodiscard]] constexpr Vec<Dim, bool> operator&&(Vec<Dim, bool>        lhs,
                                                  Vec<Dim, bool> const& rhs) noexcept
{
    std::ranges::transform(lhs, rhs, lhs.begin(), std::logical_and{});
    return lhs;
}
 
template <std::size_t Dim>
[[nodiscard]] constexpr Vec<Dim, bool> operator||(Vec<Dim, bool>        lhs,
                                                  Vec<Dim, bool> const& rhs) noexcept
{
    std::ranges::transform(lhs, rhs, lhs.begin(), std::logical_or{});
    return lhs;
}
 
/**************************************************************************************
|                                                                                     |
|                              Scalar binary operators                                |
|                                                                                     |
**************************************************************************************/
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator+(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs += rhs;
    return lhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator-(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs -= rhs;
    return lhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator*(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs *= rhs;
    return lhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator/(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs /= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator%(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs %= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator&(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs &= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator|(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs |= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator^(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs ^= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator<<(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs <<= rhs;
    return lhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator>>(Vec<Dim, T> lhs, T rhs) noexcept
{
    lhs >>= rhs;
    return lhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator+(T lhs, Vec<Dim, T> rhs) noexcept
{
    rhs += lhs;
    return rhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator-(T lhs, Vec<Dim, T> rhs) noexcept
{
    std::ranges::transform(rhs, rhs.begin(), [lhs](T v) { return lhs - v; });
    return rhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator*(T lhs, Vec<Dim, T> rhs) noexcept
{
    rhs *= lhs;
    return rhs;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> operator/(T lhs, Vec<Dim, T> rhs) noexcept
{
    std::ranges::transform(rhs, rhs.begin(), [lhs](T v) { return lhs / v; });
    return rhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator%(T lhs, Vec<Dim, T> rhs) noexcept
{
    std::ranges::transform(rhs, rhs.begin(), [lhs](T v) { return lhs % v; });
    return rhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator&(T lhs, Vec<Dim, T> rhs) noexcept
{
    rhs &= lhs;
    return rhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator|(T lhs, Vec<Dim, T> rhs) noexcept
{
    rhs |= lhs;
    return rhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator^(T lhs, Vec<Dim, T> rhs) noexcept
{
    rhs ^= lhs;
    return rhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator<<(T lhs, Vec<Dim, T> rhs) noexcept
{
    std::ranges::transform(rhs, rhs.begin(), [lhs](T v) { return lhs << v; });
    return rhs;
}
 
template <std::size_t Dim, std::integral T>
[[nodiscard]] constexpr Vec<Dim, T> operator>>(T lhs, Vec<Dim, T> rhs) noexcept
{
    std::ranges::transform(rhs, rhs.begin(), [lhs](T v) { return lhs >> v; });
    return rhs;
}
 
/**************************************************************************************
|                                                                                     |
|                                      Functions                                      |
|                                                                                     |
**************************************************************************************/
 
template <class T, std::size_t Dim, class U>
[[nodiscard]] constexpr Vec<Dim, T> cast(Vec<Dim, U> const& v) noexcept
{
    return Vec<Dim, T>(v);
}
 
template <class To, std::size_t Dim, class U>
[[nodiscard]] constexpr To convert(Vec<Dim, U> const& v) noexcept
{
    using T = typename To::value_type;
    To                res{};
    std::size_t const D = std::min(res.size(), v.size());
    for (std::size_t i{}; D > i; ++i) {
        res[i] = static_cast<T>(v[i]);
    }
    return res;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr T dot(Vec<Dim, T> const& a, Vec<Dim, T> const& b) noexcept
{
#if defined(__cpp_lib_ranges_zip) && defined(__cpp_lib_fold_left)
    return std::ranges::fold_left(std::views::zip_transform(std::multiplies{}, a, b), T{},
                                  std::plus{});
#else
    T res{};
    for (std::size_t i{}; Dim > i; ++i) {
        res += a[i] * b[i];
    }
    return res;
#endif
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr T length(Vec<Dim, T> const& v) noexcept
{
    return std::sqrt(dot(v, v));
}
 
template <std::size_t Dim, class T>
  requires std::floating_point<T>
[[nodiscard]] constexpr T distance(Vec<Dim, T> const& a, Vec<Dim, T> const& b) noexcept
{
    return length(a - b);
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr T distanceSquared(Vec<Dim, T> const& a,
                                          Vec<Dim, T> const& b) noexcept
{
    return dot(a - b, a - b);
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> normalize(Vec<Dim, T> const& v) noexcept
{
    return v / length(v);
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr bool isNormalized(Vec<Dim, T> const& v) noexcept
{
    if constexpr (std::floating_point<T>) {
        return std::abs(dot(v, v) - T(1)) <= std::numeric_limits<T>::epsilon() * T(8);
    } else {
        return dot(v, v) == T(1);
    }
}
 
template <class T>
[[nodiscard]] constexpr Vec<3, T> cross(Vec<3, T> const& a, Vec<3, T> const& b) noexcept
{
    return Vec<3, T>{(a[1] * b[2]) - (a[2] * b[1]), (a[2] * b[0]) - (a[0] * b[2]),
                     (a[0] * b[1]) - (a[1] * b[0])};
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> reflect(Vec<Dim, T> const& v,
                                            Vec<Dim, T> const& n) noexcept
{
    return v - T(2) * dot(v, n) * n;
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> refract(Vec<Dim, T> const& v, Vec<Dim, T> const& n,
                                            T eta) noexcept
{
    assert(isNormalized(v) && isNormalized(n) && "Input vectors should be normalized");
    Vec<Dim, T> tmp = dot(v, n);
    Vec<Dim, T> k   = T(1) - eta * eta * (T(1) - tmp * tmp);
    return k < T(0) ? Vec<Dim, T>{} : eta * v - (eta * tmp + sqrt(k)) * n;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> min(Vec<Dim, T> v1, Vec<Dim, T> const& v2) noexcept
{
    std::ranges::transform(v1, v2, v1.begin(), [](T a, T b) { return std::min(a, b); });
    return v1;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr T min(Vec<Dim, T> const& v) noexcept
{
    return *std::ranges::min_element(v);
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> max(Vec<Dim, T> v1, Vec<Dim, T> const& v2) noexcept
{
    std::ranges::transform(v1, v2, v1.begin(), [](T a, T b) { return std::max(a, b); });
    return v1;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr T max(Vec<Dim, T> const& v) noexcept
{
    return *std::ranges::max_element(v);
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr std::size_t minIndex(Vec<Dim, T> const& v) noexcept
{
    return static_cast<std::size_t>(std::ranges::min_element(v) - v.begin());
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr std::size_t maxIndex(Vec<Dim, T> const& v) noexcept
{
    return static_cast<std::size_t>(std::ranges::max_element(v) - v.begin());
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr T sum(Vec<Dim, T> const& v) noexcept
{
#ifdef __cpp_lib_fold_left
    return std::ranges::fold_left(v, T{}, std::plus{});
#else
    T res{};
    for (auto x : v) {
        res += x;
    }
    return res;
#endif
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr T product(Vec<Dim, T> const& v) noexcept
{
#ifdef __cpp_lib_fold_left
    return std::ranges::fold_left(v, T{1}, std::multiplies{});
#else
    T res{1};
    for (auto x : v) {
        res *= x;
    }
    return res;
#endif
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> abs(Vec<Dim, T> v) noexcept
{
    std::ranges::transform(v, v.begin(), [](T x) { return std::abs(x); });
    return v;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> clamp(Vec<Dim, T> v, Vec<Dim, T> const& lo,
                                          Vec<Dim, T> const& hi) noexcept
{
#ifdef __cpp_lib_ranges_zip
    for (auto [vi, li, hi_i] : std::views::zip(v, lo, hi)) {
        vi = std::clamp(vi, li, hi_i);
    }
#else
    for (std::size_t i{}; Dim > i; ++i) {
        v[i] = std::clamp(v[i], lo[i], hi[i]);
    }
#endif
    return v;
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> ceil(Vec<Dim, T> v) noexcept
{
    std::ranges::transform(v, v.begin(), [](T x) { return std::ceil(x); });
    return v;
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> floor(Vec<Dim, T> v) noexcept
{
    std::ranges::transform(v, v.begin(), [](T x) { return std::floor(x); });
    return v;
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> trunc(Vec<Dim, T> v) noexcept
{
    std::ranges::transform(v, v.begin(), [](T x) { return std::trunc(x); });
    return v;
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> round(Vec<Dim, T> v) noexcept
{
    std::ranges::transform(v, v.begin(), [](T x) { return std::round(x); });
    return v;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, bool> equal(Vec<Dim, T> const& v1,
                                             Vec<Dim, T> const& v2) noexcept
{
    Vec<Dim, bool> res;
    std::ranges::transform(v1, v2, res.begin(), std::equal_to{});
    return res;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, bool> notEqual(Vec<Dim, T> const& v1,
                                                Vec<Dim, T> const& v2) noexcept
{
    Vec<Dim, bool> res;
    std::ranges::transform(v1, v2, res.begin(), std::not_equal_to{});
    return res;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, bool> less(Vec<Dim, T> const& v1,
                                            Vec<Dim, T> const& v2) noexcept
{
    Vec<Dim, bool> res;
    std::ranges::transform(v1, v2, res.begin(), std::less{});
    return res;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, bool> lessEqual(Vec<Dim, T> const& v1,
                                                 Vec<Dim, T> const& v2) noexcept
{
    Vec<Dim, bool> res;
    std::ranges::transform(v1, v2, res.begin(), std::less_equal{});
    return res;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, bool> greater(Vec<Dim, T> const& v1,
                                               Vec<Dim, T> const& v2) noexcept
{
    Vec<Dim, bool> res;
    std::ranges::transform(v1, v2, res.begin(), std::greater{});
    return res;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, bool> greaterEqual(Vec<Dim, T> const& v1,
                                                    Vec<Dim, T> const& v2) noexcept
{
    Vec<Dim, bool> res;
    std::ranges::transform(v1, v2, res.begin(), std::greater_equal{});
    return res;
}
 
template <std::size_t Dim>
[[nodiscard]] constexpr bool all(Vec<Dim, bool> const& v) noexcept
{
    return std::ranges::all_of(v, std::identity{});
}
 
template <std::size_t Dim>
[[nodiscard]] constexpr bool any(Vec<Dim, bool> const& v) noexcept
{
    return std::ranges::any_of(v, std::identity{});
}
 
template <std::size_t Dim>
[[nodiscard]] constexpr bool some(Vec<Dim, bool> const& v) noexcept
{
    return any(v) && !all(v);
}
 
template <std::size_t Dim>
[[nodiscard]] constexpr bool none(Vec<Dim, bool> const& v) noexcept
{
    return std::ranges::none_of(v, std::identity{});
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr bool isnan(Vec<Dim, T> const& v) noexcept
{
    return std::ranges::any_of(v, [](T x) { return std::isnan(x); });
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr bool isfinite(Vec<Dim, T> const& v) noexcept
{
    return std::ranges::all_of(v, [](T x) { return std::isfinite(x); });
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr bool isnormal(Vec<Dim, T> const& v) noexcept
{
    return std::ranges::all_of(v, [](T x) { return std::isnormal(x); });
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> lerp(Vec<Dim, T> a, Vec<Dim, T> const& b,
                                         T t) noexcept
{
    std::ranges::transform(a, b, a.begin(), [t](T x, T y) { return std::lerp(x, y, t); });
    return a;
}
 
template <std::size_t Dim, class T, class U>
[[nodiscard]] constexpr Vec<Dim, T> mix(Vec<Dim, T> const& x, Vec<Dim, T> const& y,
                                        Vec<Dim, U> const& a) noexcept
{
    return Vec<Dim, T>(Vec<Dim, U>(x) * (U{1} - a) + Vec<Dim, U>(y) * a);
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> mix(Vec<Dim, T> const& x, Vec<Dim, T> const& y,
                                        Vec<Dim, bool> const& a) noexcept
{
    Vec<Dim, T> res;
    for (auto [ri, xi, yi, ai] : std::views::zip(res, x, y, a)) {
        ri = ai ? yi : xi;
    }
    return res;
}
 
template <std::size_t Dim, class T>
[[nodiscard]] constexpr Vec<Dim, T> sign(Vec<Dim, T> const& x) noexcept
{
    Vec<Dim, T> res;
    std::ranges::transform(x, res.begin(), [](T v) { return static_cast<T>(sign(v)); });
    return res;
}
 
template <std::size_t Dim, std::floating_point T>
[[nodiscard]] constexpr Vec<Dim, T> ldexp(Vec<Dim, T> v, int exp) noexcept
{
    std::ranges::transform(v, v.begin(), [exp](T x) { return std::ldexp(x, exp); });
    return v;
}
 
/**************************************************************************************
|                                                                                     |
|                               Structured bindings                                   |
|                                                                                     |
**************************************************************************************/
 
template <std::size_t I, std::size_t Dim, class T>
  requires(I < Dim)
[[nodiscard]] constexpr T& get(Vec<Dim, T>& v) noexcept
{
    return v[I];
}
 
template <std::size_t I, std::size_t Dim, class T>
  requires(I < Dim)
[[nodiscard]] constexpr T const& get(Vec<Dim, T> const& v) noexcept
{
    return v[I];
}
 
template <std::size_t I, std::size_t Dim, class T>
  requires(I < Dim)
[[nodiscard]] constexpr T&& get(Vec<Dim, T>&& v) noexcept
{
    return std::move(v[I]);
}
 
/**************************************************************************************
|                                                                                     |
|                                       Print                                         |
|                                                                                     |
**************************************************************************************/
 
template <std::size_t Dim, class T>
std::ostream& operator<<(std::ostream& out, Vec<Dim, T> const& v)
{
    static constexpr std::array names = {'x', 'y', 'z', 'w'};
#ifdef __cpp_lib_ranges_enumerate
    for (auto const [i, name] : std::views::enumerate(names)) {
        if (i) out << ' ';
        out << name << ": " << v[i];
    }
#else
    for (std::size_t i{}; std::min(Dim, names.size()) > i; ++i) {
        if (i) out << ' ';
        out << names[i] << ": " << v[i];
    }
#endif
    return out;
}
}  // namespace ufo
 
template <std::size_t Dim, class T>
struct std::tuple_size<ufo::Vec<Dim, T>> : std::integral_constant<std::size_t, Dim> {
};
 
template <std::size_t I, std::size_t Dim, class T>
struct std::tuple_element<I, ufo::Vec<Dim, T>> {
    using type = T;
};
 
template <std::size_t Dim, class T>
  requires std::formattable<T, char>
struct std::formatter<ufo::Vec<Dim, T>> {
    constexpr auto parse(std::format_parse_context& ctx) { return ctx.begin(); }
 
    auto format(ufo::Vec<Dim, T> const& v, std::format_context& ctx) const
    {
        static constexpr std::array names = {'x', 'y', 'z', 'w'};
        auto                        out   = ctx.out();
#ifdef __cpp_lib_ranges_enumerate
        for (auto const [i, name] : std::views::enumerate(names)) {
            if (i) out = std::format_to(out, " ");
            out = std::format_to(out, "{}: {}", name, v[i]);
        }
#else
        for (std::size_t i{}; std::min(Dim, names.size()) > i; ++i) {
            if (i) out = std::format_to(out, " ");
            out = std::format_to(out, "{}: {}", names[i], v[i]);
        }
#endif
        return out;
    }
};
 
#endif  // UFO_NUMERIC_VEC_HPP