00001
00002
00003
00004
00005
00006
00007 #ifndef SYNAPS_ARITHM_DIADIC_H
00008 #define SYNAPS_ARITHM_DIADIC_H
00009
00010 #include <cassert>
00011 #include <synaps/init.h>
00012 #include <synaps/base/type.h>
00013 #include <synaps/base/comparison.h>
00014 #include <synaps/arithm/let.h>
00015 #include <synaps/arithm/pow.h>
00016 #include <synaps/arithm/gmp.h>
00017
00018 __BEGIN_NAMESPACE_SYNAPS
00019
00020 template <class NT_>
00021 NT_ shift_2(const NT_& z, int d)
00022 {
00023 if(d>0)
00024 return z*pow(NT_(2),d);
00025 else
00026 return z/pow(NT_(2),-d);
00027 }
00028
00029
00030 template< typename NT_ > class diadic;
00031 namespace let
00032 {
00033 template<typename NT_> void assign(diadic<NT_> & x, int i);
00034 }
00035
00040 template < typename NT_ >
00041 class diadic
00042 {
00043 public:
00044 typedef NT_ NT;
00045
00046 diadic() : num_(0), exp_(0) { }
00047
00048 diadic(const NT& n,int d=0) : num_(n), exp_(d) {normalize();}
00049
00050 template < typename RT>
00051 diadic(const RT& n,int d=0)
00052 {
00053 let::assign(*this,n);
00054 exp_+=d;
00055 }
00056
00057 diadic<NT>& operator+= (const diadic<NT>& r);
00058 diadic<NT>& operator-= (const diadic<NT>& r);
00059 diadic<NT>& operator*= (const diadic<NT>& r);
00060 diadic<NT>& operator/= (const diadic<NT>& r);
00061 diadic<NT>& operator/= (int r);
00062 diadic<NT>& operator+= (const NT& r);
00063 diadic<NT>& operator-= (const NT& r);
00064 diadic<NT>& operator*= (const NT& r);
00065 diadic<NT>& operator/= (const NT& r);
00066
00067 diadic<NT>& normalize();
00068
00069 const NT& numerator() const { return num_; }
00070 NT denominator() const { return shift_2(1,exp_); }
00071 int valuation () const { return exp_; }
00072
00073 public:
00074 NT num_;
00075 int exp_;
00076 };
00077
00078 template < typename NT > inline
00079 diadic<NT>& diadic<NT>::normalize()
00080 {
00081 if( num_==0)
00082 exp_=0;
00083 else
00084 while( (num_%2) ==0 )
00085 {
00086 exp_++;
00087 num_/=2;
00088 }
00089 return *this;
00090 }
00091
00092
00093
00095 template < typename NT > inline
00096 diadic<NT>&
00097 diadic<NT>::operator+= (const diadic<NT>& r)
00098 {
00099 if(r.exp_<exp_)
00100 {
00101 int d=exp_-r.exp_;
00102 num_ = shift_2(num_ ,d);
00103 num_ += r.num_;
00104 exp_ = r.exp_;
00105 }
00106 else if(r.exp_>exp_)
00107 {
00108 num_ += shift_2(r.num_,r.exp_-exp_);
00109 }
00110 else
00111 {
00112 num_ += r.num_;
00113 normalize();
00114 }
00115 return *this;
00116 }
00117
00119 template < typename NT > inline
00120 diadic<NT>&
00121 diadic<NT>::operator-= (const diadic<NT>& r)
00122 {
00123 if(r.exp_<exp_)
00124 {
00125 num_ = shift_2(num_,exp_-r.exp_);
00126 num_ -= r.num_;
00127 exp_ = r.exp_;
00128 }
00129 else if(r.exp_>exp_)
00130 {
00131 num_ -= shift_2(r.num_ ,r.exp_-exp_);
00132 }
00133 else
00134 {
00135 num_ -= r.num_;
00136 normalize();
00137 }
00138 return *this;
00139 }
00140
00142 template < typename NT > inline
00143 diadic<NT>&
00144 diadic<NT>::operator*= (const diadic<NT>& r)
00145 {
00146 num_ *= r.num_;
00147 exp_ += r.exp_;
00148 return *this;
00149 }
00150
00155 template < typename NT > inline
00156 diadic<NT>&
00157 diadic<NT>::operator/= (const diadic<NT>& r)
00158 {
00159 assert( r.num_ != 0 );
00160 exp_ -= r.exp_;
00161 num_ /= r.num_;
00162 normalize();
00163 return *this;
00164 }
00165
00169 template < typename NT > inline
00170 diadic<NT>&
00171 diadic<NT>::operator/= (const NT& r)
00172 {
00173 assert( r != 0 );
00174 NT d=r;
00175 while( (d%2) ==0 )
00176 {
00177 exp_--;
00178 d/=2;
00179 }
00180 num_/=d;
00181 return *this;
00182 }
00183
00187 template < typename NT > inline
00188 diadic<NT>&
00189 diadic<NT>::operator/= (int r)
00190 {
00191 assert( r != 0 );
00192 int d=r;
00193 while( d%2==0 ) { exp_--; d/=2; }
00194 num_/=d;
00195 return *this;
00196 }
00197
00199 template < typename NT > inline
00200 diadic<NT>&
00201 diadic<NT>::operator+= (const NT& r)
00202 {
00203 if(exp_<0)
00204 num_ += shift_2(r,exp_);
00205 else if(exp_>0)
00206 {
00207 num_ = shift_2(num_,exp_);
00208 num_+=r;
00209 }
00210 else
00211 num_+=r;
00212 normalize();
00213 return *this;
00214 }
00215
00217 template < typename NT > inline
00218 diadic<NT>&
00219 diadic<NT>::operator-= (const NT& r)
00220 {
00221 if(exp_<0)
00222 num_ -= shift_2(r,exp_);
00223 else if(exp_>0)
00224 {
00225 num_ = shift_2(num_,exp_);
00226 num_-= r;
00227 }
00228 else
00229 num_-= r;
00230 normalize();
00231 return *this;
00232 }
00233
00235 template < typename NT > inline
00236 diadic<NT>&
00237 diadic<NT>::operator*= (const NT& r)
00238 {
00239 num_ *= r;
00240 normalize();
00241 return *this;
00242 }
00243
00244
00245 template < typename NT > inline
00246 int
00247 quotient_cmp(const diadic<NT>& a, const diadic<NT>& b)
00248 {
00249 int asign = sign(a.num_);
00250 int bsign = sign(b.num_);
00251
00252 if ((asign == 0) && (bsign == 0)) return EQUAL;
00253 if (asign < bsign) return SMALLER;
00254 if (asign > bsign) return LARGER;
00255
00256
00257 NT ta=a.num_, tb=b.num_;
00258 int d=a.exp_-b.exp_;
00259
00260 if(d>0)
00261 ta = shift_2(ta,d);
00262 else
00263 tb = shift_2(tb,-d);
00264
00265 if(ta<tb)
00266 return SMALLER;
00267 else if(tb<ta)
00268 return LARGER;
00269 else
00270 return EQUAL;
00271 }
00272
00273 template < typename NT >
00274 inline
00275 int
00276 compare(const diadic<NT>& x, const diadic<NT>& y)
00277 { return quotient_cmp(x, y); }
00278
00279 template < typename NT >
00280 std::ostream&
00281 operator<<(std::ostream& s, const diadic<NT>& r)
00282 {
00283 s << r.numerator();
00284 if(r.valuation()) s<< "*2^(" << r.valuation()<<")";
00285 return s;
00286 }
00287
00288
00289
00291 template < typename NT > inline
00292 diadic<NT>
00293 operator+(const diadic<NT>& x, const diadic<NT>& y)
00294 {
00295 diadic<NT> z = x;
00296 return z += y;
00297 }
00298
00300 template < typename NT > inline
00301 diadic<NT>
00302 operator-(const diadic<NT>& x, const diadic<NT>& y)
00303 { return (diadic<NT>(x) -= y); }
00304
00305
00307 template < typename NT > inline
00308 diadic<NT>
00309 operator*(const diadic<NT>& x, const diadic<NT>& y)
00310 {
00311 diadic<NT> z = x;
00312 return z *= y;
00313 }
00314
00316 template < typename NT > inline
00317 diadic<NT>
00318 operator/(const diadic<NT>& x, const diadic<NT>& y)
00319 {
00320 diadic<NT> z = x;
00321 return z /= y;
00322 }
00323
00325 template < typename NT > inline
00326 diadic<NT>
00327 operator/(const diadic<NT>& x, int d)
00328 {
00329 diadic<NT> z(x);
00330 return z /= d;
00331 }
00332
00334 template < typename NT > inline
00335 diadic<NT>
00336 operator-(const diadic<NT>& x)
00337 {
00338 return diadic<NT>(-x.num_,x.exp_);
00339 }
00340
00342 template < typename NT > inline
00343 diadic<NT>
00344 operator+(const NT& x, const diadic<NT>& y)
00345 {
00346 diadic<NT> z(x);
00347 return z += y;
00348 }
00349
00351 template < typename NT > inline
00352 diadic<NT>
00353 operator+(const diadic<NT>& x, const NT& y)
00354 {
00355 diadic<NT> z = x;
00356 return z += y;
00357 }
00358
00360 template < typename NT > inline
00361 diadic<NT>
00362 operator-(const NT& x, const diadic<NT>& y)
00363 {
00364 diadic<NT> z(x);
00365 return z -= y;
00366 }
00367
00369 template < typename NT > inline
00370 diadic<NT>
00371 operator-(const diadic<NT>& x, const NT& y)
00372 {
00373 diadic<NT> z = x;
00374 return z -= y;
00375 }
00376
00378 template < typename NT > inline
00379 diadic<NT>
00380 operator*(const NT& x, const diadic<NT>& y)
00381 {
00382 diadic<NT> z(x);
00383 return z *= y;
00384 }
00385
00386 template < typename NT > inline
00387 diadic<NT>
00388 operator*(const diadic<NT>& x, const NT& y)
00389 {
00390 diadic<NT> z = x;
00391 return z *= y;
00392 }
00393
00394
00395 #ifdef SYNAPS_WITH_GMP
00396 QQ operator*(const diadic<ZZ>& x, const QQ& y)
00397 {
00398 return (y*x.numerator()/shift_2(QQ(1),x.valuation()));
00399 }
00400
00401 QQ& operator*=(QQ& y, const diadic<ZZ>& x)
00402 {
00403 y *=x.numerator();
00404 y /=shift_2(QQ(1),x.valuation());
00405 return y;
00406 }
00407 #endif //SYNAPS_WITH_GMP
00408
00410 template < typename NT > inline
00411 diadic<NT>
00412 operator/(const NT& x, const diadic<NT>& y)
00413 {
00414 diadic<NT> z(x) ;
00415 return z /= y;
00416 }
00417
00418 template < typename NT > inline
00419 diadic<NT>
00420 operator/(const diadic<NT>& x, const NT& y)
00421 {
00422 diadic<NT> z = x;
00423 return z /= y;
00424 }
00425
00427 template < typename NT > inline
00428 int
00429 sign(const diadic<NT>& x)
00430 {
00431 return sign(x.num_);
00432 }
00433
00434
00435
00437 template < typename NT > inline
00438 bool
00439 operator==(const diadic<NT>& x, const diadic<NT>& y)
00440 {
00441 return shift_2(x.num_ ,y.exp_) == shift_2(y.num_,x.exp_);
00442 }
00443
00445 template < typename NT > inline
00446 bool operator==(const diadic<NT>& x, const NT& y)
00447 { return shift_2(y, x.exp_) == x.num_; }
00448
00450 template < typename NT > inline
00451 bool
00452 operator==(const NT& x, const diadic<NT>& y)
00453 { return y == x; }
00454
00456 template < typename NT > inline
00457 bool
00458 operator!=(const diadic<NT>& x, const diadic<NT>& y)
00459 { return ! (x == y); }
00460
00462 template < typename NT > inline
00463 bool
00464 operator!=(const diadic<NT>& x, const NT& y)
00465 { return ! (x == y); }
00466
00468 template < typename NT > inline
00469 bool
00470 operator!=(const NT& x, const diadic<NT>& y)
00471 { return ! (x == y); }
00472
00474 template < typename NT > inline
00475 bool
00476 operator<(const diadic<NT>& x, const diadic<NT>& y)
00477 {
00478 return quotient_cmp(x,y) == SMALLER;
00479 }
00480
00482 template < typename NT > inline
00483 bool
00484 operator<(const diadic<NT>& x, const NT& y)
00485 {
00486 return quotient_cmp(x,diadic<NT>(y)) == SMALLER;
00487 }
00488
00490 template < typename NT > inline
00491 bool
00492 operator<(const NT& x, const diadic<NT>& y)
00493 {
00494 return quotient_cmp(diadic<NT>(x),y) == SMALLER;
00495 }
00496
00498 template < typename NT > inline
00499 bool
00500 operator>(const diadic<NT>& x, const diadic<NT>& y)
00501 {
00502 return(y < x);
00503 }
00504
00506 template < typename NT > inline
00507 bool
00508 operator>(const diadic<NT>& x, const NT& y)
00509 {
00510 return y < x;
00511 }
00512
00514 template < typename NT > inline
00515 bool
00516 operator>(const NT& x, const diadic<NT>& y)
00517 { return y < x; }
00518
00520 template < typename NT > inline
00521 bool
00522 operator<=(const diadic<NT>& x, const diadic<NT>& y)
00523 { return ! (y < x); }
00524
00526 template < typename NT > inline
00527 bool
00528 operator<=(const diadic<NT>& x, const NT& y)
00529 { return ! (y < x); }
00530
00531 template < typename NT > inline
00532 bool
00533 operator<=(const NT& x, const diadic<NT>& y)
00534 { return ! (y < x); }
00535
00537 template < typename NT > inline
00538 bool
00539 operator>=(const diadic<NT>& x, const diadic<NT>& y)
00540 { return ! (x < y); }
00541
00543 template < typename NT > inline
00544 bool
00545 operator>=(const diadic<NT>& x, const NT& y)
00546 { return ! (x < y); }
00547
00549 template < typename NT > inline
00550 bool
00551 operator>=(const NT& x, const diadic<NT>& y)
00552 { return ! (x < y); }
00553
00555 template < typename NT > inline
00556 NT denominator(const diadic<NT>& q)
00557 { return q.denominator(); }
00558
00560 template < typename NT > inline
00561 const NT&
00562 numerator(const diadic<NT>& q)
00563 {
00564 return q.numerator();
00565 }
00566
00567 template<class NT,class X>
00568 bool operator< (const diadic<NT>& d, const X& q)
00569 {
00570 NT r;
00571 if(d.valuation()>=0)
00572 {
00573 r= shift_2(d.numerator(),d.valuation());
00574 return (r<q);
00575 }
00576 else
00577 {
00578 r=NT(q);
00579 r= shift_2(r,-d.valuation());
00580 return (d.numerator()<r);
00581 }
00582 }
00583
00584 template<class NT,class X>
00585 bool operator> (const diadic<NT>& d, const X& q)
00586 {
00587 NT r;
00588 if(d.valuation()>=0)
00589 {
00590 r= shift_2(d.numerator(),d.valuation());
00591 return (r>q);
00592 }
00593 else
00594 {
00595 r=NT(q);
00596 r= shift_2(r,-d.valuation());
00597 return (d.numerator()>r);
00598 }
00599 }
00600
00601 template<class NT,class X>
00602 bool operator> (const X& q, const diadic<NT>& d)
00603 {
00604 return (d<q);
00605 }
00606
00607 template<class NT,class X>
00608 bool operator< (const X& q, const diadic<NT>& d)
00609 {
00610 return (d>q);
00611 }
00612
00613
00614
00615
00616 namespace let {
00618 template < typename NT >
00619 double
00620 convert(const diadic<NT>& q, type::As<double> )
00621 {
00622 return (convert(q.numerator(),type::As<double>())/shift_2(1.0,q.valuation()));
00623 }
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637 #ifdef SYNAPS_WITH_GMP
00638 void assign(QQ& r, const diadic<ZZ> & x)
00639 {
00640 assign(r,x.numerator());
00641 for(unsigned i=0;i<x.exp_;i++) r/=2;
00642 }
00643
00644 template<>
00645 void assign(diadic<ZZ> & x, int n)
00646 {
00647 assign(x.num_,n);
00648 x.exp_=0;
00649 x.normalize();
00650 }
00651 #endif //SYNAPS_WITH_GMP
00652
00653 }
00654
00655 __END_NAMESPACE_SYNAPS
00656 #endif // SYNAPS_ARITHM_DIADIC_H
00657
