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 data[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)
                data[i] = val;
        }
 
 
        vec(unsigned int size, Type val) {
 
            resize(size);
 
            for (unsigned int i = 0; i < N; ++i)
                data[i] = val;
        }
 
 
        template<unsigned int M>
        vec(const vec<Type, M>& other) {
            algebra::vec_copy(*this, other);
        }
 
 
        template<typename Vector>
        vec<Type, N>& operator=(const Vector& other) {
            return algebra::vec_copy(*this, other);
        }
 
 
        vec(std::initializer_list<Type> l) {
 
            if(l.size() != N) {
                TH_MATH_ERROR("vec::vec(initializer_list<Type>)", l.size(),
                    INVALID_ARGUMENT);
                // Set all elements to NaN
                *this = vec<Type, N>(Type(nan()));
                return;
            }
 
            std::copy(l.begin(), l.end(), &data[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.data[i] = scalar * data[i];
 
            return result;
        }
 
 
        inline vec<Type, N> operator/(Type scalar) const {
            vec<Type, N> result;
 
            for (unsigned int i = 0; i < N; ++i)
                result.data[i] = data[i] / scalar;
 
            return result;
        }
 
 
        template<typename Vector>
        inline Type dot(const Vector& other) const {
            return algebra::dot(*this, other);
        }
 
 
        template<typename Vector>
        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 {
        
            if(other.size() != 3) {
                TH_MATH_ERROR("vec::cross", other.size(), INVALID_ARGUMENT);
                return vec<Type, N>(Type(nan()));
            }
 
            return algebra::cross(*this, other);
        }
 
 
        template<typename Vector>
        inline vec<Type, N>& operator+=(const Vector& other) {
 
            for (unsigned int i = 0; i < N; ++i)
                data[i] += other.data[i];
        
            return *this;
        }
 
 
        template<typename Vector>
        inline vec<Type, N>& operator-=(const Vector& other) {
 
            for (unsigned int i = 0; i < N; ++i)
                data[i] -= other.data[i];
        
            return *this;
        }
 
 
        inline vec<Type, N>& operator*=(Type scalar) {
 
            for (unsigned int i = 0; i < N; ++i)
                data[i] *= scalar;
        
            return *this;
        }
 
 
        inline vec<Type, N>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("vec::operator/=", scalar, DIV_BY_ZERO);
                vec_error(*this);
                return *this;
            }
 
            for (unsigned int i = 0; i < N; ++i)
                data[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 data[i];
        }
 
 
        inline const Type& operator[](unsigned int i) const {
            return data[i];
        }
 
 
        inline Type& at(unsigned int i) {
            return data[i];
        }
 
 
        inline Type get(unsigned int i) const {
            return data[i];
        }
 
 
        inline void set(unsigned int i, Type x) {
            data[i] = x;
        }
 
 
        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 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(data[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, INVALID_ARGUMENT);
            }
        }
 
 
        inline static vec<Type, N> euclidean_base(
            unsigned int i, unsigned int n = N) {
 
            if(i >= n) {
                TH_MATH_ERROR("vec::euclidean_base", i, INVALID_ARGUMENT);
                return vec<Type, N>(Type(nan()));
            }
 
            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 << data[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> data;
 
    public:
 
        // Vector size template argument
        static constexpr size_t size_argument = 0;
 
        vec() {}
 
 
        vec(unsigned int n) {
            resize(n);
            algebra::vec_zeroes(*this);
        }
 
 
        vec(unsigned int n, Type a) {
            data = std::vector<Type>(n, a);
        }
 
 
        template<unsigned int M>
        vec(const vec<Type, M>& other) {
            algebra::vec_copy(*this, other);
        }
 
 
        template<typename Vector>
        vec<Type>& operator=(const Vector& other) {
            return algebra::vec_copy(*this, other);
        }
 
 
        vec(std::initializer_list<Type> l) {
 
            resize(l.size());
            std::copy(l.begin(), l.end(), &data[0]);
        }
 
        ~vec() = default;
 
 
        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.data[i] = scalar * data[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.data[i] = data[i] / scalar;
 
            return result;
        }
 
 
        template<typename Vector>
        inline Type dot(const Vector& other) const {
            return algebra::dot(*this, other);
        }
 
 
        template<typename Vector>
        inline Type operator*(const Vector& other) const {
            return dot(other);
        }
 
 
        template<typename Vector>
        inline vec<Type> cross(const Vector& other) const {
        
            if(other.size() != 3) {
                TH_MATH_ERROR("vec::cross", other.size(), INVALID_ARGUMENT);
                return vec<Type>(3, nan());
            }
 
            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(), INVALID_ARGUMENT);
                return (*this = vec<Type>(max(size(), 1), nan()));
            }
 
            for (unsigned int i = 0; i < size(); ++i)
                data[i] += other.data[i];
        
            return *this;
        }
 
 
        template<typename Vector>
        inline vec<Type>& operator-=(const Vector& other) {
 
            if(size() != other.size()) {
                TH_MATH_ERROR("vec::operator-=", size(), INVALID_ARGUMENT);
                return (*this = vec<Type>(max(size(), 1), nan()));
            }
 
            for (unsigned int i = 0; i < size(); ++i)
                data[i] -= other.data[i];
        
            return *this;
        }
 
 
        inline vec<Type>& operator*=(Type scalar) {
 
            for (unsigned int i = 0; i < size(); ++i)
                data[i] *= scalar;
        
            return *this;
        }
 
 
        inline vec<Type>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("vec::operator/=", scalar, DIV_BY_ZERO);
                *this = vec<Type>(max(size(), 1), nan());
                return *this;
            }
 
            for (unsigned int i = 0; i < size(); ++i)
                data[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 data[i];
        }
 
 
        inline const Type& operator[](unsigned int i) const {
            return data[i];
        }
 
 
        inline Type& at(unsigned int i) {
            return data[i];
        }
 
 
        inline Type get(unsigned int i) const {
            return data[i];
        }
 
 
        inline void set(unsigned int i, Type x) {
            data[i] = x;
        }
 
 
        using iterator = typename Container<Type>::iterator;
 
 
        inline auto begin() {
            return data.begin();
        }
 
 
        inline auto begin() const {
            return data.cbegin();
        }
 
 
        inline auto end() {
            return data.end();
        }
 
 
        inline auto end() const {
            return data.cend();
        }
 
 
        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(data[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 data.size();
        }
 
 
        inline void resize(size_t n) {
            data.resize(n);
        }
 
 
        inline void push(const Type& x) {
            data.push_back(x);
        }
 
 
        inline void push(Type&& x) {
            data.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, INVALID_ARGUMENT);
                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 << data[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
 
    };
 
 
    template<typename ElementType, typename Type, typename ...Args>
    void make_vec(vec<ElementType>& v, size_t index, Type last) {
        v[index] = last;
    }
 
 
    template<typename ElementType, typename Type, typename ...Args>
    void make_vec(vec<ElementType>& v, size_t index, Type first, Args... elements) {
 
        v[index] = first;
        make_vec<ElementType>(v, index + 1, elements...);
    }
 
 
    template<typename Type, typename ...Args>
    vec<Type> make_vec(Type first, Args... elements) {
 
        vec<Type> v;
        v.resize(sizeof...(elements) + 1);
 
        v[0] = first;
        make_vec<Type>(v, 1, elements...);
 
        return v;
    }
 
}
 
 
#endif