Interval_fcts.hpp

Go to the documentation of this file.

#ifndef realroot_ARITHM_INTERVAL_C
#define realroot_ARITHM_INTERVAL_C

#include <string.h>
#include <string>
#include <realroot/texp_sup.hpp>
namespace mmx {

#ifndef MMX_SIGN_FCT
#define MMX_SIGN_FCT
  template<typename T> inline
  int sign(const T& x) {return ( x < T(0) ? -1 : (T(0) < x ? 1 : 0) ); }
#endif

template<class C, int r> inline bool 
with_plus_sign(const Interval<C,r> & I) { return true;}

namespace numerics
{
  template<class T, int r>
  struct interval_base
  {
    struct rnd{};
    inline static T dwmul( const T& a, const T& b ) { return a*b; };
    inline static T upmul( const T& a, const T& b ) { return a*b; };
    inline static T dwdiv( const T& a, const T& b ) { return a/b; };
    inline static T updiv( const T& a, const T& b ) { return a/b; };
    inline static T dwadd( const T& a, const T& b ) { return a+b; };
    inline static T upadd( const T& a, const T& b ) { return a+b; };
    inline static T dwsub( const T& a, const T& b ) { return a-b; };
    inline static T upsub( const T& a, const T& b ) { return a-b; };
  };
  
  template<class T>
  struct interval_base<T,2> : interval_base<T,0> 
  {
    struct rnd{};
    inline static T dwmul( const T& a, const T& b ) { return rup::dwmul(a,b); };
    inline static T dwdiv( const T& a, const T& b ) { return rup::dwdiv(a,b); };
    inline static T dwadd( const T& a, const T& b ) { return rup::dwadd(a,b); };
    inline static T dwsub( const T& a, const T& b ) { return rup::dwsub(a,b); };
  };
  
  template<class T>
  struct interval_base<T,1> : interval_base<T,0>
  {
    struct rnd{};
    inline static T upmul( const T& a, const T& b ) { return rdw::upmul(a,b); };
    inline static T updiv( const T& a, const T& b ) { return rdw::updiv(a,b); };
    inline static T upadd( const T& a, const T& b ) { return rdw::upadd(a,b); };
    inline static T upsub( const T& a, const T& b ) { return rdw::upsub(a,b); };
  };

  template<class T>
  struct interval_base<T,3> : interval_base<T,1>
  {
    struct rnd : numerics::rounding<T>
    {
      rnd() : numerics::rounding<T>( numerics::rounding<T>::rnd_dw() ) {};
    };
  };
};

template<class T, int r> struct Interval;  
namespace let
{
  template<class T, int r>
  void assign( Interval<T,r>& i, const char * s );
};

template<class C, int r> inline
Interval<C,r>::Interval() : m(0), M(0) {};
template<class C, int r> inline
Interval<C,r>::Interval( int n ) :m(n), M(n) {};
template<class C, int r> inline
Interval<C,r>::Interval( unsigned n ) :m(n), M(n) {};
template<class C, int r> inline
Interval<C,r>::Interval( const C& x ): m(x), M(x) {};
template<class C, int r> inline
Interval<C,r>::Interval( const C& a, const C& b ) {
  if ( a > b ) m = b, M = a;
  else         m = a, M = b;
};
template<class C, int r> inline
Interval<C,r>::Interval( const char * s ) { let::assign(*this,s); };
template<class T, int r> inline
T lower ( const Interval<T,r>& x ) { return x.lower(); };
template<class T, int r> inline 
T upper ( const Interval<T,r>& x ) { return x.upper(); };
template<class T, int r> inline
T median( const Interval<T,r>& x ) { return (lower(x)+upper(x))/T(2); };
template<class T, int r> inline
T width( const Interval<T,r>& x ) { return x.width(); };
template<class T, int r> inline
bool singleton( const Interval<T,r>& x ) { return x.lower() == x.upper(); };
template<class T, int r> inline
bool contain_zero( const Interval<T,r>& x ) 
{ return x.lower() <= static_cast<T>(0) && x.upper() >= static_cast<T>(0); };
template<class T, int r> inline
bool in( const T& x, const Interval<T,r>& y ) 
{ return y.lower() <= x && x <= y.upper(); };
template<class T, int r> inline
std::pair< Interval<T,r>, Interval<T,r> > 
bissect(const Interval<T,r>& x)
{
  typedef Interval<T,r> I;
  const T m(median(x));
  return std::pair<I,I>(I(x.lower(), m), I(m, x.upper()));
}

template<class T, int r> inline void 
bissect( Interval<T,r>& g, Interval<T,r>& d, const Interval<T,r>& x )
{
  const T m(median(x));
  g.define( x.lower(), m );
  d.define( m, x.upper() );
};

template<class T, int r> inline void 
hull( Interval<T,r>& v, const Interval<T,r>& a, const Interval<T,r>& b ) {
  v.define( std::min( a.lower(), b.lower() ), std::max( a.upper(), b.upper() ) );
};

template<class T, int r> inline Interval<T,r> 
hull( const Interval<T,r>& a, const Interval<T,r>& b ) {
  return Interval<T,r>( std::min(a.lower(),b.lower()), std::max(a.upper(),b.upper()) );
};

template<class T, int r> inline Interval<T,r>
min( const Interval<T,r>& a, const Interval<T,r>& b ) {
  return Interval<T,r>( std::min(a.lower(),b.lower()), std::min(a.upper(),b.upper()) );
};

template<class T, int r> inline Interval<T,r>
max( const Interval<T,r>& a, const Interval<T,r>& b ) {
  return Interval<T,r>( std::max(a.lower(),b.lower()), std::max(a.upper(),b.upper()) );
};

template<class T,  int r> inline bool 
intersectp( const Interval<T,r>& a, const Interval<T,r>& b ) {
  return lower(a) <= upper(b) && lower(b) <= upper(a);
};

template<class T, int r > inline Interval<T,r> 
intersection( const Interval<T,r>& a, const Interval<T,r>& b ) {
  return Interval<T,r>(std::max(lower(a),lower(b)),std::min(upper(a),upper(b)));
};

template<class T,  int r> inline bool 
intersect( Interval<T,r>& result, 
       const Interval<T,r>& a, const Interval<T,r>& b ) {
  if ( intersectp( a, b ) )
    {
      result.define(std::max(lower(a),lower(b)),std::min(upper(a),upper(b)));
      return true;
    };
  return false;
};

/*
template<typename T,class X> inline
void distance( T& dmin, T& dmax, const X& x, const Interval<T,r>& i )
{
  if ( x < i.min() ) { dmin = i.min()-x; dmax = i.max()-x;  return; };
  if ( x > i.max() ) { dmin = x-i.max(); dmax = x-i.min();  return; };
  dmin = x-i.min();
  dmax = i.max()-x;
  if ( dmax < dmin ) dmax = dmin;
  dmin = 0;
};
*/

template<class T,  int r>

T size( const Interval<T,r>& i  ) { return i.size(); };

template<class T,  int r>
std::ostream& operator<<( std::ostream& o, const Interval<T,r>& i )
{
  o << "[" << (i.lower()) << ", " << (i.upper()) << "]";
  return o;
};

template<class C,  int r >
void abs( Interval<C,r>& x, const Interval<C,r>& a )
{
  if ( a.upper() > -a.lower() )
    {
      x.upper() =  a.upper();
      x.lower() =  std::max(-a.lower(),0);
    }
  else
    {
      x.upper() = -a.lower();
      x.lower() =  a.upper();
    };
};

template<class C, int r> inline void 
neg( Interval<C,r>& a )
{
  C sv(a.M);
  a.M = -a.m;
  a.m = -sv;
};

template<class C, int r> inline void
neg( Interval<C,r>& a, const Interval<C,r>& b )
{
  a.M = -b.m;
  a.m = -b.M;
};

template<class C, int r> inline void 
add( Interval<C,r>& a, const C& x )
{
  a.m = a.dwadd(a.m,x);
  a.M = a.upadd(a.M,x);
};


template<class C, int r> inline void 
add( Interval<C,r>& a, const Interval<C,r>& b, const C& x ) { 
  a.m = a.dwadd(b.m,x); a.M = a.upadd(b.M,x); };

template<class C, int r> inline void 
add( Interval<C,r>& a , const C& x , const Interval<C,r>& b ) { 
  add(a,b,x); };

template<class C, int r> inline void 
sub( Interval<C,r>& a, const C& x ) { 
  a.m = a.dwsub(a.m,x); a.M = a.upsub(a.M,x); };

template<class C, int r> inline void 
sub( Interval<C,r>& a, const Interval<C,r>& b, const C& x ) { 
  a.m = a.dwsub(b.m,x); a.M = a.upsub(b.M,x); };

template<class C, int r> inline void 
sub( Interval<C,r>& a,  const C& x , const Interval<C,r>& b ) { 
  a.m = a.dwsub(x,b.M); a.M = a.upsub(x,b.m); };

template<class C,  int r> inline void 
add( Interval<C,r>& a, const Interval<C,r>& x ) { 
  a.m = a.dwadd(a.m,x.m); a.M = a.upadd(a.M,x.M); };

template<class C, int r> inline void 
add( Interval<C,r>& a, const Interval<C,r>& b, const Interval<C,r>& x ) {  
  a.m = a.dwadd(b.m,x.m); a.M = a.upadd(b.M,x.M); };

template<class C, int r> inline void 
sub( Interval<C,r>& a, const Interval<C,r>& x ) { 
a.m = a.dwsub(a.m,x.M); a.M = a.upsub(a.M,x.m); };

template<class C, int r> inline void 
sub( Interval<C,r>& a, const Interval<C,r>& b, const Interval<C,r>& x ) { 
  a.m = a.dwsub(b.m,x.M); a.M = a.upsub(b.M,x.m); };

template<class C, int r> inline void 
mul( Interval<C,r>& a, const C& x ) {
  if ( x > 0 )
    { 
      a.m = a.dwmul(a.m,x); 
      a.M = a.upmul(a.M,x); }
  else 
    { 
      C sv(a.m); 
      a.m = a.dwmul(a.M,x); 
      a.M = a.dwmul(sv,x);
    };
};


template<class C, int r> inline void 
mul( Interval<C,r>& a, const Interval<C,r>& b, const C& x ) {
  if ( &a == &b ) { mul(a,x); return; };
  if ( x > 0 )
    {
      a.m = a.dwmul(b.m,x);
      a.M = a.upmul(b.M,x);
    }
  else 
    {
      a.m = a.dwmul(b.M,x);
      a.M = a.dwmul(b.m,x);
    };
};

template<class C, int r> inline void 
mul( Interval<C,r>& a, const C& x, const Interval<C,r>& b  ) { 
  mul(a,b,x); };


template<class C, int r> inline void 
mul( Interval<C,r>& a, const Interval<C,r>& b ) {
  if ( a.m > 0 ) 
    {
      if ( b.m > 0 )
    {
      a.m = a.dwmul(a.m,b.m);
      a.M = a.upmul(a.M,b.M);
      return;
    };
      if ( b.M < 0 ) 
    {
      C sv(a.m);
      a.m = a.dwmul(a.M,b.m);
      a.M = a.upmul(sv,b.M);
      return;
    };
      a.m = a.dwmul(a.M,b.m);
      a.M = a.upmul(a.M,b.M);
      return;
    };
  
  if ( a.M < 0 )
    {
      if ( b.m > 0 )
    {
      a.m = a.dwmul(b.M,a.m);
      a.M = a.upmul(b.m,a.M);
      return;
    };
      if ( b.M < 0 )
    {
      C sv(a.m);
      a.m = a.dwmul(a.M,b.M);
      a.M = a.upmul(sv,b.m);
      return;
    };
      C sv(a.m);
      a.m = a.dwmul(a.m,b.M);
      a.M = a.upmul(sv,b.m);
      return;
    };
  
  if ( b.m > 0 ) 
    {
      a.m = a.dwmul(a.m,b.M);
      a.M = a.upmul(a.M,b.M);
      return;
    }
  if ( b.M < 0 )
    {
      C sv(a.m);
      a.m = a.dwmul(a.M,b.m);
      a.M = a.upmul(sv,b.m);
      return;
    };
  //  std::cout << "last case X\n";
  C m0(a.dwmul(a.M,b.m));
  C m1(a.dwmul(a.m,b.M));
  if ( m0 > m1 ) m0 = m1;
  C M0(a.upmul(a.M,b.M));
  C M1(a.upmul(a.m,b.m));
  if ( M0 < M1 ) M0 = M1;
  a.m = m0;
  a.M = M0;

};

template<class C, int r> inline void 
mul( Interval<C,r>& x, const Interval<C,r>& a, const Interval<C,r>& b ) {
  if ( &x == &a ) { mul(x,b); return; };
  if ( &x == &b ) { mul(x,a); return; };
  
  if ( a.m > 0 ) 
    {
      if ( b.m > 0 )
    {
      x.m = a.dwmul(a.m,b.m);
      x.M = a.upmul(a.M,b.M);
      return;
    };
      if ( b.M < 0 ) 
    {
      x.m = a.dwmul(a.M,b.m);
      x.M = a.upmul(a.m,b.M);
      return;
    };
      x.m = a.dwmul(a.M,b.m);
      x.M = a.upmul(a.M,b.M);
      return;
    };
  //std::cout << "popo\n";
  if ( a.M < 0 )
    {
      if ( b.m > 0 )
    {
      x.m = a.dwmul(b.M,a.m);
      x.M = a.upmul(b.m,a.M);
      return;
    };
      if ( b.M < 0 )
    {
      x.m = a.dwmul(a.M,b.M);
      x.M = a.upmul(a.m,b.m);
      return;
    };
       x.m = a.dwmul(a.m,b.M);
      x.M = a.upmul(a.m,b.m);
      return;
    };
  
  if ( b.m > 0 ) 
    {
      x.m = a.dwmul(a.m,b.M);
      x.M = a.upmul(a.M,b.M);
      return;
    }
  if ( b.M < 0 )
    {
      x.m = a.dwmul(a.M,b.m);
      x.M = a.upmul(a.m,b.m);
      return;
    };
  
  C m0(a.dwmul(a.M,b.m));
  C m1(a.dwmul(a.m,b.M));
  if ( m0 > m1 ) m0 = m1;
  
  C M0(a.upmul(a.M,b.M));
  C M1(a.upmul(a.m,b.m));
  if ( M0 < M1 ) M0 = M1;
  
  x.m = m0;
  x.M = M0;
};

template<class C,  int r> inline void
div( Interval<C,r>& a, const C& x ) {
  if ( x > 0 )
    {
      a.m = a.dwdiv(a.m,x);
      a.M = a.updiv(a.M,x);
    }
  else
    {
      C sv(a.m);
      a.m = a.dwdiv(a.M,x);
      a.M = a.updiv(sv,x);
    };
};

template<class C, int r> inline void 
div( Interval<C,r>& a, const Interval<C,r>& b, const C& x ) {
  if ( &a == &b ) { div(a,x); return; };
  if ( x > 0 )
    {
      a.m = a.dwdiv(b.m,x);
      a.M = a.updiv(b.M,x);
    }
  else
    {
      a.m = a.dwdiv(b.M,x);
      a.M = a.updiv(b.m,x);
    };
};

template<class C, int r> inline void 
div( Interval<C,r>& a, const C& x, const Interval<C,r>& b ) {
  if ( x > 0 )
    {
      a.m = a.dwdiv(x,b.M);
      a.M = a.updiv(x,b.m);
    }
  else
    {
      a.m = a.dwdiv(x,b.m);
      a.M = a.updiv(x,b.M);
    };
};

template<class C,  int r> inline void
div( Interval<C,r>& a, const Interval<C,r>& b ) {
  if ( a.m > 0 )
    {
      if ( b.m > 0 )
    {
      a.m = a.dwdiv(a.m,b.M);
      a.M = a.updiv(a.M,b.m);
      return;
    };
      if ( b.M < 0 )
    {
      C sv(a.m);
      a.m = a.dwdiv(a.M,b.M);
      a.M = a.updiv(sv,b.m);
      return;
    };
      if ( r == 3 )
    {
      throw typename Interval<C,r>::extended(a.updiv(a.m,b.m),a.dwdiv(a.M,b.M));
    };
    };
  
  if ( a.M < 0 )
    {
      if ( b.m > 0 )
    {
      a.m = a.dwdiv(a.m,b.m);
      a.M = a.updiv(a.M,b.M);
      return;
    };
      
      if ( b.M < 0 ) 
    {
      C sv(a.m);
      a.m = a.dwdiv(a.M,b.m);
      a.M = a.updiv(sv,b.M);
      return;
    };
      
      if ( r == 3 )
    {
      throw typename Interval<C,r>::extended(a.updiv(a.m,b.M),a.dwdiv(a.M,b.m));
    };
    }
    
  if ( b.m > 0 )
    {
      a.m = a.dwdiv(a.m,b.m);
      a.M = a.updiv(a.M,b.m);
      return ;
    };
  
  if ( b.M < 0 )
    {
      C sv(a.m);
      a.m = a.dwdiv(a.M,b.m);
      a.M = a.updiv(sv,b.m);
      return;
    };
  if ( r == 3 ) throw typename Interval<C,r>::extended(0,0);
};

template<class C, int r> inline void 
div( Interval<C,r>& x, const Interval<C,r>& a, const Interval<C,r>& b ) {
  /*
  if ( b.m > 0 )
    {
      if ( a.m > 0 ) { 
    x.m = a.dwdiv(a.m,b.M);
    x.M = a.updiv(a.M,b.m);
      }
      else           x.m = a.dwdiv(a.m,b.m);
      return;
    }
  if ( b.M < 0 )
    {
      
    };

  */
  switch( sign(a.m) )
    {
    case 1:
      if ( b.m > 0 )
    {
      x.m = a.dwdiv(a.m,b.M);
      x.M = a.updiv(a.M,b.m);
    } else 
      if ( b.M < 0 )
    {
      x.m = a.dwdiv(a.M,b.M);
      x.M = a.updiv(a.m,b.m);
    }
      else throw typename Interval<C,r>::extended(a.updiv(a.m,b.m),a.dwdiv(a.M,b.M));
      break;
    case -1:
      if ( b.m > 0 )
    {
      x.m = a.dwdiv(a.m,b.m);
      x.M = a.updiv(a.M,b.M);
    }else 
      if ( b.M < 0 ) 
    {
      //      C sv(a.M);
      x.m = a.dwdiv(a.M,b.m);
      x.M = a.updiv(a.m,b.M);
    } 
      else throw typename Interval<C,r>::extended(a.updiv(a.m,b.M),a.dwdiv(a.M,b.m));
      break;
    case 0:
      if ( b.m > 0 )
    {
      x.m = 0;
      x.M = a.updiv(a.M,b.m);
    }
      else 
    if ( b.M < 0 )
      {
        x.m = a.dwdiv(a.M,b.m);
        x.M = a.updiv(a.m,b.m);
      }
    else throw typename Interval<C,r>::extended(0,0);
      break;
    };
};




template<class C>
struct Intervals
{
  typedef Interval<C,0> simple_t;
  typedef Interval<C,1> rdw_t;
  typedef Interval<C,2> rup_t;
  typedef Interval<C,3> autoround_t;
};

namespace let 
{
  template<class T, int r> void 
  assign( Interval<T,r>& i, const char * s ) {
    int n = strlen(s)+1;
    if ( *s == '[' )
    {
      char  _s[ n ];
      std::copy(s,s+n,_s);
      char * sm = _s+1;
      char * sM = _s+1;
      while ( *sM != ',' ) sM++;
      *sM = 0;
      sM++;
      char * e = sM;
      while ( *e && *e != ']' ) e++;
      *e = 0;
      let::assign(i.m,(char*)sm);
      let::assign(i.M,(char*)sM);
    }
    else
    {
      T tmp;
      let::assign(tmp,(char*)s);
      i.m = i.M = tmp;
    };
  };
}

template<class C, int R> void 
split( Interval<C,R>& l, Interval<C,R>& r ) {
  r.M = l.M;
  r.m = (l.M+l.m)/C(2);
  l.M = r.m;
};

template<class C, int R> Interval<C,R> 
operator-(const Interval<C,R>& I) {
  Interval<C,R> res; neg(res,I); return res;
}
#define TMPL  template<class Ca, int Na, class Cb, int Nb>
#define ARG0 Interval<Ca,Na>
#define ARG1 Interval<Cb,Nb>
TMPL typename texp::sup<ARG0,ARG1>::T operator+(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; add(res,Ia,Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator-(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; sub(res,Ia,Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator*(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; mul(res,Ia,Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator/(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; div(res,Ia,Ib); return res;
}

#undef  ARG1
#undef  TMPL
#define TMPL  template<class Ca, int Na, class Cb>
#define ARG1 Cb 

TMPL typename texp::sup<ARG0,ARG1>::T operator+(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; add(res,Ia,(typename texp::sup<Ca,Cb>::T)Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator-(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; sub(res,Ia,(typename texp::sup<Ca,Cb>::T)Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator*(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; mul(res,Ia,(typename texp::sup<Ca,Cb>::T)Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator/(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; div(res,Ia,(typename texp::sup<Ca,Cb>::T)Ib); return res;
}
//----------------------------------------------------------------------
TMPL typename texp::sup<ARG0,ARG1>::T operator+(const ARG1& Ia, const ARG0& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; add(res,Ib,(typename texp::sup<Ca,Cb>::T)Ia); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator-(const ARG1& Ia, const ARG0& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; add(res,-Ib,(typename texp::sup<Ca,Cb>::T)Ia); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator*(const ARG1& Ia, const ARG0& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; mul(res,Ib,(typename texp::sup<Ca,Cb>::T)Ia); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator/(const ARG1& Ia, const ARG0& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; div(res,Ib,(typename texp::sup<Ca,Cb>::T)Ia); return res;
}

//  declare_binary_operator(TMPL,ARG1,ARG0,_add_,operator+);
//  declare_binary_operator(TMPL,ARG1,ARG0,_sub_,operator-);
//  declare_binary_operator(TMPL,ARG1,ARG0,_mul_,operator*);
//  declare_binary_operator(TMPL,ARG1,ARG0,_div_,operator/);
#undef TMPL
#undef ARG1
#define TMPL template<class Ca, int Na, class K>
#define ARG1 typename K::integer
TMPL typename texp::sup<ARG0,ARG1>::T operator+(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; add(res,Ia,(Ca)Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator-(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; sub(res,Ia,(Ca)Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator*(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; mul(res,Ia,Ib); return res;
}
TMPL typename texp::sup<ARG0,ARG1>::T operator/(const ARG0& Ia, const ARG1& Ib) {
  typename texp::sup<ARG0,ARG1>::T res; div(res,Ia,Ib); return res;
}
//   declare_binary_operator(TMPL,ARG0,ARG1,_add_,operator+);
//   declare_binary_operator(TMPL,ARG0,ARG1,_sub_,operator-);
//   declare_binary_operator(TMPL,ARG0,ARG1,_mul_,operator*);
//   declare_binary_operator(TMPL,ARG0,ARG1,_div_,operator/);

  declare_binary_operator(TMPL,ARG1,ARG0,_add_,operator+);
  declare_binary_operator(TMPL,ARG1,ARG0,_sub_,operator-);
  declare_binary_operator(TMPL,ARG1,ARG0,_mul_,operator*);
  declare_binary_operator(TMPL,ARG1,ARG0,_div_,operator/);
#undef ARG1
#undef TMPL
#define TMPL template<class Ca, int Na, class K>
#define ARG1 typename K::floating
  declare_binary_operator(TMPL,ARG0,ARG1,_add_,operator+);
  declare_binary_operator(TMPL,ARG0,ARG1,_sub_,operator-);
  declare_binary_operator(TMPL,ARG0,ARG1,_mul_,operator*);
  declare_binary_operator(TMPL,ARG0,ARG1,_div_,operator/);
  declare_binary_operator(TMPL,ARG1,ARG0,_add_,operator+);
  declare_binary_operator(TMPL,ARG1,ARG0,_sub_,operator-);
  declare_binary_operator(TMPL,ARG1,ARG0,_mul_,operator*);
  declare_binary_operator(TMPL,ARG1,ARG0,_div_,operator/);
#undef ARG1
#undef TMPL
#define TMPL template<class Ca, int Na, class K>
#define ARG1 typename K::rational
  declare_binary_operator(TMPL,ARG0,ARG1,_add_,operator+);
  declare_binary_operator(TMPL,ARG0,ARG1,_sub_,operator-);
  declare_binary_operator(TMPL,ARG0,ARG1,_mul_,operator*);
  declare_binary_operator(TMPL,ARG0,ARG1,_div_,operator/);
  declare_binary_operator(TMPL,ARG1,ARG0,_add_,operator+);
  declare_binary_operator(TMPL,ARG1,ARG0,_sub_,operator-);
  declare_binary_operator(TMPL,ARG1,ARG0,_mul_,operator*);
  declare_binary_operator(TMPL,ARG1,ARG0,_div_,operator/);
#undef ARG1
#undef ARG0
#undef TMPL

//======================================================================
} //end namespace mmx
#endif