/* ============================================================================
 * I B E X - Interval definition
 * ============================================================================
 * Copyright   : Ecole des Mines de Nantes (FRANCE)
 * License     : This program can be distributed under the terms of the GNU LGPL.
 *               See the file COPYING.LESSER.
 *
 * Author(s)   : Gilles Chabert
 * Bug fixes   : Gilles Trombettoni, Bertrand Neveu
 * Created     : Dec 05, 2011
 * ---------------------------------------------------------------------------- */

#ifndef _IBEX_INTERVAL_H_
#define _IBEX_INTERVAL_H_

#include "ibex_Exception.h"

/* ========================================================*/
/* The following header file is automatically generated by
 * the compilation. It only contains the definition of
 * _IBEX_WITH_GAOL_ or _IBEX_WITH_BIAS_ */
#include "ibex_Setting.h"
/* ======================================================= */




#ifdef _IBEX_WITH_GAOL_
	extern "C" {
	  //#include "gdtoaimp.h"
	}
	#include "gaol/gaol.h"
	#include "gaol/gaol_interval.h"
	/** \brief NEG_INFINITY: double representation of -oo */
	#define NEG_INFINITY (-GAOL_INFINITY)
	/** \brief POS_INFINITY: double representation of +oo */
	#define POS_INFINITY GAOL_INFINITY
#else
#ifdef _IBEX_WITH_BIAS_
	#include "Interval.h"
	/** \brief NEG_INFINITY: double representation of -oo */
	#define NEG_INFINITY BiasNegInf
	/** \brief POS_INFINITY: double representation of +oo */
	#define POS_INFINITY BiasPosInf
#else
#ifdef _IBEX_WITH_FILIB_
	#include "interval/interval.hpp"
//	#include <fp_traits/fp_traits.hpp>
	/* simplify instantiation */
	#define FI_BASE double
	//#define FI_ROUNDING filib::no_rounding
	//#define FI_ROUNDING filib::multiplicative
	#define FI_ROUNDING filib::native_switched
	//#define FI_MODE filib::i_mode_extended
	#define FI_MODE filib::i_mode_extended_flag
	//#define FI_MODE filib::i_mode_normal
	/** \brief NEG_INFINITY: <double> representation of -oo */
	#define NEG_INFINITY filib::primitive::compose(1,0x7FF,0,0)
//	#define NEG_INFINITY filib::primitive::compose(1,0x7FE,(1 << 21)-1,0xffffffff)
	/** \brief POS_INFINITY: <double> representation of +oo */
	#define POS_INFINITY filib::primitive::compose(0,0x7FF,0,0)
//	#define POS_INFINITY filib::primitive::compose(0,0x7FE,(1 << 21)-1,0xffffffff)
	/** \brief IBEX_NAN: <double> representation of NaN */
	#define IBEX_NAN filib::primitive::compose(0,0x7FF,1 << 19,0)
#endif
#endif
#endif

namespace ibex {

template <class T>  class Affine2Main;

class IntervalVector;
class IntervalMatrix;
class IntervalMatrixArray;
class ExprConstant;

/** \defgroup arithmetic Interval Arithmetic */

/** \ingroup arithmetic */
/*@{*/

/**
 * \brief Sets the rounding direction mode of the FPU towards -oo.
 */
void fpu_round_down();

/**
 * \brief Sets the rounding direction mode of the FPU towards +oo.
 */
void fpu_round_up();

/**
 * \brief Sets the rounding direction mode of the nearest
 */
void fpu_round_near();

/*
 * Sets the rounding direction mode of the FPU towards zero.
 */
void fpu_round_zero();

/**
 * \brief Return the previous float
 */
double previous_float(double x);

/**
 * \brief Return the next float
 */
double next_float(double x);


/*@}*/

/**
 * \ingroup arithmetic
 *
 * \brief Interval
 *
 * This class defines the interval interface of IBEX and encapsulates an interval "itv" whose
 * type depends on the chosen implementation (currently: Gaol, Bias of filib).
 *
 * Note that some functions of the Gaol interval interface do not appear here (like "possibly relations")
 * because there are not used by ibex; while other have been introduced (like "ratio_delta"). Some
 * functions are also renamed to match more conventional use.
 *
 * Note that with filib several precision and mode are available. We choose :
 * base type = double
 * rounding_strategy = native_switched
 * interval_mode = i_mode_extended_flag
 *
 */
class Interval {
  public:
	/** \brief Create (-oo,+oo). */
    Interval();

    /** \brief Create [a,b].  */
    Interval(double a, double b);

    /** \brief Create [a,a]. */
    Interval(double a);

    /** \brief True iff *this and x are exactly the same intervals. */
    bool operator==(const Interval& x) const;

    /** \brief True iff *this and x are not exactly the same intervals. */
    bool operator!=(const Interval& x) const;

    /** \brief Set this interval to the empty set. */
    void set_empty();

    /** \brief Set *this to x.
     */
    Interval& operator=(const Interval& x);

    /** \brief Set *this to x.
     */
    template<class T>
    Interval& operator=(const Affine2Main<T>& x);

    /** \brief Set *this to d.
     */
    Interval& operator=(double x);

    /** \brief Intersection of *this and x.
     * \param x - the interval to compute the intersection with.*/
    Interval& operator&=(const Interval& x);

    /** \brief Intersection of *this and x.
     * \param x - the affine form to compute the intersection with.*/
    template<class T>
    Interval& operator&=( const Affine2Main<T>& x);

    /** \brief Union of *this and x.
     * \param x - the interval to compute the hull with.*/
    Interval& operator|=(const Interval& x);

    /** \brief Union of *this and x.
     * \param x - the affine form to compute the hull with.*/
    template<class T>
    Interval& operator|=(const Affine2Main<T>& x);

    /**
     * \brief Add [-rad,+rad] to *this.
     *
     * Return a reference to *this.
     */
    Interval& inflate(double rad);

    /** \brief Lower bound.
     *
     * Return the lower bound of *this. */
    double lb() const;

    /** \brief Upper bound.
     *
     * Return the upper bound of *this. */
    double ub() const;

    /** \brief Midpoint.
     *
     * Returns the midpoint of *this.
     * The return point is guaranteed to be included in *this
     * but not necessarily to be the closest floating point
     * from the real midpoint.
     *
     * Cases are:
     * - \emptyset  -> Quiet NaN
     * - [-oo, +oo] -> midP = 0.0
     * - [-oo, b]   -> midP = -MAXREAL
     * - [a, +oo]   -> midP = MAXREAL
     * - [a, b]     -> midP ~ a + .5*(b-a) */
    double mid() const;

    /**
     * \brief Radius.
     *
     * Return the diameter of *this.
     * By convention, 0 if *this is empty.*/
    double rad() const;

    /**
     * \brief Diameter.
     *
     * Return the diameter of *this.
     * By convention, 0 if *this is empty.*/
    double diam() const;

    /**
     * \brief Mignitude.
	 *
     * Returns the mignitude of *this:
     * <lu>
     * <li> +(lower bound)  if *this > 0
     * <li> -(upper bound) if *this < 0
     * <li> 0 otherwise.
     * </lu> */
    double mig() const;

    /**
     * \brief Magnitude.
	 *
     * Returns the magnitude of *this:
     * mag(*this)=max(|lower bound|, |upper bound|). */
    double mag() const;

    /**
     * \brief True iff this interval is a subset of \a x.
     *
     * \note Always return true if *this is empty.
     */
    bool is_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is a subset of \a x and not \a x itself.
     *
     * \note In particular, (-oo,oo) is not a strict subset of (-oo,oo)
     * and the empty set is not a strict subset of the empty set although
     * in both cases, the first is inside the interior of the second.
     */
    bool is_strict_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is in the interior of \a x.
     *
    * \note In particular, (-oo,oo) is in the interior of (-oo,oo)
     * and the empty set is in the interior of the empty set.
     * \note Always return true if *this is empty.
     */
    bool is_interior_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is in the interior of \a x and different from x.
     *
     * \note In particular, (-oo,oo) is not "strictly" in the interior of (-oo,oo)
     * and the empty set is not "strictly" in the interior of the empty set.
     */
    bool is_strict_interior_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is a superset of \a x.
     *
     * \note Always return true if x is empty.
     */
    bool is_superset(const Interval& x) const;

    /**
     * \brief True iff this interval is a superset of \a x different from x.
     *
     * \see #is_strict_subset(const Interval&) const.
     */
    bool is_strict_superset(const Interval& x) const;

    /**
     * \brief True iff *this contains \a d.
     *
     * \note d can also be an "open bound", i.e., infinity.
     * So this function is not restricted to a set-membership
     * interpretation.
     */
    bool contains(const double& d) const;

    /**
     * \brief True iff the interior of *this contains \a d.
     *
     */
    bool interior_contains(const double& d) const;

    /**
     * \brief True iff *this and \a x intersect.
     */
    bool intersects(const Interval &x) const;

    /**
     * \brief True iff *this and \a x intersect and the intersection has a non-null volume.
     *
     * Equivalently, some interior points (of this or x) must belong to the intersection.
     */
    bool overlaps(const Interval &x) const;

    /**
     * \brief True iff *this and \a x do not intersect.
     *
     */
    bool is_disjoint(const Interval &x) const;

    /**
     * \brief True iff *this is empty.
     */
    bool is_empty() const;

    /**
     * \brief True iff *this is degenerated.
     *
     * An interval is degenerated if it is of the form [a, a]
     *
     * \note An empty interval is considered here as degenerated. */
    bool is_degenerated() const;

    /**
     * \brief True if one bound of *this is infinite.
     *
     * \note An empty interval is always bounded.
     */
    bool is_unbounded() const;

    /**
     * \brief True iff *this can be bisected into two non-degenerated intervals.
     *
     * Examples of non bisectable intervals are [0,next_float(0)] or [DBL_MAX,+oo).
     */
    bool is_bisectable() const;

    /**
     * \brief Relative Hausdorff distance between *this and x.
     *
     * The relative distance is basically distance(x)/diam(*this).
     * \see #ibex::distance (const ibex::Interval &x1, const ibex::Interval &x2).
     */
    double rel_distance(const Interval& x) const;

    /*
     * \brief The complementary of x.
     */
    int complementary(Interval& c1, Interval& c2) const;

    /**
     * \brief x\y
     */
    int diff(const Interval& y, Interval& c1, Interval& c2) const;

    /** \brief Return -*this. */
    Interval operator-() const;

    /** \brief Add \a d to *this and return the result.  */
    Interval& operator+=(double d);

    /** \brief Subtract \a d to *this and return the result. */
    Interval& operator-=(double d);

    /** \brief Multiply *this by \a d and return the result. */
    Interval& operator*=(double d);

    /** \brief Divide *this by \a d and return the result. */
    Interval& operator/=(double d);

    /** \brief Add \a x to *this and return the result. */
    Interval& operator+=(const Interval& x);

    /** \brief Subtract \a x to *this and return the result. */
    Interval& operator-=(const Interval& x);

    /** \brief Multiply *this by \a x and return the result. */
    Interval& operator*=(const Interval& x);

    /**
     * \brief Divide *this by \a x and return the result.
     *
     * Does better than *this=*this/x: because calculates
     * the union of *this/x as intermediate result. */
    Interval& operator/=(const Interval& x);

    /**
     * \brief Set this interval to the intersection of itself with the division of two others.
     *
     * \param x - the numerator
     * \param y - the divisor
     * \param out2 - In return, *this and out2 contains the lower and upper part respectively
     * of the division. If the result of the generalized division and intersection
     * is a single interval, out2 is set to the empty interval.
     * \code
     * Interval intv(-10,10);
     * Interval out2;
     * intv.div2_inter(Interval(2,3), Interval(-1,2), out2);
     * cout << intv << " " << out2 << endl;  // will display: [-10,-2] [1,10]
     * \endcode
	 *
     * \return \c true if the intersection is non empty.
     * \note Contrary to the "cset" theory, the result is empty if \a y=[0,0] (whatever \a x is).
     */
    bool div2_inter(const Interval& x, const Interval& y, Interval& out2);

    /**
     * \brief Set this interval to the intersection of itself with the division of two others.
     */
    Interval& div2_inter(const Interval& x, const Interval& y);

    /**
     * \brief Return diam(*this)-diam(x), for x\subseteq *this [deprecated]
     *
     * Deprecated. Kept for compatibility with ibex 1.xx.
     *
     * \pre \a x must be included in this interval.
     * \note The result may be +oo (if the set difference is infinite).
     * \note An empty interval is considered here to have a null diamater (as a degenerated interval). <br>
     * If either \a x or this interval is empty, then the method returns the diameter of this interval
     * (which is 0 if the latter is empty).
     */
    double delta(const Interval& x) const;

    /**
     * \brief Compute the ratio of the diameter to #delta(x) [deprecated].
     *
     * Deprecated. Kept for compatibility with ibex 1.xx.
     *
     * \pre \a x must be included in this interval.
     * \note An empty interval is considered to have a null diamater (as a degenerated interval). <br>
     * <ul><li>If either \a x or this interval is empty, then
     * <ul><li>the method returns 1 (100% of reduction) if this diameter is not null,
     *     <li>0 otherwise (as if 0/0=0).</ul>
     * <li>As a pure convention, the method returns \c 1 if one bound of this interval is infinite and the corresponding bound of \a x
     * is not (in particular if this interval is unbounded and \a x not). </ul>
     */
    double ratiodelta(const Interval& x) const;

    /**
     * \brief Bisect *this into two subintervals.
     *
     * \param ratio - says where to split (0.5=middle)
     * \pre is_bisectable() must be true.
     * \pre 0<ratio<1.
     */
    std::pair<Interval,Interval> bisect(double ratio=0.5) const;

    /** \brief pi. */
    static const Interval PI;
    /** \brief 2*pi. */
    static const Interval TWO_PI;
    /** \brief pi/2. */
    static const Interval HALF_PI;
    /** \brief the empty interval. */
    static const Interval EMPTY_SET;
    /** \brief (-oo,oo). */
    static const Interval ALL_REALS;
    /** \brief [0,0]. */
    static const Interval ZERO;
    /** \brief [1,1]. */
    static const Interval ONE;
    /** \brief [0,+oo). */
    static const Interval POS_REALS;
    /** \brief (-oo,0]. */
    static const Interval NEG_REALS;

//    friend class IntervalVector;

    typedef Interval SCALAR;
    typedef IntervalVector VECTOR;
    typedef IntervalMatrix MATRIX;
    typedef IntervalMatrixArray MATRIX_ARRAY;

    /**
     * \brief Cast the interval to an expression
     */
    operator const ExprConstant&() const;

    //private:
#ifdef _IBEX_WITH_GAOL_
	/* \brief Wrap the gaol-interval [x]. */
    Interval(const gaol::interval& x);
    /* \brief Assign this to the gaol-interval [x]. */
    Interval& operator=(const gaol::interval& x);
	gaol::interval itv;
#else
#ifdef _IBEX_WITH_BIAS_
	/* \brief Wrap the bias-interval [x]. */
    Interval(const INTERVAL& x);
    /* \brief Assign this to the bias-interval [x]. */
    Interval& operator=(const INTERVAL& x);

    INTERVAL itv;
#else
#ifdef _IBEX_WITH_FILIB_
    //typedef filib::interval<double, filib::native_switched,filib::i_mode_extended_flag> INTERVAL;
    typedef filib::interval<FI_BASE,FI_ROUNDING,FI_MODE> FI_INTERVAL;
	/* \brief Wrap the filib-interval [x]. */
    Interval(const FI_INTERVAL& x);
    /* \brief Assign this to the filib-interval [x]. */
    Interval& operator=(const FI_INTERVAL& x);

    FI_INTERVAL itv;

#endif
#endif
#endif
};

/** \ingroup arithmetic */
/*@{*/

/** \brief Stream out \a x. */
std::ostream& operator<<(std::ostream& os, const Interval& x);

/** \brief \f$[x]_1\cap [x]_2\f$.
 * \return Interval::EMPTY if the intersection is empty. */
Interval operator&(const Interval& x1, const Interval& x2);

/** \brief \f$\square([x]_1\cup [x]_2)\f$. */
Interval operator|(const Interval& x1, const Interval& x2);

/** \brief [x]+d. */
Interval operator+(const Interval& x, double d);

/** \brief [x]-d. */
Interval operator-(const Interval& x, double d);

/** \brief [x]*d. */
Interval operator*(const Interval& x, double d);

/** \brief [x]/d. */
Interval operator/(const Interval& x, double d);

/** \brief Generalized division
 *
 * The result of x/y is stored in out1 and out2. These intervals contain the lower and upper part respectively
 * of the division. If the result of the generalized division is a single interval, out2
 * is set to the empty interval. If the result is the empty set, both out1 and out2 are set to the
 * empty interval.
 *
 * \note Contrary to the "cset" theory, the result is empty if \a y=[0,0] (whatever \a x is).
 */
void div2(const Interval& x, const Interval& y, Interval& out1, Interval& out2);

/** \brief d+[x]. */
Interval operator+(double d,const Interval& x);

/** \brief d-[x]. */
Interval operator-(double d, const Interval& x);

/** \brief d*[x]. */
Interval operator*(double d, const Interval& x);

/** \brief d/[x]. */
Interval operator/(double d, const Interval& x);

/** \brief [x]_1+[x]_2. */
Interval operator+(const Interval& x1, const Interval& x2);

/** \brief [x]_1-[x]_2. */
Interval operator-(const Interval& x1, const Interval& x2);

/** \brief [x]_1*[x]_2. */
Interval operator*(const Interval& x1, const Interval& x2);

/** \brief [x]_1/[x]_2. */
Interval operator/(const Interval& x1, const Interval& x2);

/** \brief Hausdorff distance of [x]_1 and [x]_2. */
double distance(const Interval &x1, const Interval &x2);

/** \brief [x]^2 */
Interval sqr(const Interval& x);

/** \brief sqrt{[x]} */
Interval sqrt(const Interval& x);

/** \brief [x]^n. */
Interval pow(const Interval& x, int n);

/** \brief [x]^d. */
Interval pow(const Interval& x, double d);

/** \brief [x]^[y]. */
Interval pow(const Interval &x, const Interval &y);

/** \brief n^{th} root of [x]. */
Interval root(const Interval& x, int n);

/** \brief exp([x]). */
Interval exp(const Interval& x);

/** \brief log([x]). */
Interval log(const Interval& x);

/** \brief cos([x]). */
Interval cos(const Interval& x);

/** \brief sin([x]). */
Interval sin(const Interval& x);

/** \brief tan([x]). */
Interval tan(const Interval& x);

/** \brief acos([x]). */
Interval acos(const Interval& x);

/** \brief asin([x]). */
Interval asin(const Interval& x);

/** \brief atan([x]). */
Interval atan(const Interval& x);

/** \brief atan2([x],[y]). */
Interval atan2(const Interval& x, const Interval& y);

/** \brief cosh([x]). */
Interval cosh(const Interval& x);

/** \brief sinh([x]). */
Interval sinh(const Interval& x);

/** \brief tanh([x]). */
Interval tanh(const Interval& x);

/** \brief acosh([x]). */
Interval acosh(const Interval& x);

/** \brief asinh([x]). */
Interval asinh(const Interval& x);

/** \brief atanh([x]). */
Interval atanh(const Interval& x);

/** \brief \f$abs([x]) = \{|x|, x\in[x]\}.\f$. */
Interval abs(const Interval &x);

/** \brief Maximum of two intervals.
 *
 * Return \f$\max([x],[y]) = \{ \max\{x,y\}, x\in[x], y\in[y] \}\f$. */
Interval max(const Interval& x, const Interval& y);

/** \brief Minimum of two intervals.
 *
 *  Return \f$\min([x],[y]) = \{ \min\{x,y\}, x\in[x], y\in[y] \}\f$. */
Interval min(const Interval& x, const Interval& y);

/** \brief Sign of [x].
 *
 *  Return \f$sign([x]) = hull \{ sign{x}, x\in[x] \}\f$.
 * \remark By convention, \f$ 0\in[x] \Longrightarrow sign[x]=[-1,1]\f$. */
Interval sign(const Interval& x);

/** \brief Chi of [a], [b] and [c].
 *
 *  Return \f$chi([a],[b],[c]) = 0.5*(1-\sign([a]))*[b] + 0.5*(\sign([a])+1)*[c]. \f$
 * \remark  chi([a],[b],[c]) =[b] if [a]<=0, [c] if [a]>0, hull \{[b], [c]\} else.  */
Interval chi(const Interval& a, const Interval& b, const Interval& c);

/**
 * \brief Return the largest integer interval included in x.
 */
Interval integer(const Interval& x);

/** \brief Projection of y=x_1+x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1+x_2\}\f$. */
bool bwd_add(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x_1-x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1-x_2\}\f$. */
bool bwd_sub(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x_1*x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1\times x_2\}\f$. */
bool bwd_mul(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x_1/x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1/x_2\}\f$. */
bool bwd_div(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x^2.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=x^2 \}\f$. */
bool bwd_sqr(const Interval& y, Interval& x);

/** \brief Projection of y=sqrt{x}.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=sqrt{x} \}\f$. */
bool bwd_sqrt(const Interval& y, Interval& x);

/** \brief Projection of y=x^n.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=x^n \}\f$. */
bool bwd_pow(const Interval& y, int n, Interval& x);

/** \brief Projection of y=x_1^{x_2}.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1^{x_2}\}\f$. */
bool bwd_pow(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of the \f$y=x^{\frac{1}{n}}\f$.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=x^{\frac{1}{n}} \}\f$. */
bool bwd_root(const Interval& y, int n, Interval& x);

/** \brief Projection of y=exp(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\exp(x) \}\f$. */
bool bwd_exp(const Interval& y,  Interval& x);

/** \brief Projection of y=log(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\log(x) \}\f$. */
bool bwd_log(const Interval& y,  Interval& x);

/** \brief Projection of y=cos(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\cos(x) \}\f$. */
bool bwd_cos(const Interval& y,  Interval& x);

/** \brief Projection of y=sin(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\sin(x) \}\f$. */
bool bwd_sin(const Interval& y,  Interval& x);

/** \brief Projection of y=tan(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\tan(x) \}\f$. */
bool bwd_tan(const Interval& y,  Interval& x);

/** \brief Projection of y=acos(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\arccos(x) \}\f$. */
bool bwd_acos(const Interval& y,  Interval& x);

/** \brief Projection of y=asin(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\arcsin(x) \}\f$. */
bool bwd_asin(const Interval& y,  Interval& x);

/** \brief Projection of y=atan(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\arctan(x) \}\f$. */
bool bwd_atan(const Interval& y,  Interval& x);

/** \brief Projection of y=atan2(x_1,x_2).
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=atan2(x_1,x_2)\f$. */
bool bwd_atan2(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=cosh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\cosh(x) \}\f$. */
bool bwd_cosh(const Interval& y,  Interval& x);

/** \brief Projection of y=sinh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\sinh(x) \}\f$. */
bool bwd_sinh(const Interval& y,  Interval& x);

/** \brief Projection of y=tanh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\tanh(x) \}\f$. */
bool bwd_tanh(const Interval& y,  Interval& x);

/** \brief Projection of y=acosh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=arccosh(x) \}\f$. */
bool bwd_acosh(const Interval& y,  Interval& x);

/** \brief Projection of y=asinh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=arcsinh(x) \}\f$. */
bool bwd_asinh(const Interval& y,  Interval& x);

/** \brief Projection of y=atanh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=arctanh(x) \}\f$. */
bool bwd_atanh(const Interval& y,  Interval& x);

/** \brief Projection of y=|x|.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=|x| \}\f$. */
bool bwd_abs(const Interval& y,  Interval& x);

/** \brief Projection of y=max(x_1,x_2).
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=\max(x_1,x_2)\}\f$. */
bool bwd_max(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=min(x_1,x_2).
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=\min(x_1,x_2)\}\f$. */
bool bwd_min(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=sign(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=sign(x) \}\f$. */
bool bwd_sign(const Interval& y, Interval& x);

/** \brief Projection of f=chi(a,b,c). */
bool bwd_chi(const Interval& f, Interval& a, Interval& b, Interval& c);

/** \brief Contract x w.r.t. the fact that it must be integral.
 *
 */
bool bwd_integer(Interval& x);

} // end namespace ibex

#ifdef _IBEX_WITH_GAOL_
#include "ibex_gaol_Interval.h_"
#else
#ifdef _IBEX_WITH_BIAS_
#include "ibex_bias_Interval.h_"
#else
#ifdef _IBEX_WITH_FILIB_
#include "ibex_filib_Interval.h_"
#endif
#endif
#endif

/*@}*/

/*============================================ features with common implementation ============================================ */

namespace ibex {

inline Interval::Interval() : itv(NEG_INFINITY, POS_INFINITY) {

}

inline Interval::Interval(double a, double b) : itv(a,b) {
	if (a==POS_INFINITY || b==NEG_INFINITY || a>b) *this=EMPTY_SET;
}

inline Interval::Interval(double a) : itv(a,a) {
	if (a==NEG_INFINITY || a==POS_INFINITY) *this=EMPTY_SET;
}

inline bool Interval::operator==(const Interval& x) const {
	return (is_empty() && x.is_empty()) || (lb()==x.lb() && ub()==x.ub());
}

inline Interval& Interval::operator=(double x) {
	if (x==NEG_INFINITY || x==POS_INFINITY)
		*this=EMPTY_SET;
	else
		itv = x;
	return *this;
}


inline Interval& Interval::operator=(const Interval& x) {
	itv = x.itv;
	return *this;
}

inline Interval& Interval::inflate(double radd) {
	(*this) += Interval(-radd,radd);
	return *this;
}

inline bool Interval::operator!=(const Interval& x) const {
	return !(*this==x);
}

inline double Interval::rad() const {
	if (is_empty()) return 0;
	else if (is_unbounded()) return POS_INFINITY;
	else {
		double t = mid();
		double t1 =(t-*this).ub();
		double t2= (*this-t).ub();
		return (t1>t2) ? t1 : t2;
	}
}

inline bool Interval::is_bisectable() const {
	if (is_empty()) return false;
	double m=mid();
	return (lb()<m && m<ub());
}


inline double Interval::rel_distance(const Interval& x) const {
	  double d=distance(*this,x);
	  if (d==POS_INFINITY) return 1;
	  double D=diam();
	  return (D==0 || D==POS_INFINITY) ? 0.0 : (d/D); // if diam(this)=infinity here, necessarily d=0
}

inline double distance(const Interval &x1, const Interval &x2) {

    if (x1.is_empty()) return x2.rad();

    if (x2.is_empty()) return x1.rad();

    if (x1.lb()==NEG_INFINITY) {
    	if (x2.lb()!=NEG_INFINITY)
    		return POS_INFINITY;
    	else if (x1.ub()==POS_INFINITY) {
    		if (x2.ub()==POS_INFINITY) return 0.0;
    		else return POS_INFINITY;
    	}
    	else if (x2.ub()==POS_INFINITY) return POS_INFINITY;
    	else return fabs(x1.ub()-x2.ub());
    }
    else if (x1.ub()==POS_INFINITY) {
    	if (x2.ub()!=POS_INFINITY)
    		return POS_INFINITY;
    	else if (x2.lb()==NEG_INFINITY)
    		return POS_INFINITY;
    	else return fabs(x1.lb()-x2.lb());
    }
    else if (x2.is_unbounded())
    	return POS_INFINITY;
    else
#ifdef _IBEX_WITH_GAOL_
    	return hausdorff(x1,x2);
#else
#ifdef _IBEX_WITH_BIAS_
    	return Distance(x1.itv,x2.itv);
#else
#ifdef _IBEX_WITH_FILIB_
    	return x1.itv.dist(x2.itv);
#endif
#endif
#endif

}

inline Interval sign(const Interval& x) {
	return x.ub()<0 ? Interval(-1,-1) : x.lb()>0 ? Interval(1,1) : Interval(-1,1);
}

inline Interval chi(const Interval& a, const Interval& b, const Interval& c){
	if (a.ub()<=0) {
		return b;
	} else if (a.lb()>0) {
		return c;
	} else {
		return b|c;
	}
}

// ArcTan2 by Jordan Ninin (Sept 2013) car ni Gaol, ni Profil/Bias ni Filib n'implemente atan2
inline Interval atan2(const Interval& y, const Interval& x) {
	if (y.is_empty() || x.is_empty()) return Interval::EMPTY_SET;

	if ( x.lb()> 0) {
		return atan(y/x);
	}
	else if  (x.ub()< 0 && y.lb()>= 0) {
		return atan(y/x)+Interval::PI;
	}
	else if (x.ub() < 0 && y.ub()< 0) {
		return atan(y/x)-Interval::PI;
	}
	else if (x.ub() < 0) {
		return Interval(-1,1)*Interval::PI;
	}
	else if (x.contains(0.0) && y.contains(0.0) ) {
		return Interval(-1,1)*Interval::PI;
	}
	else {
		return 2*atan((sqrt(sqr(x)+sqr(y))-x)/y);
	}
}

/*
// ArcTan2 by Firas Khemane (feb 2010)
inline Interval atan2(const Interval& b, const Interval& a) {
	if (a.is_empty() || b.is_empty()) return Interval::EMPTY_SET;

	if (( a.lb() > 0 && b.lb() > 0) || (  a.lb() > 0 &&   b.ub() < 0 )  || ( b.lb() < 0 && b.ub()> 0  && a.lb() > 0 ) )

		return atan(b/a);

	else if  ( a.lb()<=0  &&  a.ub()>0 )

		return 2*atan(b/a);

	else if  (a.ub()< 0 && b.lb()> 0)

		return atan(b/a)+Interval::PI;

	else if (a.ub() < 0 && b.ub()< 0)

		return atan(b/a)-Interval::PI;

	else if (b.lb() < 0 && b.ub() > 0 && a.ub()< 0 ) {

		if (  (atan(b/a)+Interval::PI).lb() <= (atan(b/a)-Interval::PI).ub() )

			return 2*INTERVAL((atan(b/a)+Interval::PI).lb(), (atan(b/a)-Interval::PI).ub());

		else
			return 2*INTERVAL((atan(b/a)-Interval::PI).lb() , (atan(b/a)+Interval::PI).ub() );
	}
}
*/

inline bool bwd_add(const Interval& y, Interval& x1, Interval& x2) {
	if ((x1 &= y-x2).is_empty()) { x2.set_empty(); return false; }
	if ((x2 &= y-x1).is_empty()) { x1.set_empty(); return false; }
	return true;
}

inline bool bwd_sub(const Interval& y, Interval& x1, Interval& x2) {
	if ((x1 &= y+x2).is_empty()) { x2.set_empty(); return false; }
	if ((x2 &= x1-y).is_empty()) { x1.set_empty(); return false; }
	return true;
}

inline bool bwd_div(const Interval& y, Interval& x1, Interval& x2) {
	if ((x1 &= y*x2).is_empty()) { x2.set_empty(); return false; }
	Interval tmp=y;
	bwd_mul(x1, tmp, x2);
	if (x2.is_empty()) { x1.set_empty(); return false; }
	return true;
}

inline bool bwd_sqrt(const Interval& y, Interval& x) {
	x &= sqr(y);
	return !x.is_empty();
}

inline bool bwd_root(const Interval& y, int n, Interval& x) {
	x &= pow(y,n);
	return !x.is_empty();
}

inline bool bwd_exp(const Interval& y,  Interval& x) {
	x &= log(y);
	return !x.is_empty();
}

inline bool bwd_log(const Interval& y,  Interval& x) {
	x &= exp(y);
	return !x.is_empty();
}

inline bool bwd_acos(const Interval& y,  Interval& x) {
	x &= cos(y);
	return !x.is_empty();
}

inline bool bwd_asin(const Interval& y,  Interval& x) {
	x &= sin(y);
	return !x.is_empty();
}

inline bool bwd_atan(const Interval& y,  Interval& x) {
	x &= tan(y);
	return !x.is_empty();
}

inline bool bwd_acosh(const Interval& y,  Interval& x) {
	if (y.ub()<0.0) {
		x.set_empty(); return false;
	}
	else {
		x &= cosh(Interval(y.lb()<0?0:y.lb(),y.ub()));
		return !x.is_empty();
	}
}

inline bool bwd_asinh(const Interval& y,  Interval& x) {
	x &= sinh(y);
	return !x.is_empty();
}

inline bool bwd_atanh(const Interval& y,  Interval& x) {
	x &= tanh(y);
	return !x.is_empty();
}

inline bool bwd_max(const Interval& y, Interval& x1, Interval& x2) {

	/* ---- Disjoint intervals ---- */
	if (x2.lb()>x1.ub() || y.lb()>x1.ub()) {
		/* then, max(x,x2) is necessarily x2 */
		x2 &= y; return !x2.is_empty();
	} else if (x1.lb()>x2.ub() || y.lb()>x2.ub()) {
		x1 &= y; return !x1.is_empty();
	}
	/*------------------------------*/

	if (y.ub()<x1.lb() || y.ub()<x2.lb()) {
		return false; // inconsistency
	}

	/* At this point, x, y and y all mutually intersect. */
	if (x1.ub()>y.ub())
		x1=Interval(x1.lb(),y.ub());
	if (x2.ub()>y.ub())
		x2=Interval(x2.lb(),y.ub());

	return true;
}

inline bool bwd_min(const Interval& y, Interval& x1, Interval& x2) {

	Interval mx1=-x1;
	Interval mx2=-x2;

	if (!bwd_max(-y,mx1,mx2)) return false;

	x1=-mx1;
	x2=-mx2;
	return true;
}

inline bool bwd_sign(const Interval& y,  Interval& x) {
	if(y.lb()>0)
		x &= Interval::POS_REALS;
	else if(y.ub()<0)
		x &= Interval::NEG_REALS;
	return !x.is_empty();
}


inline bool bwd_atan2(const Interval& theta, Interval& y, Interval& x) {
	bool b=true;
	if (theta.is_empty()) {
		x.set_empty(); y.set_empty();
		b=false;
	}
	not_implemented("bwd_atan2 non implemented yet");
	//  <=> cos(Theta)*X - sin(Theta)*Y > 0;
	//      sin(Theta)*X + cos(Theta)*Y = 0
/*
	Interval costheta, sintheta, xcostheta, xsintheta, ycostheta, ysintheta, equ1, equ2;
	costheta = cos(theta);
	sintheta = sin(theta);
	for (int i=0; b && (i<1);i++) {
		xcostheta = x*costheta;
		xsintheta = x*sintheta;
		ycostheta = y*costheta;
		ysintheta = y*sintheta;
		equ1 = (xcostheta - ysintheta)& Interval::POS_REALS;
		equ2 = (xsintheta + ycostheta)& Interval::ZERO;
		b= b && (!equ1.is_empty()||!equ2.is_empty());
		b= b && bwd_sub(equ1, xcostheta, ysintheta);
		b= b && bwd_add(equ2, xsintheta, ycostheta);
		b= b && bwd_mul(xcostheta, x, costheta);
		b= b && bwd_mul(xsintheta, x, sintheta);
		b= b && bwd_mul(ysintheta, y, sintheta);
		b= b && bwd_mul(ycostheta, y, costheta);
	}
*/
	return b;
}


inline bool bwd_chi(const Interval& f, Interval& a, Interval& b, Interval& c){
	if      (a.ub()<=0) {if ((b &= f).is_empty()) { a.set_empty(); c.set_empty(); return false; } }
	else if (a.lb()>0)  {if ((c &= f).is_empty()) { a.set_empty(); b.set_empty(); return false; } }

	if (f.is_disjoint(b)) {
		if ((a &= Interval::POS_REALS).is_empty()) { b.set_empty(); c.set_empty(); return false; }
		if ((c &= f).is_empty()) { a.set_empty(); b.set_empty(); return false; }
	}
	if (f.is_disjoint(c)) {
		if ((a &= Interval::NEG_REALS).is_empty()) { b.set_empty(); c.set_empty(); return false; }
		if ((b &= f).is_empty()) { a.set_empty(); c.set_empty(); return false; }
	}
	return true;
}


inline bool bwd_integer(Interval& x) {
	return !(x = integer(x)).is_empty();
}

} // end namespace ibex

#endif // _IBEX_INTERVAL_H_
