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

Go to the documentation of this file.

#ifndef numerix_kernel_hpp
#define numerix_kernel_hpp

#include <math.h>
#include <basix/basix.hpp>
#include <basix/glue.hpp>
#include <numerix/integer.hpp>
#include <numerix/integer_meta.hpp>
#include <numerix/rational.hpp>
#include <numerix/rational_meta.hpp>
#include <numerix/floating.hpp>
#include <numerix/floating_meta.hpp>
#include <numerix/ball.hpp>
#include <numerix/interval.hpp>

//#include <realroot/mpfr_rounding_mode.hpp>
//#include <realroot/Interval.hpp>


namespace mmx {
  template<class OSTREAM> inline void
  print(OSTREAM& os, const integer& x) { os << as_charp(as_string(x)); }

  template<class OSTREAM> inline void
  print(OSTREAM& os, const rational& x) { 
    os   << as_charp(as_string(numerator(x))); 
    if (denominator(x) != 1){ 
      os <<"/";
      os << as_charp(as_string(denominator(x))); 
    }
  }


  template<class OSTREAM, class V> inline void
  print(OSTREAM& os, const floating<V>& x) { os << as_charp(as_string(x)); }



  template<class X> struct scalar_set {
    typedef X T;
    scalar_set(){}
    scalar_set(const scalar_set<X>& t) {}
    scalar_set(const X& x) {}
  };
  
  template<class X> 
  bool operator== (const scalar_set<X>& r1, const scalar_set<X>& r2) { return true; }
  template<class X> 
  bool operator!= (const scalar_set<X>& r1, const scalar_set<X>& r2) { return !(r1==r2); }
  template<class C> inline int hash(const scalar_set<C>& R){ return 1;}
  template<class C> inline int exact_hash(const scalar_set<C>& R){ return 1;}
  template<class C> inline bool 
  exact_eq(const scalar_set<C>& R1, const scalar_set<C>& R2) {
    return R1==R2;
  }
  
  template<class C> inline scalar_set<C> make_set(const C& c) { 
    return scalar_set<C>(); 
  }

  template<typename U,typename T>
  struct as_helper<scalar_set<U>,scalar_set<T> > {
    static inline 
    scalar_set<U> cv (const scalar_set<T>& x) { return scalar_set<U>(); }
  };

#define set_of(T) scalar_set< T > ()

#define DECLARE_SCALAR_SET(X,T)                                               \
inline syntactic flatten(const scalar_set<T>& Z){return syntactic(X);} \

#define bigfloat floating<>

DECLARE_SCALAR_SET("DD",double);
DECLARE_SCALAR_SET("ZZ",integer);
DECLARE_SCALAR_SET("QQ",rational);
DECLARE_SCALAR_SET("RR",bigfloat);
DECLARE_SCALAR_SET("CD",complex< double >);
DECLARE_SCALAR_SET("CC",complex< bigfloat >);

#ifndef WITH_AS
#define WITH_AS
template<typename T,typename F>
struct as_helper { static inline T cv (const F& x) { return x; } };

template<typename T,typename F> inline T as (const F& x) { return as_helper<T,F>::cv (x); }
#endif


  struct Numerix {

    typedef ::mmx::integer    integer;
    typedef ::mmx::rational   rational;
    typedef ::mmx::floating<> floating;
    typedef double            ieee;

    typedef ::mmx::interval<rational> interval_rational;
  };


  namespace texp {
    template<class R> struct kernelof_;
    //    template<> struct kernelof_<Numerix::ieee>     { typedef Numerix T; };
    template<> struct kernelof_<Numerix::integer>  { typedef Numerix T; };
    template<> struct kernelof_<Numerix::rational> { typedef Numerix T; };
    template<> struct kernelof_<Numerix::floating> { typedef Numerix T; };
  }

  //  inline int floor(double x) { return (int)std::floor(x);}

 template<class T> inline
 bool singleton(const interval<T>& x ) { return lower(x) == upper(x); };

  namespace meta {
    template<class C> struct kernelof_;
    template<> struct kernelof_< integer  >   { typedef Numerix T; };
    template<> struct kernelof_< rational >   { typedef Numerix T; };
    template<> struct kernelof_< floating<> > { typedef Numerix T; };
    
    //  template<class C> struct structureof_;
    //  template<> struct structureof_< integer  >   { typedef structure::scalar T; };
    //  template<> struct structureof_< rational >   { typedef structure::scalar T; };
    //  template<> struct structureof_< floating<> > { typedef structure::scalar T; };
    
  }
 
 

  namespace let {
    template<class A,class B> void assign(A& a, const B&b);

    inline void 
    assign(rational& q, char* s)  {q= rational(integer(s)); }

    //    template<> inline void assign(double& x, const integer& z)   {x = as_double(z);}
    template<> inline void 
    assign(double& x, const rational& q)  {x = as_double(q);}
    //    template<> inline void assign(double& x, const floating<>& r){x = as_double(r);}


    template<> inline void 
    assign(floating<>& r, const int& n) {r = floating<>(n);}

//    template<class C, class T> inline void 
//    assign( Interval<C>& a, const T & b ) {
//       C mn,mx;
//       {
//         numerics::rounding<C> rnd( numerics::rounding<C>::rnd_up() );
//         assign(mx,b);
//       }
//       
//       {
//         numerics::rounding<C> rhd( numerics::rounding<C>::rnd_dw() );
//         assign(mn,b);
//       }
//       if( mx < mn ) std::swap(mn,mx);
//       new (&a) Interval<C>(mn,mx);
//     };
  }


  inline bool less_prec(const rational& q, unsigned prec)
  {
    integer n;
    if(q<0)
      n=-(numerator(q)<<prec);
    else
      n=numerator(q)<<prec;
    return (n<denominator(q));
  }

  inline rational nearest(const floating<>& x, int prec) 
  {
    integer er = integer(1)<<prec;

    floating<> r=(x>0?x:-x), f = floor(r);
    integer a=as_integer(f);
    
    integer n0=1, n1= a, tn=a;
    integer d0=0, d1= 1, td=1;
    
    for(;d0*d1<er && r!=f;)
      {
        r=1/(r-f); f = floor(r); 
        a=as_integer(f);
        n1*=a;n1+=n0;n0=tn;tn=n1;
        d1*=a;d1+=d0;d0=td;td=d1;
      }
    return (x>0? rational(n1)/d1: rational(-n1)/d1);
  }
  template<> struct as_helper<rational,floating<> >
  {
    static inline 
    rational cv (const floating<>& x) 
    {
      return nearest(x,precision(x));
    }
  };

 template<> struct as_helper<double,floating<> > 
  {
    static inline double cv (const floating<>& x) { return as_double(x); }
  };

 //  template<> struct as_helper< Interval<double>, rational>

//   { 
//     static inline Interval<double> cv(const rational & b )
//     {
//       double mn,mx;
//       {
//         numerics::rounding<double> rnd( numerics::rounding<double>::rnd_up() );
//      let::assign(mx,b);
//       }
      
//       {
//         numerics::rounding<double> rhd( numerics::rounding<double>::rnd_dw() );
//      let::assign(mn,b);
//       }
//       if( mx < mn ) std::swap(mn,mx);
//      return Interval<double>(mn,mx);
//     }
//   };

}
#endif //numerix_kernel_hpp