00001 #ifndef numerix_kernel_hpp
00002 #define numerix_kernel_hpp
00003
00004 #include <math.h>
00005 #include <basix/basix.hpp>
00006
00007 #include <numerix/integer.hpp>
00008 #include <numerix/integer_meta.hpp>
00009 #include <numerix/rational.hpp>
00010 #include <numerix/rational_meta.hpp>
00011 #include <numerix/floating.hpp>
00012 #include <numerix/floating_meta.hpp>
00013 #include <numerix/ball.hpp>
00014 #include <numerix/interval.hpp>
00015
00016
00017
00018
00019
00020 namespace mmx {
00021 template<class OSTREAM> inline void
00022 print(OSTREAM& os, const integer& x) { os << as_charp(as_string(x)); }
00023
00024 template<class OSTREAM> inline void
00025 print(OSTREAM& os, const rational& x) {
00026 os << as_charp(as_string(numerator(x)));
00027 if (denominator(x) != 1){
00028 os <<"/";
00029 os << as_charp(as_string(denominator(x)));
00030 }
00031 }
00032
00033
00034 template<class OSTREAM, class V> inline void
00035 print(OSTREAM& os, const floating<V>& x) { os << as_charp(as_string(x)); }
00036
00037
00038
00039 template<class X> struct scalar_set {
00040 typedef X T;
00041 scalar_set(){}
00042 scalar_set(const scalar_set<X>& t) {}
00043 scalar_set(const X& x) {}
00044 };
00045
00046 template<class X>
00047 bool operator== (const scalar_set<X>& r1, const scalar_set<X>& r2) { return true; }
00048 template<class X>
00049 bool operator!= (const scalar_set<X>& r1, const scalar_set<X>& r2) { return !(r1==r2); }
00050 template<class C> inline int hash(const scalar_set<C>& R){ return 1;}
00051 template<class C> inline int exact_hash(const scalar_set<C>& R){ return 1;}
00052 template<class C> inline bool
00053 exact_eq(const scalar_set<C>& R1, const scalar_set<C>& R2) {
00054 return R1==R2;
00055 }
00056
00057 template<class C> inline scalar_set<C> make_set(const C& c) {
00058 return scalar_set<C>();
00059 }
00060
00061 template<typename U,typename T>
00062 struct as_helper<scalar_set<U>,scalar_set<T> > {
00063 static inline
00064 scalar_set<U> cv (const scalar_set<T>& x) { return scalar_set<U>(); }
00065 };
00066
00067 #define set_of(T) scalar_set< T > ()
00068
00069 #define DECLARE_SCALAR_SET(X,T) \
00070 inline syntactic flatten(const scalar_set<T>& Z){return syntactic(X);} \
00071
00072 #define bigfloat floating<>
00073
00074 DECLARE_SCALAR_SET("DD",double);
00075 DECLARE_SCALAR_SET("ZZ",integer);
00076 DECLARE_SCALAR_SET("QQ",rational);
00077 DECLARE_SCALAR_SET("RR",bigfloat);
00078 DECLARE_SCALAR_SET("CD",complex< double >);
00079 DECLARE_SCALAR_SET("CC",complex< bigfloat >);
00080
00081 #ifndef WITH_AS
00082 #define WITH_AS
00083 template<typename T,typename F>
00084 struct as_helper { static inline T cv (const F& x) { return x; } };
00085
00086 template<typename T,typename F> inline T as (const F& x) { return as_helper<T,F>::cv (x); }
00087 #endif
00088
00089
00090 struct Numerix {
00091
00092 typedef ::mmx::integer integer;
00093 typedef ::mmx::rational rational;
00094 typedef ::mmx::floating<> floating;
00095 typedef double ieee;
00096
00097 typedef ::mmx::interval<rational> interval_rational;
00098 };
00099
00100
00101 namespace texp {
00102 template<class R> struct kernelof_;
00103
00104 template<> struct kernelof_<Numerix::integer> { typedef Numerix T; };
00105 template<> struct kernelof_<Numerix::rational> { typedef Numerix T; };
00106 template<> struct kernelof_<Numerix::floating> { typedef Numerix T; };
00107 }
00108
00109
00110
00111 template<class T> inline
00112 bool singleton(const interval<T>& x ) { return lower(x) == upper(x); };
00113
00114 namespace meta {
00115 template<class C> struct kernelof_;
00116 template<> struct kernelof_< integer > { typedef Numerix T; };
00117 template<> struct kernelof_< rational > { typedef Numerix T; };
00118 template<> struct kernelof_< floating<> > { typedef Numerix T; };
00119
00120
00121
00122
00123
00124
00125 }
00126
00127
00128
00129 namespace let {
00130 template<class A,class B> void assign(A& a, const B&b);
00131
00132 inline void
00133 assign(rational& q, char* s) {q= rational(integer(s)); }
00134
00135
00136 template<> inline void
00137 assign(double& x, const rational& q) {x = as_double(q);}
00138
00139
00140
00141 template<> inline void
00142 assign(floating<>& r, const int& n) {r = floating<>(n);}
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 }
00160
00161
00162 inline bool less_prec(const rational& q, unsigned prec)
00163 {
00164 integer n;
00165 if(q<0)
00166 n=-(numerator(q)<<prec);
00167 else
00168 n=numerator(q)<<prec;
00169 return (n<denominator(q));
00170 }
00171
00172 inline rational nearest(const floating<>& x, int prec)
00173 {
00174 integer er = integer(1)<<prec;
00175
00176 floating<> r=(x>0?x:-x), f = floor(r);
00177 integer a=as_integer(f);
00178
00179 integer n0=1, n1= a, tn=a;
00180 integer d0=0, d1= 1, td=1;
00181
00182 for(;d0*d1<er && r!=f;)
00183 {
00184 r=1/(r-f); f = floor(r);
00185 a=as_integer(f);
00186 n1*=a;n1+=n0;n0=tn;tn=n1;
00187 d1*=a;d1+=d0;d0=td;td=d1;
00188 }
00189 return (x>0? rational(n1)/d1: rational(-n1)/d1);
00190 }
00191 template<> struct as_helper<rational,floating<> >
00192 {
00193 static inline
00194 rational cv (const floating<>& x)
00195 {
00196 return nearest(x,precision(x));
00197 }
00198 };
00199
00200 template<> struct as_helper<double,floating<> >
00201 {
00202 static inline double cv (const floating<>& x) { return as_double(x); }
00203 };
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 }
00226 #endif //numerix_kernel_hpp