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

    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 defintions. Ulp calculation.  

Authors
    Tomislav Maric maric@mma.tu-darmstadt.de
\*---------------------------------------------------------------------------*/

#ifndef tolerance_HPP 
#define tolerance_HPP

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

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

namespace geophase { 

    template<typename FT> 
    std::enable_if_t<std::is_floating_point_v<FT>, FT> 
    ulp(FT x)
    {
        // Alternative ulp calculation, costs much more.
        //const FT xabs = std::fabs(x); 
        //return std::nexttoward(xabs, std::numeric_limits<FT>::infinity()) - xabs;
        // TODO: Profile and speed up if possible slows down is_inside_ulp 2x. TM.
        int exp;
        std::frexp(std::fabs(x), &exp);
        return std::ldexp(std::numeric_limits<double>::epsilon(), exp - 1);
    }

    template<typename Vector, std::size_t DIM = 0> 
    enable_if_type_tagged_t<Vector, vector_tag, typename Vector::value_type>
    max_ulp(Vector const& vec)
    {
        auto vecUlp = ulp(vec[DIM]); 

        if constexpr (DIM < (std::size(vec) - 1))
            vecUlp = std::max(vecUlp, ulp(vec[DIM + 1]));

        return vecUlp; 
    }

    template<typename Geo, typename RT=double> 
    struct tolerance
    {
        static constexpr RT value = std::numeric_limits<RT>::epsilon();
    };

    template <typename T> 
    int sign(T val) 
    {
        return (val > T(0)) - (val < T(0));
    }

    template <typename T> 
    int sign(T val, T epsilon) 
    {
        return (val > epsilon) - (val < -epsilon);
    }

} // End namespace geophase

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

#endif

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