numeric.h

Go to the documentation of this file.

#ifndef FRECHET_VIEW_NUM_ERROR_H
#define FRECHET_VIEW_NUM_ERROR_H

#include <limits>
#include <math.h>
#include <qdebug.h>
#include <types.h>

#include <boost/config.hpp>
#include <boost/math/special_functions/next.hpp>

#ifdef QUAD_ARITHMETIC
/*  Quadruple precision floats, supported by the quadmath library.
    Depends on platform/compiler:
        Windows+MSC = NO
        Windows+INTEL = YES
        Linux+gcc = probably?
        Linux+INTEL = probably?
        MacOS+clang = NO
        MacOS+gcc = NO (but should?)
        MacOS+INTEL = probably?
*/
# if defined(BOOST_INTEL) || defined(__GNUC__)
#   include <boost/multiprecision/float128.hpp>
    typedef boost::multiprecision::float128 float128; 
#else
    //  fall-back to long double, which in turn may fall back to double
#   pragma message ("no float128 support. falling back to double precision.")
    typedef long double float128;
# endif
#endif

typedef long double accurate_float;


namespace frechet {

using namespace data;

namespace numeric {

static const constexpr double double_eps = std::numeric_limits<double>::epsilon();

template<typename Number>
inline Number sq(const Number& x)   { return x*x; }

template<typename Float>
inline Float sqrt(const Float& x);

template<>
inline double sqrt(const double& x) { return ::sqrt(x); }

template<>
inline long double sqrt(const long double& x) { return ::sqrtl(x); }

template<typename Number>
inline Number abs(const Number& x) { return (x<0) ? -x : x; }

#if defined(BOOST_MP_USE_FLOAT128) || defined(BOOST_MP_USE_QUAD)

template<>
inline float128 sqrt(const float128& x) { return sqrt(x); }

inline const float128& min(const float128& a, const float128& b)
{
    return (a<b) ? a:b;
}
#endif

inline double round_up(double x, int ulps=1) {
    return boost::math::float_advance(x,ulps);
}

template<typename Float>
inline Float euclidean_distance_sq(
        const Float& px, const Float& py,
        const Float& qx, const Float& qy)
{
    return sq<Float>(px-qx) + sq<Float>(py-qy);
}

template<typename Float>
inline Float euclidean_distance(
        const Float& px, const Float& py,
        const Float& qx, const Float& qy)
{
    return sqrt<Float>(euclidean_distance_sq<Float>(px,py,qx,qy));
}

template<typename Float>
inline Float euclidean_distance(const Point& p, const Point& q)
{
    return euclidean_distance<Float>(p.x(),p.y(),q.x(),q.y());
}

template<typename Float>
inline Float euclidean_distance_sq(const Point& p, const Point& q)
{
    return euclidean_distance_sq<Float>(p.x(),p.y(),q.x(),q.y());
}


template<typename Float>
inline Float euclidean_segment_distance(
        const Float& ax, const Float& ay,
        const Float& bx, const Float& by,
        const Float& px, const Float& py)
{
    Float n = (by-ay)*px - by*ax - (bx-ax)*py + bx*ay;
    Float d = euclidean_distance(ax,ay, bx,by);
    Q_ASSERT(d!=0.0);
    return abs(n) / d;
}

template<typename Float>
inline Float euclidean_segment_distance(const QLineF& s, const Point& p)
{
    return euclidean_segment_distance<Float>(
            s.p1().x(), s.p1().y(),
            s.p2().x(), s.p2().y(),
            p.x(),p.y());
}

inline double min(double a, double b)
{
    return (std::isnan(a) || (b < a)) ? b:a;
}

inline double max(double a, double b)
{
    return (std::isnan(a) || (b > a)) ? b:a;
}
inline int lower_floor(double x, int min=0) {
    int i = (int)x;
    if (i==x && i > min) --i;
    return i;
}
inline int compare(double a, double b)
{
    if (a > b) return +1;
    if (a < b) return -1;
    return 0;
}
inline Point normalized(Point p)
{
    double len = sqrt( sq(p.x()) + sq(p.y()) );
    return Point(p.x()/len, p.y()/len);
}
inline Point intersection(Point p1, Point p2, Point p3, Point p4)
{
    double det = (p1.x()-p2.x())*(p3.y()-p4.y()) - (p1.y()-p2.y())*(p3.x()-p4.x());
    double Px = (p1.x()*p2.y()-p1.y()*p2.x())*(p3.x()-p4.x()) - (p1.x()-p2.x())*(p3.x()*p4.y()-p3.y()*p4.x());
    double Py = (p1.x()*p2.y()-p1.y()*p2.x())*(p3.y()-p4.y()) - (p1.y()-p2.y())*(p3.x()*p4.y()-p3.y()*p4.x());

    return Point(Px/det,Py/det);
}
inline Point between(Point p, Point q, double w)
{
    Q_ASSERT(w>=0 && w<=1.0);
    return p*(1-w) + q*w;
}
inline double segment_length(const Curve& c, int i)
{
    return euclidean_distance<double> (c[i],c[i+1]);
}
inline double max_euclidean_distance(const Curve& P, const Curve& Q)
{
    double max=0.0;
    for(int i=0; i < P.size(); ++i)
    {
        const Point& p = P[i];
        for(int j=0; j < Q.size(); ++j)
        {
            double dsq = euclidean_distance_sq<double> (p,Q[j]);
            if (dsq > max) max = dsq;
        }
    }
    return sqrt(max);
}



}}  //  namespace frechet::numeric

#endif //FRECHET_VIEW_NUM_ERROR_H