vec.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_VECTOR_H
#define THEORETICA_VECTOR_H
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include "../core/error.h"
#include "../core/real_analysis.h"
#include "./algebra.h"
#include <vector>
 
 
namespace theoretica {
 
    
    template<typename Vector, typename ReturnType = vector_element_t<Vector>&>
    class vec_iterator {
 
        private:
 
            Vector& vector;
 
            size_t i;
        
        public:
 
            using iterator_category = std::forward_iterator_tag;
            using value_type = vector_element_t<Vector>;
            using pointer = value_type*;
            using reference = value_type&;
 
            vec_iterator(Vector& vector, size_t index) : vector(vector), i(index) {}
 
            ReturnType operator*() {
                return vector[i];
            }
 
            vec_iterator& operator++() {
                ++i;
                return *this;
            }
 
            // vec_iterator& operator--() {
            //  --i;
            //  return *this;
            // }
 
 
            size_t index() {
                return i;
            }
 
 
            bool operator==(const vec_iterator& other) const {
                return i == other.i;
            }
 
            bool operator!=(const vec_iterator& other) const {
                return !(*this == other);
            }
    };
 
 
    template<typename Type = real, unsigned int N = 0>
    class vec {
        
        private:
            Type elements[N];
 
        public:
 
        // Vector size template argument
        static constexpr size_t size_argument = N;
 
        vec() {
            algebra::vec_zeroes(*this);
        }
 
 
        vec(Type val) {
 
            for (unsigned int i = 0; i < N; ++i)
                elements[i] = val;
        }
 
 
        vec(unsigned int size, Type val) {
 
            if (N != size) {
                TH_MATH_ERROR("vec::vec(size, val)", size, MathError::InvalidArgument);
                algebra::vec_error(*this);
                return;
            }
 
            for (unsigned int i = 0; i < N; ++i)
                elements[i] = val;
        }
 
 
        template <
            typename Vector,
            enable_vector<Vector> = true
        >
        vec(const Vector& other) {
            algebra::vec_copy(*this, other);
        }
 
 
        template<typename... Args>
        vec(Type x1, Type x2, Args... args) {
 
            static_assert(
                2 + sizeof...(args) == N,
                "Number of arguments must match vector size"
            );
 
            *this = {x1, x2, args...};
        }
 
 
        template<typename Vector>
        vec<Type, N>& operator=(const Vector& other) {
 
            if (N == other.size())
                return algebra::vec_copy(*this, other);
            else
                return algebra::vec_error(*this);
        }
 
 
        vec(std::initializer_list<Type> l) {
 
            if(l.size() != N) {
                TH_MATH_ERROR("vec::vec(initializer_list<Type>)", l.size(), MathError::InvalidArgument);
                algebra::vec_error(*this);
                return;
            }
 
            std::copy(l.begin(), l.end(), &elements[0]);
        }
 
        ~vec() = default;
 
 
        inline vec<Type, N> operator+() const {
            return *this;
        }
 
 
        inline vec<Type, N> operator+(const vec<Type, N>& other) const {
            
            vec<Type, N> result;
            algebra::vec_sum(result, *this, other);
            return result;
        }
 
 
        inline vec<Type, N> operator-() const {
            return *this * (Type) -1;
        }
 
 
        inline vec<Type, N> operator-(const vec<Type, N>& other) const {
            
            vec<Type, N> result;
            algebra::vec_diff(result, *this, other);
            return result;
        }
 
 
        inline vec<Type, N> operator*(Type scalar) const {
            vec<Type, N> result;
 
            for (unsigned int i = 0; i < N; ++i)
                result.elements[i] = scalar * elements[i];
 
            return result;
        }
 
 
        inline vec<Type, N> operator/(Type scalar) const {
            vec<Type, N> result;
 
            for (unsigned int i = 0; i < N; ++i)
                result.elements[i] = elements[i] / scalar;
 
            return result;
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        inline Type dot(const Vector& other) const {
            return algebra::dot(*this, other);
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        inline Type operator*(const Vector& other) const {
            return dot(other);
        }
 
 
        inline vec<Type, N> cross(const vec<Type, N>& other) const {
            static_assert(N == 3, "The vector must be three dimensional");
            return algebra::cross(*this, other);
        }
 
 
        template<typename Vector>
        inline vec<Type, N> cross(const Vector& other) const {
            return algebra::cross(*this, other);
        }
 
 
        template<typename Vector>
        inline vec<Type, N>& operator+=(const Vector& other) {
 
            for (unsigned int i = 0; i < N; ++i)
                elements[i] += other.elements[i];
        
            return *this;
        }
 
 
        template<typename Vector>
        inline vec<Type, N>& operator-=(const Vector& other) {
 
            for (unsigned int i = 0; i < N; ++i)
                elements[i] -= other.elements[i];
        
            return *this;
        }
 
 
        inline vec<Type, N>& operator*=(Type scalar) {
 
            for (unsigned int i = 0; i < N; ++i)
                elements[i] *= scalar;
        
            return *this;
        }
 
 
        inline vec<Type, N>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("vec::operator/=", scalar, MathError::DivByZero);
                algebra::vec_error(*this);
                return *this;
            }
 
            for (unsigned int i = 0; i < N; ++i)
                elements[i] /= scalar;
        
            return *this;
        }
 
 
        inline Type norm() const {
            return algebra::norm(*this);
        }
 
 
        inline Type sqr_norm() const {
            return algebra::sqr_norm(*this);
        }
 
 
        inline Type& operator[](unsigned int i) {
            return elements[i];
        }
 
 
        inline const Type& operator[](unsigned int i) const {
            return elements[i];
        }
 
 
        inline Type& at(unsigned int i) {
 
            if (i >= N) {
                throw std::out_of_range(
                    "The element index in vec::at() is out of bounds"
                );
            }
 
            return elements[i];
        }
 
 
        inline Type at(unsigned int i) const {
 
            if (i >= N) {
                throw std::out_of_range(
                    "The element index in vec::at() is out of bounds"
                );
            }
 
            return elements[i];
        }
 
 
        using iterator = vec_iterator<vec<Type, N>, Type&>;
 
        using const_iterator = vec_iterator<const vec<Type, N>, const Type&>;
 
 
        inline auto begin() {
            return iterator(*this, 0);
        }
 
 
        inline auto end() {
            return iterator(*this, size());
        }
 
 
        inline auto begin() const {
            return const_iterator(*this, 0);
        }
 
 
        inline auto end() const {
            return const_iterator(*this, size());
        }
 
 
        inline Type* data() {
            return elements;
        }
 
 
        inline const Type* data() const {
            return elements;
        }
 
 
        inline void normalize() {
            algebra::make_normalized(*this);
        }
 
 
        inline vec<Type, N> normalized() const {
            return algebra::normalize(*this);
        }
 
 
        template<typename Vector>
        inline bool operator==(const Vector& other) const {
 
            if(size() != other.size())
                return false;
 
            for (unsigned int i = 0; i < N; ++i)
                if(elements[i] != other[i])
                    return false;
 
            return true;
        }
 
 
        template<typename Vector>
        inline bool operator!=(const Vector& other) const {
            return !(*this == other);
        }
 
 
        inline TH_CONSTEXPR unsigned int size() const {
            return N;
        }
 
 
        inline void resize(size_t n) const {
            
            if(N != n) {
                TH_MATH_ERROR("vec::resize", N, MathError::InvalidArgument);
            }
        }
 
 
        inline void resize(size_t n, const Type& value) const {
            
            if(N != n) {
                TH_MATH_ERROR("vec::resize", N, MathError::InvalidArgument);
            }
        }
 
 
        inline static vec<Type, N> euclidean_base(
            unsigned int i, unsigned int n = N) {
 
            if(i >= n) {
                TH_MATH_ERROR("vec::euclidean_base", i, MathError::InvalidArgument);
                return make_error<vec<Type, N>>(n);
            }
 
            vec<Type, N> e_i = vec<Type, N>(n, Type(0.0));
            e_i.resize(n);
            e_i[i] = 1;
 
            return e_i;
        }
 
 
        inline friend vec<Type, N> operator*(Type a, const vec<Type, N>& v) {
            return v * a;
        }
 
 
#ifndef THEORETICA_NO_PRINT
 
        inline std::string to_string(
            const std::string& separator = ", ",
            bool parenthesis = true) const {
 
            std::stringstream res;
 
            if(parenthesis)
                res << "(";
            
            for (unsigned int i = 0; i < N; ++i) {
                res << elements[i];
                if(i != N - 1)
                    res << separator;
            }
 
            if(parenthesis)
                res << ")";
 
            return res.str();
        }
 
 
        inline operator std::string() {
            return to_string();
        }
 
 
        inline friend std::ostream& operator<<(
            std::ostream& out, const vec<Type, N>& obj) {
            return out << obj.to_string();
        }
 
#endif
 
    };
 
 
    template<typename Type>
    class vec<Type, 0> {
 
        // Container type for storage (alias for std::vector)
        template<typename T>
        using Container = std::vector<T>;
 
    private:
            Container<Type> elements;
 
    public:
 
        // Vector size template argument
        static constexpr size_t size_argument = 0;
 
        vec() = default;
 
        vec(unsigned int n) {
            resize(n);
            algebra::vec_zeroes(*this);
        }
 
        vec(unsigned int n, Type a) {
            elements = std::vector<Type>(n, a);
        }
 
        template <
            typename Vector,
            enable_vector<Vector> = true
        >
        vec(const Vector& other) {
            algebra::vec_copy(*this, other);
        }
 
        template<typename... Args>
        vec(Type x1, Type x2, Args... args) {
 
            elements = {x1, x2, args...};
        }
 
        
        vec(std::initializer_list<Type> l) : elements(l) {}
 
 
        inline vec<Type> operator+() const {
            return *this;
        }
 
 
        template<typename Vector>
        inline vec<Type> operator+(const Vector& other) const {
            
            vec<Type> result;
            result.resize(size());
            algebra::vec_sum(result, *this, other);
            return result;
        }
 
 
        inline vec<Type> operator-() const {
            return *this * (Type) -1;
        }
 
 
        template<typename Vector>
        inline vec<Type> operator-(const Vector& other) const {
            
            vec<Type> result;
            result.resize(size());
            algebra::vec_diff(result, *this, other);
            return result;
        }
 
 
        inline vec<Type> operator*(Type scalar) const {
 
            vec<Type> result;
            result.resize(size());
 
            for (unsigned int i = 0; i < size(); ++i)
                result.elements[i] = scalar * elements[i];
 
            return result;
        }
 
 
        inline vec<Type> operator/(Type scalar) const {
 
            vec<Type> result;
            result.resize(size());
 
            for (unsigned int i = 0; i < size(); ++i)
                result.elements[i] = elements[i] / scalar;
 
            return result;
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        inline Type dot(const Vector& other) const {
            return algebra::dot(*this, other);
        }
 
 
        template<typename Vector, enable_vector<Vector> = true>
        inline Type operator*(const Vector& other) const {
            return dot(other);
        }
 
 
        template<typename Vector>
        inline vec<Type> cross(const Vector& other) const {
            return algebra::cross(*this, other);
        }
 
 
        template<typename Vector>
        inline vec<Type>& operator+=(const Vector& other) {
 
            // If the vector is uninitialized,
            // initialize it to be a zero vector
            // with the target size
            if(size() == 0)
                resize(other.size());
 
            if(size() != other.size()) {
                TH_MATH_ERROR("vec::operator+=", size(), MathError::InvalidArgument);
                return (*this = vec<Type>(max(size(), 1), nan()));
            }
 
            for (unsigned int i = 0; i < size(); ++i)
                elements[i] += other.elements[i];
        
            return *this;
        }
 
 
        template<typename Vector>
        inline vec<Type>& operator-=(const Vector& other) {
 
            if(size() != other.size()) {
                TH_MATH_ERROR("vec::operator-=", size(), MathError::InvalidArgument);
                return (*this = vec<Type>(max(size(), 1), nan()));
            }
 
            for (unsigned int i = 0; i < size(); ++i)
                elements[i] -= other.elements[i];
        
            return *this;
        }
 
 
        inline vec<Type>& operator*=(Type scalar) {
 
            for (unsigned int i = 0; i < size(); ++i)
                elements[i] *= scalar;
        
            return *this;
        }
 
 
        inline vec<Type>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("vec::operator/=", scalar, MathError::DivByZero);
                *this = vec<Type>(max(size(), 1), nan());
                return *this;
            }
 
            for (unsigned int i = 0; i < size(); ++i)
                elements[i] /= scalar;
        
            return *this;
        }
 
 
        inline Type norm() const {
            return algebra::norm(*this);
        }
 
 
        inline Type sqr_norm() const {
            return algebra::sqr_norm(*this);
        }
 
 
        inline Type& operator[](unsigned int i) {
            return elements[i];
        }
 
 
        inline const Type& operator[](unsigned int i) const {
            return elements[i];
        }
 
 
        inline Type& at(unsigned int i) {
            return elements.at(i);
        }
 
 
        inline Type at(unsigned int i) const {
            return elements.at(i);
        }
 
 
        using iterator = typename Container<Type>::iterator;
 
 
        inline auto begin() {
            return elements.begin();
        }
 
 
        inline auto begin() const {
            return elements.cbegin();
        }
 
 
        inline auto end() {
            return elements.end();
        }
 
 
        inline auto end() const {
            return elements.cend();
        }
 
 
        inline Type* data() {
            return elements.data();
        }
 
 
        inline const Type* data() const {
            return elements.data();
        }
 
 
        inline void normalize() {
            algebra::make_normalized(*this);
        }
 
 
        inline vec<Type> normalized() const {
            return algebra::normalize(*this);
        }
 
 
        template<typename Vector>
        inline bool operator==(const Vector& other) const {
 
            if(size() != other.size())
                return false;
 
            for (unsigned int i = 0; i < size(); ++i)
                if(elements[i] != other[i])
                    return false;
 
            return true;
        }
 
 
        template<typename Vector>
        inline bool operator!=(const Vector& other) const {
            return !(*this == other);
        }
 
 
        inline TH_CONSTEXPR unsigned int size() const {
            return elements.size();
        }
 
 
        inline void resize(size_t n) {
            elements.resize(n);
        }
 
 
        inline void resize(size_t n, const Type& value) {
            elements.resize(n, value);
        }
 
 
        inline void append(const Type& x) {
            elements.push_back(x);
        }
 
 
        inline void append(Type&& x) {
            elements.push_back(x);
        }
 
 
        inline static vec<Type> euclidean_base(
            unsigned int i, unsigned int n) {
 
            if(i >= n) {
                TH_MATH_ERROR("vec::euclidean_base", i, MathError::InvalidArgument);
                return vec<Type>(n, nan());
            }
 
            vec<Type> e_i = vec<Type>(n, Type(0.0));
            e_i.resize(n);
            e_i[i] = 1;
 
            return e_i;
        }
 
 
        inline friend vec<Type> operator*(Type a, const vec<Type>& v) {
            return v * a;
        }
 
 
#ifndef THEORETICA_NO_PRINT
 
        inline std::string to_string(
            const std::string& separator = ", ",
            bool parenthesis = true) const {
 
            std::stringstream res;
 
            if(parenthesis)
                res << "(";
            
            for (unsigned int i = 0; i < size(); ++i) {
                res << elements[i];
                if(i != size() - 1)
                    res << separator;
            }
            
            if(parenthesis)
                res << ")";
 
            return res.str();
        }
 
 
        inline operator std::string() {
            return to_string();
        }
 
 
        inline friend std::ostream& operator<<(std::ostream& out, const vec<Type>& obj) {
            return out << obj.to_string();
        }
 
#endif
 
    };
 
}
 
 
#endif