/*---------------------------------------------------------------------------*\

    geophase: A C++ template library for geometrical operations on non-convex
    polyhedra with non-planar faces
     
    Copyright (C) 2019 Tomislav Maric, TU Darmstadt 

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

Description
    Arithmetic operation for vectors.    

Author
    Tomislav Maric maric@mma.tu-darmstadt.de

\*---------------------------------------------------------------------------*/

#ifndef VectorArithmetic_HPP 
#define VectorArithmetic_HPP 

#include "tag.hpp"
#include <cstddef>
#include <cmath>

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace geophase { 

// TODO: If the operations on builtins can't be abstracted, write a macro. TM
// TODO: Use rvalue refs, don't overload on them. TM.

template <auto N=0, typename Vector, typename Vector::value_type> 
enable_if_type_tagged_t<void, Vector, vector_tag>
equals(Vector& result, typename Vector::value_type value)
{
    result[N] = value; 
    if constexpr (N + 1 < result.size())
        equals<N + 1>(result, value);
}


template <auto N=0, typename FirstVector, typename SecondVector> 
enable_if_types_tagged_with_tag_t<void, vector_tag, FirstVector, SecondVector>
plus_equals(FirstVector& result, SecondVector const& input)
{
    static_assert(result.size() == input.size(), 
                  "Trying to add vectors with different dimensions."); 
    result[N] += input[N];
    if constexpr (N + 1 < result.size())
        plus_equals<N + 1>(result, input);
}

template <auto N=0, typename FirstVector, typename SecondVector> 
enable_if_types_tagged_with_tag_t<void, vector_tag, FirstVector, SecondVector>
minus_equals(FirstVector& result, SecondVector const& input)
{
    static_assert(result.size() == input.size(), 
                  "Trying to add vectors with different dimensions."); 

    result[N] -= input[N];

    if constexpr (N + 1 < result.size())
        minus_equals<N + 1>(result, input);
}

template <auto N=0, typename Vector, typename AT> 
std::enable_if_t
<   
    type_tagged<Vector , vector_tag>() && 
    std::is_same<typename Vector::value_type, AT>(),
    void 
>
times_equals(Vector& result, AT s)
{
    result[N] *= s;

    if constexpr (N + 1 < result.size())
        times_equals<N + 1>(result, s);
}

template <auto N=0, typename Vector, typename AT> 
std::enable_if_t
<   
    type_tagged<Vector , vector_tag>() && 
    std::is_same<typename Vector::value_type, AT>(),
    void 
>
div_equals(Vector& result, AT s)
{
    result[N] /= s;

    if constexpr (N + 1 < result.size())
        div_equals<N + 1>(result, s);
}

template <auto N=0, typename FirstVector, typename SecondVector> 
enable_if_types_tagged_with_tag_t<FirstVector, vector_tag, FirstVector, SecondVector>
operator+(FirstVector const& first, SecondVector const& second)
{
    FirstVector result {first}; 
    result += second; 
    return result; 
}

template <auto N=0, typename FirstVector, typename SecondVector> 
enable_if_types_tagged_with_tag_t<FirstVector, vector_tag, FirstVector, SecondVector>
operator-(FirstVector const& first, SecondVector const& second)
{
    FirstVector result {first}; 
    result -= second; 
    return result; 
}

template <auto N=0, typename Vector> 
enable_if_type_tagged_t<Vector, vector_tag, Vector>
operator-(Vector const& vec)
{
    return vec * typename Vector::value_type(-1);  
}

template <auto N=0, typename Vector, typename AT> 
std::enable_if_t
<   
    type_tagged<Vector, vector_tag>() &&
    std::is_same<typename Vector::value_type, AT>(), 
    Vector 
>
operator*(Vector const& vec, AT s) 
{
    Vector result {vec}; 
    result *= s; 
    return result; 
}

template <auto N=0, typename Vector, typename AT> 
std::enable_if_t
<   
    type_tagged<Vector, vector_tag>() &&
    std::is_same<typename Vector::value_type, AT>(), 
    Vector 
>
operator*(AT s, Vector const& vec) 
{
    return vec * s; 
}

template <auto N=0, typename Vector, typename AT> 
std::enable_if_t
<   
    type_tagged<Vector, vector_tag>() &&
    std::is_same<typename Vector::value_type, AT>(), 
    Vector 
>
operator/(Vector const& vec, AT s) 
{
    Vector result {vec}; 
    result /= s; 
    return result; 
}

template <auto N=0, typename Vector, typename AT> 
std::enable_if_t
<   
    type_tagged<Vector, vector_tag>() &&
    std::is_same<typename Vector::value_type, AT>(), 
    Vector 
>
operator/(AT s, Vector const& vec) 
{
    return vec / s; 
}

template <auto N=0, typename Vector> 
enable_if_type_tagged_t<Vector, vector_tag, typename Vector::value_type>
sum_sqr(Vector const& vec) 
{
    typename Vector::value_type result {0}; 

    result += vec[N] * vec[N]; 

    if constexpr (N + 1 < vec.size())
        result += sum_sqr<N + 1>(vec);

    return result; 
}

template <auto N=0, typename Vector> 
enable_if_type_tagged_t<Vector, vector_tag, typename Vector::value_type>
mag(Vector const& vec) 
{
    return std::sqrt(sum_sqr(vec)); 
}

template <auto N=0, typename Vector> 
enable_if_type_tagged_t<Vector, vector_tag, typename Vector::value_type>
dot(Vector const& vec1, Vector const& vec2) 
{
    typename Vector::value_type result {0}; 

    static_assert(std::size(vec1) == std::size(vec2)); 

    result += vec1[N] * vec2[N]; 

    if constexpr (N + 1 < vec1.size())
        result += dot<N + 1>(vec1, vec2);

    return result; 
}

template<std::size_t DIM, typename VectorType> 
struct cross_dispatch
{
    static VectorType apply(VectorType const& v1, VectorType const& v2)
    {
        static_assert(std::size(v1) < 3, "Cross product not implemented for 2D and 1D vectors."); 
        static_assert(std::size(v1) > 3, "Cross product not implemented for ND vectors."); 
        return VectorType();
    }
};

template<typename VectorType> 
struct cross_dispatch<3, VectorType>
{
    static VectorType apply(VectorType const& v1, VectorType const& v2)
    {
        return VectorType(v1[1]*v2[2] - v1[2]*v2[1], 
                          v1[2]*v2[0] - v1[0]*v2[2], 
                          v1[0]*v2[1] - v1[1]*v2[0]);

    }
};

template <typename VectorType> 
enable_if_type_tagged_t<VectorType, vector_tag, VectorType>
cross(VectorType const& vec1, VectorType const& vec2) 
{
    static_assert(std::size(vec1) == std::size(vec2)); 
    return cross_dispatch<std::size(vec1), VectorType>::apply(vec1, vec2);  
}

template<typename Vector, auto N = 0>
enable_if_type_tagged_t<Vector, vector_tag>
assign_less(Vector& result, Vector const& input)
{
    static_assert(result.size() == input.size()); 

    if constexpr (result.size() > 0)
    {
        if (input[N] < result[N])
            result[N] = input[N];

        if constexpr (N + 1 < result.size()) 
            assign_less<Vector, N + 1>(result, input); 
    }
}


template<typename Vector, auto N = 0>
enable_if_type_tagged_t<Vector, vector_tag>
assign_greater(Vector& result, Vector const& input)
{
    static_assert(result.size() == input.size()); 

    if constexpr (result.size() > 0)
    {
        if (input[N] > result[N])
            result[N] = input[N];

        if constexpr (N + 1 < result.size()) 
            assign_greater<Vector, N + 1>(result, input); 
    }
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace geophase 

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
