#ifndef _MATRIX_h_

#define _MATRIX_h_

#ifdef _MSC_VER
#pragma once
#endif

#include <array>
#include <string>
#include <type_traits>
#include <iostream>

template < typename T, size_t ROWS, size_t COLS >
class Matrix
{
private:
    std::array< T, ROWS* COLS > m_cValues;

public:
    Matrix( void )
    {
        static_assert( std::is_arithmetic< T >::value, "Type T has to be arithmetic!" );
        static_assert( ROWS >= 1 && COLS >= 1, "Minimal row and col dimensions are 1." );
        m_cValues.fill( static_cast<T>( 0 ) );
    }

    ~Matrix( )
    {
    }

    explicit Matrix( std::array< T, ( ROWS* COLS ) > cValues )
        : m_cValues( cValues )
    {
        static_assert( std::is_arithmetic< T >::value, "Type T has to be arithmetic!" );
    }

    template < typename... K >
    Matrix( K ... args )
    {
        static_assert( std::is_arithmetic< T >::value, "Type T has to be arithmetic!" );
        m_cValues = std::array< T, ROWS* COLS > { static_cast<T>( args )... };
    }

    constexpr static auto rows( void ) -> size_t
    {
        return ROWS;
    }

    constexpr static auto cols( void ) -> size_t
    {
        return COLS;
    }

    constexpr static auto size( void ) -> size_t
    {
        return ROWS * COLS;
    }

    constexpr static auto isVector( void ) -> bool
    {
        return ROWS == 1 || COLS == 1;
    }

    auto row( size_t iRow ) const -> Matrix< T, 1, COLS >
    {
        std::array< T, COLS > values;

        for( size_t i = 0; i < COLS; ++i )
            values.at( i ) = at( iRow, i );

        return Matrix< T, 1, COLS >{ values };
    }

    auto col( size_t iCol ) const -> Matrix< T, ROWS, 1 >
    {
        std::array< T, ROWS > values;

        for( size_t i = 0; i < ROWS; ++i )
            values.at( i ) = at( i, iCol );

        return Matrix< T, ROWS, 1 >{ values };
    }

    auto row( size_t iRow, Matrix< T, 1, COLS > mat ) -> void
    {
        for( size_t i = 0; i < COLS; ++i )
            at( iRow, i ) = mat( i );
    }

    auto col( size_t iCol, Matrix< T, ROWS, 1 > mat ) -> void
    {
        for( size_t i = 0; i < ROWS; ++i )
            at( i, iCol ) = mat( i );
    }

    auto at( size_t i ) -> T&
    {
        return m_cValues.at( i );
    }

    auto at( size_t i ) const -> const T&
    {
        return m_cValues.at( i );
    }

    auto at( size_t iRow, size_t iCol ) -> T&
    {
        return m_cValues.at( iRow * COLS + iCol );
    }

    auto at( size_t iRow, size_t iCol ) const -> const T&
    {
        return m_cValues.at( iRow * COLS + iCol );
    }

    auto operator()( size_t i ) -> T&
    {
        return at( i );
    }

    auto operator()( size_t i ) const -> const T&
    {
        return at( i );
    }

    auto operator()( size_t iRow, size_t iCol ) -> T&
    {
        return at( iRow, iCol );
    }

    auto operator()( size_t iRow, size_t iCol ) const -> const T&
    {
        return at( iRow, iCol );
    }

    auto operator+( T other ) const -> Matrix< T, ROWS, COLS >
    {
        auto buf = *this;
        buf += other;
        return buf;
    }

    auto operator-( T other ) const -> Matrix< T, ROWS, COLS >
    {
        auto buf = *this;
        buf -= other;
        return buf;
    }

    auto operator*( T other ) const -> Matrix< T, ROWS, COLS >
    {
        auto buf = *this;
        buf *= other;
        return buf;
    }

    auto operator/( T other ) const -> Matrix< T, ROWS, COLS >
    {
        auto buf = *this;
        buf /= other;
        return buf;
    }

    auto operator+=( T other ) -> Matrix< T, ROWS, COLS >&
    {
        for( auto& val : m_cValues )
            val += other;

        return *this;
    }

    auto operator-=( T other ) -> Matrix< T, ROWS, COLS >&
    {
        for( auto& val : m_cValues )
            val -= other;

        return *this;
    }

    auto operator*=( T other ) -> Matrix< T, ROWS, COLS >&
    {
        for( auto& val : m_cValues )
            val *= other;

        return *this;
    }

    auto operator/=( T other ) -> Matrix< T, ROWS, COLS >&
    {
        for( auto& val : m_cValues )
            val /= other;

        return *this;
    }

    auto operator+( const Matrix< T, ROWS, COLS >& other ) const -> Matrix< T, ROWS, COLS >
    {
        auto buf = *this;

        for( size_t i = 0; i < size( ); ++i )
            buf( i ) += other( i );

        return buf;
    }

    auto operator-( const Matrix< T, ROWS, COLS >& other ) const -> Matrix< T, ROWS, COLS >
    {
        auto buf = *this;

        for( size_t i = 0; i < size( ); ++i )
            buf( i ) -= other( i );

        return buf;
    }

    auto operator+=( const Matrix< T, ROWS, COLS >& other ) -> Matrix< T, ROWS, COLS >&
    {
        for( size_t i = 0; i < size( ); ++i )
            at( i ) += other( i );

        return ( *this );
    }

    auto operator==( const Matrix< T, ROWS, COLS >& other ) const -> bool
    {
        for( size_t i = 0; i < size( ); ++i )
            if( at( i ) != other.at( i ) )
                return false;

        return true;
    }

    auto operator!=( const Matrix< T, ROWS, COLS >& other ) const -> bool
    {
        return !( ( *this ) == other );
    }

    auto operator-=( const Matrix< T, ROWS, COLS >& other ) -> Matrix< T, ROWS, COLS >&
    {
        for( size_t i = 0; i < size( ); ++i )
            at( i ) -= other( i );

        return ( *this );
    }

    auto norm( void ) const -> T
    {
        static_assert( Matrix< T, ROWS, COLS >::isVector( ), "Matrix::norm can only be used on vectors!" );

        T buf = static_cast<T>( 0 );

        for( auto v : m_cValues )
            buf += ( v * v );

        return sqrt( buf );
    }

    auto normalize( void ) -> void
    {
        static_assert( Matrix< T, ROWS, COLS >::isVector( ), "Matrix::normalize can only be used on vectors!" );
        ( *this ) /= this->norm( );
    }

    auto normalized( void ) const -> Matrix< T, ROWS, COLS >
    {
        static_assert( Matrix< T, ROWS, COLS >::isVector( ), "Matrix::normalized can only be used on vectors!" );

        auto buf = *this;
        buf.normalize( );
        return buf;
    }

    friend std::ostream& operator <<( std::ostream& os, const Matrix< T, ROWS, COLS >& v )
    {
        for( size_t i = 0; i < ROWS; ++i )
        {
            for( size_t j = 0; j < COLS; j++ )
                os << v.at( i, j ) << "\t";

            os << "\n";
        }

        return os;
    }

    template < typename K, size_t OROWS, size_t OCOLS >
    auto dot( const Matrix< K, OROWS, OCOLS >& other ) const -> T
    {
        static_assert( Matrix< T, ROWS, COLS >::isVector( ) && Matrix< K, OROWS, OCOLS >::isVector( ), "Matrix::dot can only be used on vectors!" );

        T buf = static_cast<T>( 0 );

        for( size_t i = 0; i < min( size( ), other.size( ) ); ++i )
            buf += at( i ) * other.at( i );

        return buf;
    }

    auto cross( const Matrix< T, ROWS, COLS >& other ) const -> Matrix< T, ROWS, COLS >
    {
        static_assert( Matrix< T, ROWS, COLS >::size( ) == 3, "Matrix::cross can only be called for 3 dimensional vectors!" );

        return{
            at( 1 )* other.at( 2 ) - at( 2 )* other.at( 1 ),
            at( 2 )* other.at( 0 ) - at( 0 )* other.at( 2 ),
            at( 0 )* other.at( 1 ) - at( 1 )* other.at( 0 ),
        };
    }

    auto ncross( const Matrix< T, ROWS, COLS >& other ) const -> Matrix< T, ROWS, COLS >
    {
        static_assert( Matrix< T, ROWS, COLS >::size( ) == 3, "Matrix::ncross can only be called for 3 dimensional vectors!" );

        return cross( other ).normalized( );
    }
};

#endif