vector.hpp

Go to the documentation of this file.

 
#ifndef NUMERICAL_VECTOR_HPP
#define NUMERICAL_VECTOR_HPP
 
#include <algorithm>
#include <cmath>
#include <complex> // so std::abs resolves for a complex scalar, regardless of include order
#include <cstddef>
#include <functional>
#include <initializer_list>
#include <limits>
#include <memory>
#include <numeric>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
 
#include "buffer.hpp"
 
namespace numerical {
 
  class dimension_error : public std::runtime_error
  {
  public:
 
    explicit dimension_error(const std::string& message)
      : std::runtime_error(message)
    {}
  };
 
  namespace detail {
 
    [[nodiscard]] inline std::size_t checked_mul(std::size_t a,
                                                 std::size_t b)
    {
      if (a != 0 && b > (std::numeric_limits<std::size_t>::max() / a)) {
        throw dimension_error("dimension product overflows std::size_t");
      }
      return a * b;
    }
 
    template <typename It>
    [[nodiscard]] std::size_t element_count(It first,
                                            It last)
    {
      if (std::find(first, last, std::size_t {0}) != last) {
        return 0;
      }
      std::size_t count {1};
      for (; first != last; ++first) {
        count = checked_mul(count, *first);
      }
      return count;
    }
  } // namespace detail
 
  template <typename T>
  class vector
  {
  public:
 
    vector() = default;
 
    vector(std::initializer_list<T> elements)
      : data_(elements)
    {}
 
    explicit vector(buffer<T> elements)
      : data_(std::move(elements))
    {}
 
    explicit vector(const std::vector<T>& elements)
      : data_(elements.begin(), elements.end())
    {}
 
    [[nodiscard]] std::size_t size() const noexcept
    {
      return data_.size();
    }
 
    [[nodiscard]] const T& operator[](std::size_t index) const
    {
      return data_[index];
    }
 
    [[nodiscard]] typename buffer<T>::const_iterator begin() const noexcept
    {
      return data_.begin();
    }
 
    [[nodiscard]] typename buffer<T>::const_iterator end() const noexcept
    {
      return data_.end();
    }
 
    [[nodiscard]] bool operator==(const vector& other) const
    {
      return data_ == other.data_;
    }
 
  private:
 
    buffer<T> data_; 
  };
 
  template <typename T>
  [[nodiscard]] vector<T> operator+(const vector<T>& a,
                                    const vector<T>& b)
  {
    if (a.size() != b.size()) {
      throw dimension_error("vector sizes differ in '+'");
    }
    buffer<T> out(a.size());
    std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::plus<T> {});
    return vector<T> {std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] vector<T> operator-(const vector<T>& a,
                                    const vector<T>& b)
  {
    if (a.size() != b.size()) {
      throw dimension_error("vector sizes differ in '-'");
    }
    buffer<T> out(a.size());
    std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::minus<T> {});
    return vector<T> {std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] vector<T> operator*(const T&         scalar,
                                    const vector<T>& v)
  {
    buffer<T> out(v.size());
    std::transform(v.begin(), v.end(), out.begin(),
                   [&scalar](const T& element) { return scalar * element; });
    return vector<T> {std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] vector<T> operator*(const vector<T>& v,
                                    const T&         scalar)
  {
    return scalar * v;
  }
 
  template <typename T>
  [[nodiscard]] T dot(const vector<T>& a,
                      const vector<T>& b)
  {
    if (a.size() != b.size()) {
      throw dimension_error("vector sizes differ in 'dot'");
    }
    // std::inner_product is sequential (left-to-right), so the reduction order —
    // and thus the floating-point result — stays deterministic.
    return std::inner_product(a.begin(), a.end(), b.begin(), T {});
  }
 
  template <typename T>
  [[nodiscard]] bool approx_equal(const vector<T>& a,
                                  const vector<T>& b,
                                  double           tolerance = 1e-9)
  {
    if (a.size() != b.size()) {
      return false;
    }
    for (std::size_t i {0}; i < a.size(); ++i) {
      // Equal, or within the tolerance: a NaN is neither, so it is never close, and two equal
      // infinities are equal although their difference is NaN.
      if (!(a[i] == b[i]) && !(std::abs(a[i] - b[i]) <= tolerance)) {
        return false;
      }
    }
    return true;
  }
} // namespace numerical
 
#endif // NUMERICAL_VECTOR_HPP