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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : ball.hpp
* DESCRIPTION: Balls
* COPYRIGHT  : (C) 2006  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_BALL_HPP
#define __MMX_BALL_HPP
#include <numerix/ball_rounded.hpp>
#include <numerix/ball_rough.hpp>
#include <numerix/ball_infinities.hpp>
namespace mmx {
#define Ball_variant(C) typename ball_variant_helper<C>::BV
#define TMPL_DEF template<typename C, \
  typename R=Default_radius_type(C),  \
  typename V=Ball_variant(C) >
#define TMPL template<typename C,typename R,typename V>
#define Ball ball<C,R,V>
#define Abs_ball ball<Abs_type(C),R,Abs_type(V) >

/******************************************************************************
* The ball type
******************************************************************************/

TMPL_DEF
class ball {
MMX_ALLOCATORS
private:
  C c;
  R r;

public:
  typedef implementation<ball_rounding,V> Rnd;
  typedef Round_up(R) Up;

  inline ball () {}
  template<typename T> ball (const T& c2):
    c (as<C> (c2)), r (0) {
      Rnd::add_additive_error (*this); }
  template<typename CT, typename RT, typename VT>
  inline ball (const ball<CT,RT,VT>& z):
    c (as<C> (center (z))), r (as<R> (radius (z))) {
      XVERIFY (r >= 0, "negative radius", r);
      // FIXME: we should pre-correct the radius of z as well
      Rnd::add_additive_error (*this); }
  inline ball (const C& x):
    c (x), r (0) {}
  inline ball (const C& c2, const R& r2):
    c (c2), r (r2) {
      XVERIFY (r >= 0, "negative radius", r); }
  inline ball (const C& c2, const R& r2, bool adjust):
    c (c2), r (r2) {
      XVERIFY (r >= 0, "negative radius", r);
      if (adjust) Rnd::add_additive_error (*this); }
  inline ball (const Real_type(Ball)& x, const Real_type(Ball)& y):
    c (gaussian (center (x), center (y))),
    r (Up::hypot (radius (x), radius (y))) {
      XVERIFY (r >= 0, "negative radius", r);
      Rnd::add_additive_error (*this); }
  friend C  center LESSGTR (const Ball& z);
  friend R  radius LESSGTR (const Ball& z);
  friend C& center LESSGTR (Ball& z);
  friend R& radius LESSGTR (Ball& z);
  inline Ball& operator <<= (const xint& shift);
  inline Ball& operator >>= (const xint& shift);
};

template<> struct symbolic_type_information<ball<floating<> > > {
  static const nat id= SYMBOLIC_BALL; };

UNARY_RETURN_TYPE(TMPL,center_op,Ball,C);
UNARY_RETURN_TYPE(TMPL,radius_op,Ball,R);
UNARY_RETURN_TYPE(TMPL,abs_op,Ball,Abs_ball);
TMPL struct rounding_helper<Ball >: public rounding_helper<R> {};
TMPL struct projective_helper<Ball >: public projective_helper<C> {};

/******************************************************************************
* Standard constructors and accessors
******************************************************************************/

template<typename C, typename R, typename V, typename C2, typename R2>
struct make_ball_helper<Ball,C2,R2> {
  static inline Ball val (const C2& c, const R2& r) {
    return Ball (c, r); }
};

template<typename C, typename R, typename V, typename C2>
struct make_interval_helper<Ball,C2> {
  static inline Ball val (const C2& l, const C2& r) {
    typedef Round_up(C2) Up;
    VERIFY (l <= r, "bad interval");
    C2 cc= decexp2 (l+r, 1);
    R  rr= next_above (as<R> (max (Up::sub (r, cc), Up::sub (cc, l))));
    return Ball (cc, rr);
  }
};

TMPL inline C center (const Ball& z) { return z.c; }
TMPL inline R radius (const Ball& z) { return z.r; }
TMPL inline C& center (Ball& z) { return z.c; }
TMPL inline R& radius (Ball& z) { return z.r; }

TMPL inline C lower (const Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  return Rnd::lower (z); }
TMPL inline C upper (const Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  return Rnd::upper (z); }
TMPL inline R abs_down (const Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  return Rnd::abs_down (z); }
TMPL inline R abs_up (const Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  return Rnd::abs_up (z); }
TMPL inline R bnd_down (const Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  return Rnd::bnd_down (z); }
TMPL inline R bnd_up (const Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  return Rnd::bnd_up (z); }

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

TMPL inline void add_rounding_error (Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  Rnd::add_rounding_error (z); }
TMPL inline void add_additive_error (Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  Rnd::add_additive_error (z); }
TMPL inline void add_multiplicative_error (Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  Rnd::add_multiplicative_error (z); }
TMPL inline void add_elementary_error (Ball& z) {
  typedef implementation<ball_rounding,V> Rnd;
  Rnd::add_elementary_error (z); }

TMPL inline Ball copy (const Ball& z) {
  return Ball (copy (center (z)), copy (radius (z))); }
TMPL inline Ball duplicate (const Ball& z) {
  return Ball (duplicate (center (z)), duplicate (radius (z))); }
TMPL inline bool is_exact_zero (const Ball& z) {
  return is_exact_zero (center (z)) && is_exact_zero (radius (z)); }

template<typename Op, typename C, typename R, typename V> nat
unary_hash (const Ball& x) {
  nat h= Op::op (center (x));
  return (h<<3) ^ (h<<11) ^ (h>>29) ^ Op::op (radius (x));
}

template<typename Op, typename C, typename R, typename V> bool
binary_test (const Ball& x1, const Ball& x2) {
  return Op::op (center (x1), center (x2)) &&
         Op::op (radius (x1), radius (x2));
}

EXACT_IDENTITY_OP_SUGAR(TMPL,Ball)
HARD_IDENTITY_OP_SUGAR(TMPL,Ball)

TMPL syntactic
flatten (const Ball& z) {
  return syn ("ball", flatten (center (z)), flatten (radius (z)));
}

TMPL syntactic
flatten (const Ball& z, xnat dd) {
  return syn ("ball", flatten (center (z), dd), flatten (radius (z), dd));
}

TMPL syntactic
flatten_range (const Ball& z) {
  mmx_local_bit_precision tmp (precision (center (z)));
  if (is_nan (z)) return Nan (syntactic);
  if (is_fuzz (z)) return Fuzz (syntactic);
  if (is_exact_zero (z)) return 0;
  return flatten_range (lower (z), upper (z));
}

TMPL syntactic
flatten_range (const Ball& z, xnat dd) {
  mmx_local_bit_precision tmp (precision (center (z)));
  if (is_nan (z)) return Nan (syntactic);
  if (is_fuzz (z)) return Fuzz (syntactic);
  if (is_exact_zero (z)) return 0;
  return flatten_range (lower (z), upper (z), dd);
}

template<typename R, typename V> inline syntactic
flatten (const ball<double,R,V>& z) {
  return flatten_range (z); }
template<typename R, typename V> inline syntactic
flatten (const ball<double,R,V>& z, xnat dd) {
  return flatten_range (z, dd); }
template<typename FV, typename R, typename V> inline syntactic
flatten (const ball<floating<FV>,R,V>& z) {
  return flatten_range (z); }
template<typename FV, typename R, typename V> inline syntactic
flatten (const ball<floating<FV>,R,V>& z, xnat dd) {
  return flatten_range (z, dd); }

TMPL
struct binary_helper<Ball >: public void_binary_helper<Ball > {
  static inline string short_type_name () {
    return "Ba" * Short_type_name (C) * Short_type_name (R); }
  static inline generic full_type_name () {
    return gen ("Ball", Full_type_name (C), Full_type_name (R)); }
  static inline generic disassemble (const Ball& z) {
    return gen_vec (as<generic> (center (z)), as<generic> (radius (z))); }
  static inline Ball assemble (const generic& v) {
    return Ball (as<C> (vector_access (v, 0)),
                 as<R> (vector_access (v, 1))); }
  static inline void write (const port& out, const Ball& z) {
    binary_write<C> (out, center (z));
    binary_write<R> (out, radius (z)); }
  static inline Ball read (const port& in) {
    C c= binary_read<C> (in);
    R r= binary_read<R> (in);
    return Ball (c, r, false); }
};

template<typename TC, typename TR, typename TV,
         typename FC, typename FR, typename FV> inline void
set_as (ball<TC,TR,TV>& r, const ball<FC,FR,FV>& z) {
  typedef implementation<ball_rounding,TV> Rnd;
  set_as (center (r), center (z));
  set_as (radius (r), radius (z));
  Rnd::add_additive_error (r);
}

template<typename TC, typename TR, typename TV, typename F> inline void
set_as (ball<TC,TR,TV>& r, const F& z) {
  typedef implementation<ball_rounding,TV> Rnd;
  set_as (center (r), z);
  set_as (radius (r), 0);
  Rnd::add_additive_error (r);
}

/******************************************************************************
* Comparisons
******************************************************************************/

TMPL inline bool
is_zero (const Ball& z) {
  return !is_nan (center (z)) && !is_nan (radius (z)) &&
         abs_down (z) == 0;
}

TMPL inline bool
is_non_zero (const Ball& z) {
  return is_nan (center (z)) || is_nan (radius (z)) ||
         abs_down (z) != 0;
}

TMPL inline bool
is_negative_or_zero (const Ball& z) {
  return !is_nan (center (z)) && bnd_down (z) <= 0;
}

TMPL inline bool
is_negative (const Ball& z) {
  return !is_nan (center (z)) && bnd_up (z) < 0;
}

TMPL inline bool
is_positive_or_zero (const Ball& z) {
  return !is_nan (center (z)) && bnd_up (z) >= 0;
}

TMPL inline bool
is_positive (const Ball& z) {
  return !is_nan (center (z)) && bnd_down (z) > 0;
}

TMPL bool
operator == (const Ball& z1, const Ball& z2) {
  if (is_nan (z1) || is_nan (z2)) return is_nan (z1) && is_nan (z2);
  Ball d= z1 - z2;
  return abs_down (d) == 0;
}

TMPL bool
operator != (const Ball& z1, const Ball& z2) {
  if (is_nan (z1) || is_nan (z2)) return !(is_nan (z1) && is_nan (z2));
  Ball d= z1 - z2;
  return abs_down (d) != 0;
}

TMPL inline bool
operator <= (const Ball& z1, const Ball& z2) {
  return is_negative_or_zero (z1 - z2);
}

TMPL inline bool
operator < (const Ball& z1, const Ball& z2) {
  return is_negative (z1 - z2);
}

TMPL inline bool
operator >= (const Ball& z1, const Ball& z2) {
  return is_positive_or_zero (z1 - z2);
}

TMPL inline bool
operator > (const Ball& z1, const Ball& z2) {
  return is_positive (z1 - z2);
}

TMPL inline bool
included (const Ball& b1, const Ball& b2) {
  // returns true if b1 is provably included in b2
  return radius (b1 - Ball (center (b2))) <= radius (b2);
}

TMPL nat
hash (const Ball& z) {
  (void) z;
  // return hash (center (z));
  return 123123123;
}

EQUAL_INT_SUGAR(TMPL,Ball);
COMPARE_INT_SUGAR(TMPL,Ball);

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

TMPL inline void set_nan (Ball& x) { x= Ball (Nan (C), 0); }
TMPL inline void set_maximal (Ball& x) { x= Ball (Maximal (C), 0); }
TMPL inline void set_minimal (Ball& x) { x= Ball (Minimal (C), 0); }
TMPL inline void set_infinity (Ball& x) { x= Ball (Infinity (C), 0); }
TMPL inline void set_fuzz (Ball& x) { x= Ball (0, Infinity (R)); }
TMPL inline void set_smallest (Ball& x) { x= Ball (Smallest (C), 0); }
TMPL inline void set_largest (Ball& x) { x= Ball (Largest (C), 0); }
TMPL inline void set_accuracy (Ball& x) { x= Ball (Accuracy (C), 0); }
TMPL inline void set_log2 (Ball& x) { x= Ball (Log2 (C), 0, true); }
TMPL inline void set_pi (Ball& x) { x= Ball (Pi (C), 0, true); }
TMPL inline void set_euler (Ball& x) { x= Ball (Euler (C), 0, true); }
TMPL inline void set_imaginary (Ball& x) {
  x= Ball (Real_type (Ball) (0), Real_type (Ball) (1)); }
TMPL inline Ball times_infinity (const Ball& x) {
  return x * infinity_cst<Ball > (); }

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

TMPL inline Ball
operator - (const Ball& z) {
  typedef implementation<ball_additive,V> Impl;
  Ball d; Impl::neg (d, z); return d;
}

TMPL inline Ball
operator + (const Ball& z1, const Ball& z2) {
  typedef implementation<ball_additive,V> Impl;
  Ball d; Impl::add (d, z1, z2); return d;
}

TMPL inline Ball
operator + (const Ball& z1, const C& z2) {
  typedef implementation<ball_additive,V> Impl;
  Ball d; Impl::add (d, z1, z2); return d;
}

TMPL inline Ball
operator + (const C& z1, const Ball& z2) {
  typedef implementation<ball_additive,V> Impl;
  Ball d; Impl::add (d, z1, z2); return d;
}

TMPL inline Ball
operator - (const Ball& z1, const Ball& z2) {
  typedef implementation<ball_additive,V> Impl;
  Ball d; Impl::sub (d, z1, z2); return d;
}

TMPL inline Ball
operator - (const Ball& z1, const C& z2) {
  typedef implementation<ball_additive,V> Impl;
  Ball d; Impl::sub (d, z1, z2); return d;
}

TMPL inline Ball
operator - (const C& z1, const Ball& z2) {
  typedef implementation<ball_additive,V> Impl;
  Ball d; Impl::sub (d, z1, z2); return d;
}

TMPL inline Ball
operator * (const Ball& z1, const Ball& z2) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::mul (d, z1, z2); return d;
}

TMPL inline Ball
operator * (const Ball& z1, const C& z2) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::mul (d, z1, z2); return d;
}

TMPL inline Ball
operator * (const C& z1, const Ball& z2) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::mul (d, z1, z2); return d;
}

TMPL inline Ball
square (const Ball& z) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::square (d, z); return d;
}

TMPL inline Ball
invert (const Ball& z) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::invert (d, z); return d;
}

TMPL inline Ball
operator / (const Ball& z1, const Ball& z2) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::div (d, z1, z2); return d;
}

TMPL inline Ball
operator / (const Ball& z1, const C& z2) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::div (d, z1, z2); return d;
}

TMPL inline Ball
operator / (const C& z1, const Ball& z2) {
  typedef implementation<ball_multiplicative,V> Impl;
  Ball d; Impl::div (d, z1, z2); return d;
}

ARITH_INT_SUGAR(TMPL,Ball);

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

template<typename V> struct ball_complex;

template<typename C, typename Real, typename V>
struct elementary_variant {
  typedef ball_complex<V> EV;
};

template<typename C, typename V>
struct elementary_variant<C,C,V> {
  typedef V EV;
};

#define Elementary_variant(C,V) \
  typename elementary_variant<C,Real_type(C),V>::EV

TMPL inline Ball
sqrt (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_root,EV> Impl;
  Ball d; Impl::sqrt (d, z); return d;
}

TMPL inline Ball
hypot (const Ball& x, const Ball& y) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_root,EV> Impl;
  Ball d; Impl::hypot (d, x, y); return d;
}

TMPL inline Ball
exp (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::exp (d, z); return d;
}

TMPL inline Ball
log (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::log (d, z); return d;
}

TMPL inline Ball
pow (const Ball& z, const Ball& y) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::pow (d, z, y); return d;
}

TMPL inline Ball
pow (const Ball& z, const int& y) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::pow (d, z, Ball (y)); return d;
}

TMPL inline Ball
pow (const int& z, const Ball& y) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::pow (d, Ball (z), y); return d;
}

TMPL inline Ball
cos (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::cos (d, z); return d;
}

TMPL inline Ball
sin (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::sin (d, z); return d;
}

TMPL inline Ball
tan (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::tan (d, z); return d;
}

TMPL inline Ball
cosh (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::cosh (d, z); return d;
}

TMPL inline Ball
sinh (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::sinh (d, z); return d;
}

TMPL inline Ball
tanh (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::tanh (d, z); return d;
}

TMPL inline Ball
acos (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::acos (d, z); return d;
}

TMPL inline Ball
asin (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::asin (d, z); return d;
}

TMPL inline Ball
atan (const Ball& z) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::atan (d, z); return d;
}

TMPL inline Ball
atan2 (const Ball& y, const Ball& x) {
  typedef Elementary_variant(C,V) EV;
  typedef implementation<ball_elementary,EV> Impl;
  Ball d; Impl::atan2 (d, y, x); return d;
}

INV_TRIGO_SUGAR(TMPL,Ball)
INV_HYPER_SUGAR(TMPL,Ball)
ARG_HYPER_SUGAR(TMPL,Ball)

/******************************************************************************
* Other order related functions
******************************************************************************/

TMPL inline Abs_ball
abs (const Ball& z) {
  typedef implementation<ball_abs,V> Impl;
  Abs_ball d; Impl::abs (d, z); return d;
}

TMPL inline Ball
min (const Ball& z1, const Ball& z2) {
  typedef implementation<ball_ordered,V> Impl;
  Ball d; Impl::min (d, z1, z2); return d;
}

TMPL inline Ball
max (const Ball& z1, const Ball& z2) {
  typedef implementation<ball_ordered,V> Impl;
  Ball d; Impl::max (d, z1, z2); return d;
}

TMPL inline Ball
floor (const Ball& z) {
  typedef implementation<ball_ordered,V> Impl;
  Ball d; Impl::floor (d, z); return d;
}

TMPL inline Ball
trunc (const Ball& z) {
  typedef implementation<ball_ordered,V> Impl;
  Ball d; Impl::trunc (d, z); return d;
}

TMPL inline Ball
ceil (const Ball& z) {
  typedef implementation<ball_ordered,V> Impl;
  Ball d; Impl::ceil (d, z); return d;
}

TMPL inline Ball
round (const Ball& z) {
  typedef implementation<ball_ordered,V> Impl;
  Ball d; Impl::round (d, z); return d;
}

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

TMPL inline bool is_finite (const Ball& z) {
  return is_finite (center (z)) && is_finite (radius (z)); }
TMPL inline bool is_infinite (const Ball& z) {
  return is_infinite (center (z)) && is_finite (radius (z)); }
TMPL inline bool is_fuzz (const Ball& z) {
  return !is_nan (center (z)) && is_infinite (radius (z)); }
TMPL inline bool is_nan (const Ball& z) {
  return is_nan (center (z)) || is_nan (radius (z)); }
TMPL inline bool is_reliable (const Ball& z) {
  (void) z; return true; }

TMPL inline Ball change_precision (const Ball& z, xnat prec) {
  return Ball (change_precision (center (z), prec), radius (z)); }
TMPL inline xnat precision (const Ball& x) {
  return precision (center (x)); }
TMPL inline Abs_ball additive_error (const Ball& x) {
  return Abs_ball (additive_error (center (x))); }
TMPL inline Abs_ball multiplicative_error (const Ball& x) {
  return Abs_ball (multiplicative_error (center (x))); }
TMPL inline Abs_ball elementary_error (const Ball& x) {
  return Abs_ball (elementary_error (center (x))); }

TMPL inline xint exponent (const Ball& z) {
  return max (exponent (center (z)), exponent (radius (z))); }
TMPL inline double magnitude (const Ball& z) {
  return max (magnitude (center (z)), magnitude (radius (z))); }
TMPL inline Ball operator << (const Ball& z, const xint& shift) {
  typedef implementation<ball_shift,V> Impl;
  Ball d; Impl::shiftl (d, z, shift); return d; }
TMPL inline Ball operator >> (const Ball& z, const xint& shift) {
  typedef implementation<ball_shift,V> Impl;
  Ball d; Impl::shiftr (d, z, shift); return d; }
TMPL inline Ball& Ball::operator <<= (const xint& shift) {
  typedef implementation<ball_shift,V> Impl;
  Impl::shiftl (*this, shift); return *this; }
TMPL inline Ball& Ball::operator >>= (const xint& shift) {
  typedef implementation<ball_shift,V> Impl;
  Impl::shiftr (*this, shift); return *this; }

TMPL inline Ball sharpen (const Ball& z) {
  return Ball (center (z)); }
TMPL inline Ball blur (const Ball& z, const R& r) {
  typedef Round_up(R) Up;
  return Ball (center (z), Up::add (radius (z), r)); }
template<typename C, typename R, typename V, typename K> inline Ball
blur (const Ball& z, const ball<K,R,V>& r) {
  typedef Round_up(R) Up;
  return Ball (center (z), Up::add (radius (z), abs_up (r))); }

/******************************************************************************
* Extras for glue
******************************************************************************/

#define mmx_ball(R,C) ball<C >

template<typename R, typename C> inline ball<C>
make_mmx_ball (const C& c, const R& r) {
  return ball<C> (c, as<Default_radius_type(C) > (r), true);
}

template<typename C> inline ball<C>
make_mmx_ball (const C& c) {
  return ball<C> (c, as<Default_radius_type(C) > (promote (0, abs (c))));
}

template<typename C> inline Abs_type(C)
mmx_radius (const ball<C>& z) {
  return as<Abs_type(C) > (radius (z));
}

#undef TMPL_DEF
#undef TMPL
#undef Ball
#undef Abs_ball
} // namespace mmx
#endif // __MMX_BALL_HPP