vector.hpp

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <complex>
#include <cstddef>
#include <initializer_list>
#include <iomanip>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <type_traits>
 
namespace firefly {
 
template <typename T1, typename T2>
struct common_type : std::common_type<T1, T2> {};
 
template <typename T1, typename T2>
struct common_type<T1, std::complex<T2>> {
  using type = std::complex<typename std::common_type_t<T1, T2>>;
};
 
template <typename T1, typename T2>
struct common_type<std::complex<T1>, T2> {
  using type = std::complex<typename std::common_type_t<T1, T2>>;
};
 
template <typename T1, typename T2>
struct common_type<std::complex<T1>, std::complex<T2>> {
  using type = std::complex<typename std::common_type_t<T1, T2>>;
};
 
template <typename T1, typename T2>
using common_type_t = typename common_type<T1, T2>::type;
 
template <typename T>
concept complex_type = std::is_same_v<std::decay_t<T>, std::complex<typename T::value_type>> &&
                       std::is_arithmetic_v<typename T::value_type>;
 
template <typename T>
concept vector_type = std::is_arithmetic_v<T> || complex_type<T>;
 
template <typename T>
struct is_complex : std::false_type {};
 
template <typename T>
struct is_complex<std::complex<T>> : std::true_type {};
 
template <typename T>
inline constexpr bool is_complex_v = is_complex<T>::value;
 
template <vector_type T, std::size_t Length> //
class vector : private std::array<T, Length> {
 
public:
  using value_type = T;
  using std::array<T, Length>::begin;
  using std::array<T, Length>::end;
  using std::array<T, Length>::cbegin;
  using std::array<T, Length>::cend;
  using std::array<T, Length>::rbegin;
  using std::array<T, Length>::rend;
  using std::array<T, Length>::crbegin;
  using std::array<T, Length>::crend;
  using std::array<T, Length>::empty;
  using std::array<T, Length>::size;
  using std::array<T, Length>::operator[];
 
  [[nodiscard]] constexpr vector() : std::array<T, Length>() {
    std::fill(begin(), end(), T{});
  }
 
  [[nodiscard]] constexpr vector(std::initializer_list<T> const &list) : std::array<T, Length>() {
    if (list.size() > Length) {
      throw std::out_of_range("Initializer list size must match vector Length");
    }
    std::copy(list.begin(), list.end(), begin());
  }
 
  [[nodiscard]] constexpr vector(T const value) : std::array<T, Length>() {
    std::fill(begin(), end(), value);
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto add(vector<U, Length> const &other) const {
    vector<common_type_t<T, U>, Length> result;
    std::transform(cbegin(), cend(), other.cbegin(), result.begin(),
                   [](auto const &a, auto const &b) { return common_type_t<T, U>(a) + common_type_t<T, U>(b); });
 
    return result;
  }
 
  template <typename U>
  [[nodiscard]] constexpr auto add(U const scalar) const {
    vector<common_type_t<T, U>, Length> result;
    std::transform(cbegin(), cend(), result.begin(),
                   [scalar](auto const &a) { return common_type_t<T, U>(a) + common_type_t<T, U>(scalar); });
    return result;
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto operator+(vector<U, Length> const &other) const {
    return add(other);
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto operator+(U const scalar) const {
    return add(scalar);
  }
 
  template <vector_type U>
  friend constexpr auto operator+(U const scalar, vector<T, Length> const &vec) {
    return vec + scalar;
  }
 
  template <vector_type U>
  constexpr auto &operator+=(vector<U, Length> const &other) {
    std::transform(cbegin(), cend(), other.cbegin(), begin(),
                   [](auto const &a, auto const &b) { return common_type_t<T, U>(a) + common_type_t<T, U>(b); });
    return *this;
  }
 
  template <vector_type U>
  constexpr auto &operator+=(U const scalar) {
    std::transform(cbegin(), cend(), begin(),
                   [scalar](T const &el) { return common_type_t<T, U>(el) + common_type_t<T, U>(scalar); });
    return *this;
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto subtract(vector<U, Length> const &other) const {
    return add(-other);
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto subtract(U const scalar) const {
    return add(-scalar);
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto operator-(vector<U, Length> const &other) const {
    return subtract(other);
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto operator-(U const scalar) const {
    return subtract(scalar);
  }
 
  template <vector_type U>
  friend constexpr auto operator-(U const scalar, vector<T, Length> const &vec) {
    return vec - scalar;
  }
 
  template <vector_type U>
  constexpr auto &operator-=(vector<U, Length> const &other) {
    std::transform(cbegin(), cend(), other.cbegin(), begin(),
                   [](auto const &a, auto const &b) { return common_type_t<T, U>(a) - common_type_t<T, U>(b); });
    return *this;
  }
 
  template <vector_type U>
  constexpr auto &operator-=(U const scalar) {
    std::transform(cbegin(), cend(), begin(),
                   [scalar](auto const &el) { return common_type_t<T, U>(el) - common_type_t<T, U>(scalar); });
    return *this;
  }
 
  template <typename U>
  [[nodiscard]] constexpr auto dot(vector<U, Length> const &other) const {
    return std::transform_reduce(
        cbegin(), cend(), other.cbegin(), common_type_t<T, U>(0), std::plus<>(),
        [](auto const &a, auto const &b) { return common_type_t<T, U>(a) * common_type_t<T, U>(b); });
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto cross(vector<U, Length> const &other) const {
    static_assert(Length == 3, "Cross product is only allowed for 3D vectors.");
    vector<common_type_t<T, U>, Length> cross;
 
    cross[0] = (*this)[1] * other[2] - (*this)[2] * other[1];
    cross[1] = (*this)[2] * other[0] - (*this)[0] * other[2];
    cross[2] = (*this)[0] * other[1] - (*this)[1] * other[0];
 
    return cross;
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto scale(U const scalar) const {
    vector<common_type_t<T, U>, Length> result;
    std::transform(cbegin(), cend(), result.begin(),
                   [scalar](auto const &el) { return common_type_t<T, U>(el) * common_type_t<T, U>(scalar); });
    return result;
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto operator*(vector<U, Length> const &other) const {
    return dot(other);
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto operator*(U const scalar) const {
    return scale(scalar);
  }
 
  template <vector_type U>
  friend constexpr auto operator*(U const scalar, vector<T, Length> const &vec) {
    return vec * scalar;
  }
 
  template <vector_type U>
  constexpr auto &operator*=(U const scalar) {
    std::transform(cbegin(), cend(), begin(),
                   [scalar](auto const &el) { return common_type_t<T, U>(el) * common_type_t<T, U>(scalar); });
    return *this;
  }
 
  template <vector_type U>
  [[nodiscard]] constexpr auto operator/(U const scalar) const {
    return scale(1 / scalar);
  }
 
  template <vector_type U>
  friend constexpr auto operator/(U const scalar, vector<T, Length> const &vec) {
    return vec / scalar;
  }
 
  template <vector_type U>
  constexpr auto &operator/=(U const scalar) {
    std::transform(cbegin(), cend(), begin(),
                   [scalar](T const &el) { return common_type_t<T, U>(el) * (1 / common_type_t<T, U>(scalar)); });
    return *this;
  }
 
  [[nodiscard]] constexpr auto operator-() const {
    return scale(-1);
  }
 
  [[nodiscard]] constexpr auto norm() const {
    if constexpr (is_complex_v<T>) {
      return std::sqrt(std::transform_reduce(cbegin(), cend(), 0.0, std::plus<>(), [](const auto &val) {
        return std::norm(val); // |val|^2 = val * conj(val)
      }));
    } else {
      return std::sqrt(dot(*this));
    }
  }
 
  [[nodiscard]] constexpr auto to_normalized() const {
    auto _norm = norm();
    if (_norm == 0) {
      throw std::logic_error("zero norm results in divide by zero");
    }
    return scale(1 / norm());
  }
 
  [[nodiscard]] constexpr auto is_equal(vector<T, Length> const &other) const {
    return std::equal(cbegin(), cend(), other.cbegin());
  }
 
  [[nodiscard]] constexpr auto operator==(vector<T, Length> const &other) const {
    return is_equal(other);
  }
 
  template <vector_type AsType>
  [[nodiscard]] vector<AsType, Length> constexpr const as_type() const {
    vector<AsType, Length> result;
 
    std::transform(cbegin(), cend(), result.begin(), [](T el) {
      if constexpr (std::is_same_v<T, std::complex<typename T::value_type>>) {
        return static_cast<AsType>(el.real() * el.real() + el.imag() * el.imag());
      } else {
        return static_cast<AsType>(el);
      }
    });
 
    return result;
  }
 
  [[nodiscard]] constexpr std::string view(int precision = 20) const {
    bool f_is_first = false;
    std::stringstream ss;
 
    ss << std::setprecision(precision) << "[";
    if (!empty()) {
      for (auto const &el : *this) {
        if (!f_is_first) {
          ss << el;
          f_is_first = true;
        } else {
          ss << ", " << el;
        }
      }
    }
    ss << "]";
 
    return ss.str();
  }
 
  friend std::ostream &operator<<(std::ostream &os, vector const &other) {
    os << other.view();
    return os;
  }
};
 
} // namespace firefly