#include <iostream>

#include <cassert>
#include <cmath>

#include <array>
#include <algorithm>
#include <numeric>
#include <initializer_list>


template<typename T, std::size_t d>
struct basic_vec : public std::array<T, d> {
    static constexpr std::size_t dimensions = d;

    static constexpr T null_element = static_cast<T>(0);
    static constexpr T unit_element = static_cast<T>(1);
    static constexpr T unit_additive_inverse_element = static_cast<T>(-1);

    basic_vec();
    basic_vec(const std::initializer_list<T> l);
    template<std::size_t n> basic_vec(const basic_vec<T, n>& other);

    T length() const;
};

template<typename T, std::size_t d>
basic_vec<T, d>::basic_vec() : std::array<T, d>() {
    this->fill(basic_vec<T,d>::null_element);
}

template<typename T, std::size_t d>
basic_vec<T, d>::basic_vec(const std::initializer_list<T> l) {
    basic_vec();

    // why can't this sh*t be a constexpr with static_assert???
    assert(l.size() <= d);
    std::copy(l.begin(), l.end(), this->begin());
}

template<typename T, std::size_t d>
template<std::size_t n>
basic_vec<T, d>::basic_vec(const basic_vec<T, n>& other) {
    basic_vec();

    static_assert(d >= n);
    for (std::size_t i = 0;  i < n; i++) {
        this->at(i) = other.at(i);
    }
}

template<typename T, std::size_t d>
T basic_vec<T, d>::length() const {
    T res = basic_vec<T, d>::null_element;
    for (const T& val : *this) {
        res += val * val;
    }

    return res;
}

template<typename T, std::size_t d>
basic_vec<T, d> operator+(const basic_vec<T, d>& rhs, const basic_vec<T, d>& lhs) {
    basic_vec<T, d> out;
    
    std::transform(rhs.begin(), rhs.end(), lhs.begin(), out.begin(),
        [](const T& r, const T& l) -> T {
            return r + l;
        }
    );

    return out;
}

template<typename T, std::size_t d>
basic_vec<T, d> operator*(const T& rhs, const basic_vec<T, d>& lhs) {
    basic_vec<T, d> out;

    std::transform(lhs.begin(), lhs.end(), out.begin(), 
        [rhs](const T& t) -> T {
            return t * rhs;
    });

    return out;
}

template<typename T, std::size_t d>
basic_vec<T, d> operator-(const basic_vec<T, d>& rhs, const basic_vec<T, d>& lhs) {
    return rhs + basic_vec<T, d>::unit_additive_inverse_element * lhs;
}

template<typename T, std::size_t d>
T operator*(const basic_vec<T, d>& rhs, const basic_vec<T, d>& lhs) {
    return std::inner_product(rhs.begin(), rhs.end(), lhs.begin(), 0);
}

template<typename T, std::size_t d>
std::ostream& operator<<(std::ostream& os, const basic_vec<T, d>& v) {
    os << "<";
    for (std::size_t i = 0; i < d -1; i++) {
        os << v[i] << ", ";
    }
    os << v[d-1] << ">";
    return os;
}


template<typename T, std::size_t d>
class vec: public basic_vec<T, d> {
public:
    vec(std::initializer_list<T> l) : basic_vec<T, d>(l) {}

    template<std::size_t n>
    vec(const basic_vec<T, n>& other) : basic_vec<T, d>(other) {}
};

template<typename T>
class vec3 : public basic_vec<T, 3> {
public:
    vec3(std::initializer_list<T> l) : basic_vec<T, 3>(l) {}

    template<std::size_t n>
    vec3(const basic_vec<T, n>& other) : basic_vec<T, 3>(other) {}

    static vec3<T> cross(const vec3<T>& rhs, const vec3<T>& lhs);
};

template<typename T>
vec3<T> vec3<T>::cross(const vec3<T>& rhs, const vec3<T>& lhs) {
    vec3<T> res;

    res[0] = (rhs[1] * lhs[2]) - (rhs[2] * lhs[1]);
    res[1] = (rhs[2] * lhs[0]) - (rhs[0] * lhs[2]);
    res[2] = (rhs[0] * lhs[1]) - (rhs[1] * lhs[0]);

    return res;
}


template<typename T>
class vec2: public basic_vec<T, 2> {
public:
    vec2(std::initializer_list<T> l) : basic_vec<T, 2>(l) {}

    template<std::size_t n>
    vec2(const basic_vec<T, n>& other) : basic_vec<T, 2>(other) {}

    T polar();
    static vec3<T> cross(const vec2<T>& rhs, const vec2<T>& lhs);
};

template<typename T>
T vec2<T>::polar() {
    return std::atan2(this->at(0), this->at(1));
}

template<typename T>
vec3<T> vec2<T>::cross(const vec2<T>& rhs, const vec2<T>& lhs) {
    return vec3<T>::cross(vec3<T>(rhs), vec3<T>(lhs));
}


int main(int argc, char *argv[]) {
    vec3<double> v{1, 2, 3};
    vec3<double> u{3, 4, 5};

    std::cout << v << std::endl;
    std::cout << v.length() << std::endl;
    std::cout << v + u << std::endl;
    std::cout << v - u << std::endl;
    std::cout << 2.0 * u << std::endl;
    std::cout << v * u << std::endl;

    return 0;
}