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

    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
    Tolerance or ulp-based equality predicates. 

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

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

#ifndef Equal_HPP 
#define Equal_HPP 

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

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

namespace geophase { 

template<typename FT> 
std::enable_if_t<std::is_floating_point_v<FT>, bool>
equal_by_ulp(FT f1, FT f2, const std::size_t N)
{
    return std::abs(f1 - f2) <= (2.0 * N * std::max(ulp(f1), ulp(f2))); 
}

template<typename Vector, std::size_t D = 0> 
std::enable_if_t<type_tagged<Vector,vector_tag>(), bool>
equal_by_ulp(Vector const& v1, Vector const& v2, const std::size_t N)
{
    bool ulpEqualCmpts = equal_by_ulp(v1[D], v2[D], N); 

    if constexpr (D < (std::size(v1) - 1))
        ulpEqualCmpts = ulpEqualCmpts && equal_by_ulp<Vector, D+1>(v1, v2, N); 

    return ulpEqualCmpts;
}

template<typename FT> 
std::enable_if_t<std::is_floating_point_v<FT>, bool>
equal_by_tolerance(FT f1, FT f2, FT tol)
{
    return std::abs(f1 - f2) <= 
        tol * std::max(std::max(fabs(f1), fabs(f2)), 1.0);
}

template<typename Vector, std::size_t D = 0, typename FT=double> 
std::enable_if_t<type_tagged<Vector,vector_tag>(), bool>
equal_by_tolerance(Vector const& v1, Vector const& v2, FT tol) 
{
    bool tolEqualCmpts = equal_by_tolerance(v1[D], v2[D], tol); 

    if constexpr (D < (std::size(v1) - 1))
        tolEqualCmpts = tolEqualCmpts && equal_by_tolerance<Vector, D+1>(v1, v2, tol); 

    return tolEqualCmpts;
}

template<typename FT=double>
struct EqualByTolerance
{
    FT tolerance_; 

    EqualByTolerance(FT tolerance)
    :
        tolerance_(tolerance)
    {}

    template<typename T>
    bool operator()(T const& t1, T const& t2) 
    {
        return equal_by_tolerance(t1, t2, tolerance_);
    }
};

} // End namespace geophase 

#endif

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