/Users/mourrain/Devel/mmx/numerix/include/numerix/affine.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : affine.hpp
* DESCRIPTION: Affine numbers
* COPYRIGHT  : (C) 2011  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_AFFINE_HPP
#define __MMX_AFFINE_HPP
#include <basix/function.hpp>
#include <numerix/ball.hpp>
namespace mmx {
#define TMPL_DEF template<typename C, typename VC= C>
#define TMPL template<typename C, typename VC>
#define R Radius_type(C)
#define VR typename affine_limits_helper<C,VC>::VT
#define Affine affine<C,VC>
#define Affine_rep affine_rep<C,VC>
#define Abs_affine affine<Abs_type(C),Abs_type(VC) >
#define Real_affine affine<Real_type(C),Real_type(VC) >
TMPL class affine;
TMPL struct affine_limits_helper;
TMPL inline C  base   (const Affine& z);
TMPL inline VC slope  (const Affine& z);
TMPL inline VR limits (const Affine& z);
TMPL inline Affine operator - (const Affine& z1, const Affine& z2);

/******************************************************************************
* Routines for domains of affine numbers
******************************************************************************/

template<typename C> inline C dot (const C& x, const C& y) { return x * y; }

TMPL
struct affine_limits_helper {
  typedef Radius_type(VC) VT;
  static inline VT infinite () {
    return Maximal (VT); }
  static inline VC get_linear (const VC& lin, const VC& dom) {
    return lin + sharpen (dom); }
  static inline VT get_limits (const VC& dom) {
    return radius (dom); }
  static inline VC get_domain (const VT& lim) {
    return blur (VC (0), lim); }
  static inline C truncate (const VC& lin, const VT& lim) {
    return dot (lin, get_domain (lim)); }
  static inline VT common (const VT& lim1, const VT& lim2) {
    if (hard_eq (lim1, lim2)) return lim1;
    return inf (radius (lim1), radius (lim2)); }
};

/******************************************************************************
* The affine type
******************************************************************************/

TMPL_DEF
class affine {
MMX_ALLOCATORS
private:
  C  cst;
  VC lin;
  VR lim;

public:
  typedef affine_limits_helper<C,VC> Helper;
  inline affine ():
    cst (0), lin (0), lim (Helper::infinite ()) {}
  inline affine (const format<C>& fm1, const format<VC>& fm2):
    cst (promote (0, fm1)), lin (promote (0, fm2)),
    lim (Helper::infinite ()) {}
  template<typename T> affine (const T& cst2):
    cst (cst2), lin (0), lim (Helper::infinite ()) {}
  inline affine (const C& z):
    cst (copy (z)), lin (0), lim (Helper::infinite ()) {}
  inline affine (const C& cst2, const VC& lin2):
    cst (cst2), lin (lin2), lim (Helper::infinite ()) {}
  inline affine (const C& cst2, const VC& lin2, const VR& lim2):
    cst (cst2), lin (lin2), lim (lim2) {}
  template<typename CT,typename VCT> inline affine (const affine<CT,VCT>& z):
    cst (as<C> (base (z))),
    lin (as<VC> (slope (z))),
    lim (as<VR > (limits (z))) {}
  friend C base LESSGTR (const Affine& z);
  friend VC slope LESSGTR (const Affine& z);
  friend VR limits LESSGTR (const Affine& z);
};
DEFINE_BINARY_FORMAT_2 (affine)

TMPL struct rounding_helper<Affine >: public rounding_helper<C> {};

STYPE_TO_TYPE(TMPL,scalar_type,Affine,C);
UNARY_RETURN_TYPE(TMPL,Re_op,Affine,Real_affine);
UNARY_RETURN_TYPE(TMPL,abs_op,Affine,Abs_affine);

/******************************************************************************
* Basic routines
******************************************************************************/

TMPL inline C base (const Affine& z) {
  return z.cst; }
TMPL inline VC slope (const Affine& z) {
  return z.lin; }
TMPL inline VR limits (const Affine& z) {
  return z.lim; }
TMPL inline format<C> CF1 (const Affine& z) {
  return get_format (base (z)); }
TMPL inline format<VC> CF2 (const Affine& z) {
  return get_format (slope (z)); }
TMPL inline C truncate (const Affine& z) {
  typedef affine_limits_helper<C,VC> Helper;
  return base (z) + Helper::truncate (slope (z), limits (z)); }
TMPL inline VR common_limits (const Affine& z1, const Affine& z2) {
  typedef affine_limits_helper<C,VC> Helper;
  return Helper::common (limits (z1), limits (z2)); }

TMPL inline Affine make_affine (const C& cst, const VC& lin, const VC& dom) {
  typedef affine_limits_helper<C,VC> Helper;
  return Affine (cst, Helper::get_linear (lin, dom),
                 Helper::get_limits (dom)); }
TMPL inline VC domain (const Affine& z) {
  typedef affine_limits_helper<C,VC> Helper;
  return Helper::get_domain (limits (z)); }


TMPL inline Affine copy (const Affine& z) {
  return Affine (copy (base (z)), copy (slope (z)), copy (limits (z))); }
TMPL inline Affine duplicate (const Affine& z) {
  return Affine (duplicate (base (z)), duplicate (slope (z)),
                 duplicate (limits (z))); }

template<typename Op, typename C, typename VC> nat
unary_hash (const Affine& x) {
  nat h= Op::op (base (x));
  h += (h<<3) ^ (h<<11) ^ (h>>29) ^ Op::op (slope (x));
  return (h<<3) ^ (h<<11) ^ (h>>29) ^ Op::op (limits (x));
}

template<typename Op, typename C, typename VC> bool
binary_test (const Affine& x1, const Affine& x2) {
  return Op::op (base (x1), base (x2)) &&
         Op::op (slope (x1), slope (x2)) &&
         Op::op (limits (x1), limits (x2));
}

template<typename CS, typename CD, typename DS, typename DD> affine<CD,DD>
map (const function_1<CD,Argument(CS) >& funC,
     const function_1<DD,Argument(DS) >& funD,
     const affine<CS,DS>& z,
     const format<CD>& fmC,
     const format<DD>& fmD) {
  return affine<CD,DD> (funC (base (z)), funD (slope (z)));
}

TRUE_IDENTITY_OP_SUGAR(TMPL,Affine)
EXACT_IDENTITY_OP_SUGAR(TMPL,Affine)
HARD_IDENTITY_OP_SUGAR(TMPL,Affine)

TMPL syntactic
flatten (const Affine& z) {
  return apply ("affine", flatten (base (z)), flatten (slope (z)),
                flatten (limits (z)));
}

TMPL
struct binary_helper<Affine >: public void_binary_helper<Affine > {
  static inline string short_type_name () {
    return "Af" * Short_type_name (C) * Short_type_name (VC); }
  static inline generic full_type_name () {
    return gen ("Affine", Full_type_name (C)); }
  static inline generic disassemble (const Affine& z) {
    return gen_vec (as<generic> (base (z)),
                    as<generic> (slope (z)),
                    as<generic> (limits (z))); }
  static inline Affine assemble (const generic& v) {
    return Affine (as<C> (vector_access (v, 0)),
                   as<VC> (vector_access (v, 1)),
                   as<VR> (vector_access (v, 2))); }
  static inline void write (const port& out, const Affine& z) {
    binary_write<C> (out, base (z));
    binary_write<VC> (out, slope (z));
    binary_write<VC> (out, limits (z)); }
  static inline Affine read (const port& in) {
    C  b= binary_read<C> (in);
    VC s= binary_read<VC> (in);
    VR d= binary_read<VR> (in);
    return Affine (b, s, d); }
};

/******************************************************************************
* Constants
******************************************************************************/

TMPL inline void set_nan (Affine& x) { x= Affine (Nan (C), 0); }
TMPL inline void set_maximal (Affine& x) { x= Affine (Maximal (C), 0); }
TMPL inline void set_minimal (Affine& x) { x= Affine (Minimal (C), 0); }
TMPL inline void set_infinity (Affine& x) { x= Affine (Infinity (C), 0); }
TMPL inline void set_fuzz (Affine& x) { x= Affine (Fuzz (C), 0); }
TMPL inline void set_smallest (Affine& x) { x= Affine (Smallest (C), 0); }
TMPL inline void set_largest (Affine& x) { x= Affine (Largest (C), 0); }
TMPL inline void set_accuracy (Affine& x) { x= Affine (Accuracy (C), 0); }
TMPL inline void set_log2 (Affine& x) { x= Affine (Log2 (C), 0); }
TMPL inline void set_pi (Affine& x) { x= Affine (Pi (C), 0); }
TMPL inline void set_euler (Affine& x) { x= Affine (Euler (C), 0); }
TMPL inline void set_imaginary (Affine& x) { x= Affine (Imaginary (C), 0); }

TMPL inline Affine times_infinity (const Affine& x) {
  return x * infinity_cst<Affine > (); }

/******************************************************************************
* Basic arithmetic
******************************************************************************/

TMPL inline Affine
operator - (const Affine& z) {
  return Affine (-base (z), -slope (z), limits (z));
}

TMPL inline Affine
operator + (const Affine& z1, const Affine& z2) {
  VR lim= common_limits (z1, z2);
  return Affine (base (z1) + base (z2), slope (z1) + slope (z2), lim);
}

TMPL inline Affine
operator + (const C& z1, const Affine& z2) {
  return Affine (z1 + base (z2), slope (z2), limits (z2));
}

TMPL inline Affine
operator + (const Affine& z1, const C& z2) {
  return Affine (base (z1) + z2, slope (z1), limits (z1));
}

TMPL inline Affine
operator - (const Affine& z1, const Affine& z2) {
  return Affine (base (z1) - base (z2), slope (z1) - slope (z2));
}

TMPL inline Affine
operator - (const C& z1, const Affine& z2) {
  return Affine (z1 - base (z2), slope (z2), limits (z2));
}

TMPL inline Affine
operator - (const Affine& z1, const C& z2) {
  return Affine (base (z1) - z2, slope (z1), limits (z1));
}

TMPL inline Affine
operator * (const Affine& z1, const Affine& z2) {
  VR lim= common_limits (z1, z2);
  return Affine (base (z1) * base (z2),
                 truncate (z1) * slope (z2) + slope (z1) * truncate (z2), lim);
}

TMPL inline Affine
operator * (const C& z1, const Affine& z2) {
  return Affine (z1 * base (z2), z1 * slope (z2), limits (z2));
}

TMPL inline Affine
operator * (const Affine& z1, const C& z2) {
  return Affine (base (z1) * z2, slope (z1) * z2, limits (z1));
}

TMPL inline Affine
square (const Affine& z) {
  return Affine (square (base (z)),
                 (base (z) + truncate (z)) * slope (z), limits (z));
}

TMPL inline Affine
invert (const Affine& z) {
  ERROR ("not yet implemented");
  //C inv= invert (base (z));
  //return Affine (inv, - square (inv) * slope (z));
}

TMPL inline Affine
operator / (const Affine& z1, const Affine& z2) {
  ERROR ("not yet implemented");
  //C inv= invert (base (z2));
  //return Affine (inv * base (z1),
  //               inv * (slope (z1) - inv * base (z1) * slope (z2)));
}

TMPL inline Affine
operator / (const C& z1, const Affine& z2) {
  ERROR ("not yet implemented");
  //C inv= invert (base (z2));
  //return Affine (inv * z1, - square (inv) * z1 * slope (z2));
}

TMPL inline Affine
operator / (const Affine& z1, const C& z2) {
  C inv= invert (z2);
  return Affine (inv * base (z1), inv * slope (z1), limits (z1));
}

ARITH_INT_SUGAR(TMPL,Affine);

/******************************************************************************
* Elementary functions
******************************************************************************/

/*
TMPL inline Affine
sqrt (const Affine& z) {
  C s= sqrt (base (z));
  return Affine (s, slope (z) / (2 * s));
}

TMPL inline Affine
exp (const Affine& z) {
  C e= exp (base (z));
  return Affine (e, e * slope (z));
}

TMPL inline Affine
log (const Affine& z) {
  return Affine (log (base (z)), slope (z) / base (z));
}

template<typename C, typename VC, typename K> inline Affine
pow (const Affine& z, const K& y) {
  return exp (Affine (y) * log (z));
}

TMPL inline Affine
cos (const Affine& z) {
  return Affine (cos (base (z)), -sin (base (z)) * slope (z));
}

TMPL inline Affine
sin (const Affine& z) {
  return Affine (sin (base (z)), cos (base (z)) * slope (z));
}

TMPL inline Affine
tan (const Affine& z) {
  return sin (z) / cos (z);
}

TMPL inline Affine
cosh (const Affine& z) {
  return Affine (cosh (base (z)), sinh (base (z)) * slope (z));
}

TMPL inline Affine
sinh (const Affine& z) {
  return Affine (sinh (base (z)), cosh (base (z)) * slope (z));
}

TMPL Affine
tanh (const Affine& z) {
  return sinh (z) / cosh (z);
}

TMPL Affine
acos (const Affine& z) {
  C df= -1 / sqrt (1 - square (base (z)));
  return Affine (acos (base (z)), df * slope (z));
}

TMPL Affine
asin (const Affine& z) {
  C df= 1 / sqrt (1 - square (base (z)));
  return Affine (asin (base (z)), df * slope (z));
}

TMPL Affine
atan (const Affine& z) {
  C df= 1 / (1 + square (base (z)));
  return Affine (asin (base (z)), df * slope (z));
}

TMPL Affine
atan2 (const Affine& y, const Affine& x) {
  ERROR ("not yet implemented");
  return Affine (atan2 (base (y), base (x)), 0);
}

HYPOT_SUGAR(TMPL,Affine)
INV_TRIGO_SUGAR(TMPL,Affine)
INV_HYPER_SUGAR(TMPL,Affine)
ARG_HYPER_SUGAR(TMPL,Affine)
*/

/******************************************************************************
* Comparisons and other order related functions
******************************************************************************/

TMPL inline bool
operator <= (const Affine& z1, const Affine& z2) {
  return base (z1) < base (z2) || z1 == z2;
}

TMPL inline bool
operator < (const Affine& z1, const Affine& z2) {
  return base (z1) < base (z2);
}

TMPL inline bool
operator >= (const Affine& z1, const Affine& z2) {
  return base (z1) > base (z2) || z1 == z2;
}

TMPL inline bool
operator > (const Affine& z1, const Affine& z2) {
  return base (z1) > base (z2);
}

TMPL inline Affine
min (const Affine& z1, const Affine& z2) {
  if (z1 <= z2) return z1; else return z2;
}

TMPL Affine
max (const Affine& z1, const Affine& z2) {
  if (z1 >= z2) return z1; else return z2;
}

EQUAL_INT_SUGAR(TMPL,Affine);
COMPARE_INT_SUGAR(TMPL,Affine);

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Affine& z) {
  return is_finite (base (z)) && is_finite (slope (z)); }
TMPL inline bool is_nan (const Affine& z) {
  return is_nan (base (z)) || is_nan (slope (z)); }
TMPL inline bool is_infinite (const Affine& z) {
  return !is_nan (z) && (is_infinite (base (z)) || is_infinite (slope (z))); }
TMPL inline bool is_fuzz (const Affine& z) {
  return !is_nan (z) && (is_fuzz (base (z)) || is_fuzz (slope (z))); }
TMPL inline bool is_reliable (const Affine& z) {
  return is_reliable (base (z)); }

TMPL inline Affine change_precision (const Affine& z, xnat prec) {
  return Affine (change_precision (base (z), prec),
                 change_precision (slope (z), prec),
                 change_precision (limits (z), prec)); }
TMPL inline xnat precision (const Affine& z) {
  return precision (base (z)); }
TMPL inline Affine additive_error (const Affine& z) {
  return Affine (additive_error (base (z)),
                 additive_error (slope (z)),
                 limits (z)); }

TMPL inline double magnitude (const Affine& z) {
  return magnitude (base (z)); }
TMPL inline xint exponent (const Affine& z) {
  return exponent (base (z)); }
TMPL inline Affine operator << (const Affine& z, const xint& s) {
  return Affine (incexp2 (base (z), s), incexp2 (slope (z), s), limits (z)); }
TMPL inline Affine operator >> (const Affine& z, const xint& s) {
  return Affine (decexp2 (base (z), s), decexp2 (slope (z), s), limits (z)); }

TMPL inline Affine sharpen (const Affine& z) {
  return Affine (sharpen (base (z)), sharpen (slope (z)), limits (z)); }
template<typename C, typename VC, typename C2, typename VC2> inline Affine
blur (const Affine& z, const affine<C2,VC2>& r) {
  return Affine (blur (base (z), base (r)), blur (slope (z), slope (r)),
                 limits (z)); }

/******************************************************************************
* Complex affine numbers
******************************************************************************/

TMPL inline Abs_affine
abs (const Affine& z) {
  // FIXME: this is not really nice
  return Abs_affine (abs (base (z)), abs (slope (z)), limits (z));
}

TMPL inline Real_affine
Re (const Affine& z) {
  return Real_affine (Re (base (z)), Re (slope (z)), limits (z));
}

TMPL inline Real_affine
Im (const Affine& z) {
  return Real_affine (Im (base (z)), Im (slope (z)), limits (z));
}

#undef TMPL_DEF
#undef TMPL
#undef R
#undef VR
#undef Affine
#undef Affine_rep
} // namespace mmx
#endif // __MMX_AFFINE_HPP