#include <cassert>
#include <iostream>
#include <cmath>

template <std::size_t N>
vec<N>::vec() {
    for ( std::size_t k{0}; k < N; ++k ) {
        entries_[k] = 0.0;
    }
}

template <std::size_t N>
vec<N>::vec( double s ) {
    for ( std::size_t k{0}; k < N; ++k ) {
        entries_[k] = s;
    }
}

template <std::size_t N>
vec<N>::vec( vec<N> const &v ) {
    for ( std::size_t k{0}; k < N; ++k ) {
        entries_[k] = v.entries_[k];
    }
}

template <std::size_t N>
vec<N>::vec( std::initializer_list<double> init ) {
    std::size_t k{0};

    for ( std::initializer_list<double>::iterator itr{init.begin()}; itr != init.end(); ++itr ) {
        entries_[k] = *itr;
        ++k;
    }

    assert( k <= N );

    for ( ; k < N; ++k ) {
        entries_[k] = 0.0;
    }
}

template <std::size_t N>
vec<N> vec<N>::operator =( vec<N> const &v ) {
    for ( std::size_t k{0}; k < N; ++k ) {
        entries_[k] = v.entries_[k];
    }

    return *this;
}

template <std::size_t N>
vec<N> &vec<N>::operator *=( double s ) {
    for ( std::size_t k{0}; k < N; ++k ) {
        entries_[k] *= s;
    }

    return *this;
}

template <std::size_t N>
vec<N> vec<N>::operator *( double s ) const {
    vec<N> ret;

    for ( std::size_t k{0}; k < N; ++k ) {
        ret.entries_[k] = entries_[k]*s;
    }

    return ret;
}

template <std::size_t N>
vec<N> &vec<N>::operator +=( vec const &v ) {
    for ( std::size_t k{0}; k < N; ++k ) {
        entries_[k] += v.entries_[k];
    }

    return *this;
}

template <std::size_t N>
vec<N> vec<N>::operator +( vec const &v ) const {
    vec<N> ret;

    for ( std::size_t k{0}; k < N; ++k ) {
        ret.entries_[k] = entries_[k] + v.entries_[k];
    }

    return ret;
}

template <std::size_t N>
vec<N> &vec<N>::operator -=( vec const &v ) {
    for ( std::size_t k{0}; k < N; ++k ) {
        entries_[k] -= v.entries_[k];
    }

    return *this;
}

template <std::size_t N>
vec<N> vec<N>::operator -( vec const &v ) const {
    vec<N> ret;

    for ( std::size_t k{0}; k < N; ++k ) {
        ret.entries_[k] = entries_[k] - v.entries_[k];
    }

    return ret;
}

template <std::size_t N>
vec<N> vec<N>::operator +() const {
    vec<N> ret;

    for ( std::size_t k{0}; k < N; ++k ) {
        ret.entries_[k] = entries_[k];
    }

    return ret;
}
template <std::size_t N>
vec<N> vec<N>::operator -() const {
    vec<N> ret;

    for ( std::size_t k{0}; k < N; ++k ) {
        ret.entries_[k] = -entries_[k];
    }

    return ret;
}

template <std::size_t N>
double vec<N>::operator * ( vec const &v ) const {
    double ret{0.0};

    for ( std::size_t k{0}; k < N; ++k ) {
        ret += entries_[k] * v.entries_[k];
    }

    return ret;
}

template <std::size_t N>
double &vec<N>::operator []( std::size_t k ) {
    return entries_[k];
}

template <std::size_t N>
double vec<N>::operator []( std::size_t k ) const {
    return entries_[k];
}

template <std::size_t N>
bool vec<N>::operator ==( vec const &v ) const {
    for ( std::size_t k{0}; k < N; ++k ) {
        if ( v.entries_[k] != entries_[k] ) {
            return false;
        }
    }

    return true;
}

template <std::size_t N>
bool vec<N>::operator !=( vec const &v ) const {
    for ( std::size_t k{0}; k < N; ++k ) {
        if ( v.entries_[k] != entries_[k] ) {
            return true;
        }
    }

    return false;
}

/***************************************************************
 * The 1-norm is the sum of the
 * absolute values of the entries:
 *
 *                 N
 *               -----
 *                \        
 *      ||u||  =   )   |u |
 *           1    /      k 
 *               -----
 *               k = 1
 *
 * The 2-norm is the square-root of the
 * sum of the squares of the entries:
 *                        _________
 *                       /  N
 *                      / -----
 *                     /   \      2
 *      ||u||  =      /     )   u 
 *           2       /     /     k 
 *               \  /     -----
 *                \/      k = 1
 *
 *                          
 * The infinity-norm is the maximum entry
 * in absolute value.
 *                         (      )
 *      ||u||   =   max    < |u | >
 *           oo   k = 1..N (  k   )
 *                          
 ***************************************************************/

template <std::size_t N>
double vec<N>::norm( int p ) const {
    switch( p ) {
        case 0: {
            // The infinity-norm or max length
            double ret{ std::abs( entries_[0] ) };

            for ( std::size_t k{1}; k < N; ++k ) {
                double absk{ std::abs( entries_[k] ) };

                if ( absk > ret ) {
                    ret = absk;
                }
            }

            return ret;
        } case 1: {
            // The 1-norm or Manhatten length
            double ret{ std::abs( entries_[0] ) };

            for ( std::size_t k{1}; k < N; ++k ) {
                ret += std::abs( entries_[k] );
            }

            return ret;
        } case 2: {
            // The 2-norm or Euclidean length
            double ret{ entries_[0]*entries_[0] };

            for ( std::size_t k{1}; k < N; ++k ) {
                ret += entries_[k]*entries_[k];
            }

            return std::sqrt( ret );
        } default:
            throw std::invalid_argument( "Supported p-norms include 1, 2 and 0 (infinity)" );
    }
}

template <std::size_t M>
vec<M> operator *( double s, vec<M> const &v ) {
    return v * s;
}

template <std::size_t M>
std::ostream &operator<<( std::ostream &out, vec<M> const &v ) {
    out << "[";

    if ( M != 0 ) {
        out << v[0];

        for ( std::size_t k{1}; k < M; ++k ) {
            out << ", " << v[k];
        }
    }
   
    out << "]";

    return out;
}

template <std::size_t N>
double vec<N>::norm( vec<N> const &v, int p ) {
    return v.norm( p );
}