core/vnl/vnl_rational.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_rational.h
00002 #ifndef vnl_rational_h_
00003 #define vnl_rational_h_
00004 //:
00005 // \file
00006 // \brief High-precision rational numbers
00007 //
00008 // The  vnl_rational  class  provides  high-precision rational numbers and
00009 // arithmetic, using the built-in type long, for the numerator and denominator.
00010 // Implicit conversion to the system defined types short, int, long, float, and
00011 // double is supported by  overloaded  operator member functions.  Although the
00012 // rational class makes judicious use of inline  functions and  deals only with
00013 // integral values, the user  is warned that  the rational  integer  arithmetic
00014 // class is still considerably slower than the built-in  integer data types. If
00015 // the range  of values  anticipated will  fit into a  built-in  type, use that
00016 // instead.
00017 //
00018 // In  addition  to  the  original  COOL Rational class, vnl_rational is able to
00019 // represent plus and minus infinity.  An  other  interesting  addition  is  the
00020 // possibility  to construct a rational from a double.  This allows for lossless
00021 // conversion from e.g. double 1.0/3.0 to the rational number 1/3, hence no more
00022 // rounding errors.  This is implemented with continued fraction approximations.
00023 //
00024 // \author
00025 // Copyright (C) 1991 Texas Instruments Incorporated.
00026 //
00027 // Permission is granted to any individual or institution to use, copy, modify,
00028 // and distribute this software, provided that this complete copyright and
00029 // permission notice is maintained, intact, in all copies and supporting
00030 // documentation.
00031 //
00032 // Texas Instruments Incorporated provides this software "as is" without
00033 // express or implied warranty.
00034 //
00035 // \verbatim
00036 // Modifications
00037 //  Peter Vanroose, 13 July 2001: Added continued fraction cnstrctr from double
00038 //  Peter Vanroose, 10 July 2001: corrected operator%=()
00039 //  Peter Vanroose, 10 July 2001: corrected ceil() and floor() for negative args
00040 //  Peter Vanroose, 10 July 2001: extended operability range of += by using gcd
00041 //  Peter Vanroose, 10 July 2001: added abs().
00042 //  Peter Vanroose, 10 July 2001: removed state data member and added Inf repres
00043 //  Peter Vanroose,  9 July 2001: ported to vnl from COOL
00044 //  Peter Vanroose, 11 June 2009: made "*" and "/" robust against int overflow
00045 //                                (actually a full re-implementation, using gcd)
00046 // \endverbatim
00047 
00048 #include <vcl_iostream.h>
00049 #include <vcl_cassert.h>
00050 
00051 //: High-precision rational numbers
00052 //
00053 // The  vnl_rational  class  provides  high-precision rational numbers and
00054 // arithmetic, using the built-in type long, for the numerator and denominator.
00055 // Implicit conversion to the system defined types short, int, long, float, and
00056 // double is supported by  overloaded  operator member functions.  Although the
00057 // rational class makes judicious use of inline  functions and  deals only with
00058 // integral values, the user  is warned that  the rational  integer  arithmetic
00059 // class is still considerably slower than the built-in  integer data types. If
00060 // the range  of values  anticipated will  fit into a  built-in  type, use that
00061 // instead.
00062 //
00063 // In  addition  to  the  original  COOL Rational class, vnl_rational is able to
00064 // represent plus and minus infinity.  An  other  interesting  addition  is  the
00065 // possibility  to construct a rational from a double.  This allows for lossless
00066 // conversion from e.g. double 1.0/3.0 to the rational number 1/3, hence no more
00067 // rounding errors.  This is implemented with continued fraction approximations.
00068 //
00069 class vnl_rational
00070 {
00071   long num_; //!< Numerator portion
00072   long den_; //!< Denominator portion
00073 
00074  public:
00075   //: Creates a rational with given numerator and denominator.
00076   //  Default constructor gives 0.
00077   //  Also serves as automatic cast from long to vnl_rational.
00078   //  The only input which is not allowed is (0,0);
00079   //  the denominator is allowed to be 0, to represent +Inf or -Inf.
00080   inline vnl_rational(long num = 0L, long den = 1L)
00081     : num_(num), den_(den) { assert(num!=0||den!=0); normalize(); }
00082   //: Creates a rational with given numerator and denominator.
00083   //  Note these are not automatic type conversions because of a bug
00084   //  in the Borland compiler.  Since these just convert their
00085   //  arguments to long anyway, there is no harm in letting
00086   //  the long overload be used for automatic conversions.
00087   explicit inline vnl_rational(int num, int den = 1)
00088     : num_(num), den_(den) { assert(num!=0||den!=0); normalize(); }
00089   explicit inline vnl_rational(unsigned int num, unsigned int den = 1)
00090     : num_((long)num), den_((long)den) { assert(num!=0||den!=0); normalize(); }
00091   //: Creates a rational from a double.
00092   //  This is done by computing the continued fraction approximation for d.
00093   //  Note that this is explicitly \e not an automatic type conversion.
00094   explicit vnl_rational(double d);
00095   //  Copy constructor
00096   inline vnl_rational(vnl_rational const& from)
00097     : num_(from.numerator()), den_(from.denominator()) {}
00098   //  Destructor
00099   inline ~vnl_rational() {}
00100   //  Assignment: overwrite an existing vnl_rational
00101   inline void set(long num, long den) { assert(num!=0||den!=0); num_=num; den_=den; normalize(); }
00102 
00103   //: Return the numerator of the (simplified) rational number representation
00104   inline long numerator() const { return num_; }
00105   //: Return the denominator of the (simplified) rational number representation
00106   inline long denominator() const { return den_; }
00107 
00108   //: Copies the contents and state of rhs rational over to the lhs
00109   inline vnl_rational& operator=(vnl_rational const& rhs) {
00110     num_ = rhs.numerator(); den_ = rhs.denominator(); return *this; }
00111 
00112   //: Returns true if the two rationals have the same representation
00113   inline bool operator==(vnl_rational const& rhs) const {
00114     return num_ == rhs.numerator() && den_ == rhs.denominator(); }
00115   inline bool operator!=(vnl_rational const& rhs) const { return !operator==(rhs); }
00116   inline bool operator==(long rhs) const { return num_ == rhs && den_ == 1; }
00117   inline bool operator!=(long rhs) const { return !operator==(rhs); }
00118   inline bool operator==(int rhs) const { return num_ == rhs && den_ == 1; }
00119   inline bool operator!=(int rhs) const { return !operator==(rhs); }
00120 
00121   //: Unary minus - returns the negation of the current rational.
00122   inline vnl_rational operator-() const { return vnl_rational(-num_, den_); }
00123   //: Unary plus - returns the current rational.
00124   inline vnl_rational operator+() const { return *this; }
00125   //: Unary not - returns true if rational is equal to zero.
00126   inline bool operator!() const { return num_ == 0L; }
00127   //: Returns the absolute value of the current rational.
00128   inline vnl_rational abs() const { return vnl_rational(num_<0?-num_:num_, den_); }
00129   //: Replaces rational with 1/rational and returns it.
00130   //  Inverting 0 gives +Inf, inverting +-Inf gives 0.
00131   vnl_rational& invert() {
00132     long t = num_; num_ = den_; den_ = t; normalize(); return *this; }
00133 
00134   //: Plus/assign: replace lhs by lhs + rhs
00135   //  Note that +Inf + -Inf and -Inf + +Inf are undefined.
00136   inline vnl_rational& operator+=(vnl_rational const& r) {
00137     if (den_ == r.denominator()) num_ += r.numerator();
00138     else { long c = vnl_rational::gcd(den_,r.denominator()); if (c==0) c=1;
00139            num_ = num_*(r.denominator()/c) + (den_/c)*r.numerator();
00140            den_ *= r.denominator()/c; }
00141     assert(num_!=0 || den_ != 0); // +Inf + -Inf is undefined
00142     normalize(); return *this;
00143   }
00144   inline vnl_rational& operator+=(long r) { num_ += den_*r; return *this; }
00145   //: Minus/assign: replace lhs by lhs - rhs
00146   //  Note that +Inf - +Inf and -Inf - -Inf are undefined.
00147   inline vnl_rational& operator-=(vnl_rational const& r) {
00148     if (den_ == r.denominator()) num_ -= r.num_;
00149     else { long c = vnl_rational::gcd(den_,r.denominator()); if (c==0) c=1;
00150            num_ = num_*(r.denominator()/c) - (den_/c)*r.numerator();
00151            den_ *= r.denominator()/c; }
00152     assert(num_!=0 || den_ != 0); // +Inf - +Inf is undefined
00153     normalize(); return *this;
00154   }
00155   inline vnl_rational& operator-=(long r) { num_ -= den_*r; return *this; }
00156   //: Multiply/assign: replace lhs by lhs * rhs
00157   //  Note that 0 * Inf and Inf * 0 are undefined.
00158   //  Also note that there could be integer overflow during this calculation!
00159   //  In that case, an approximate result will be returned.
00160   vnl_rational& operator*=(vnl_rational const& r);
00161   //: Multiply/assign: replace lhs by lhs * rhs
00162   //  Note that there could be integer overflow during this calculation!
00163   //  In that case, an approximate result will be returned.
00164   vnl_rational& operator*=(long r);
00165   //: Divide/assign: replace lhs by lhs / rhs
00166   //  Note that 0 / 0 and Inf / Inf are undefined.
00167   //  Also note that there could be integer overflow during this calculation!
00168   //  In that case, an approximate result will be returned.
00169   vnl_rational& operator/=(vnl_rational const& r);
00170   //: Divide/assign: replace lhs by lhs / rhs
00171   //  Note that 0 / 0 is undefined.
00172   //  Also note that there could be integer overflow during this calculation!
00173   //  In that case, an approximate result will be returned.
00174   vnl_rational& operator/=(long r);
00175   //: Modulus/assign: replace lhs by lhs % rhs
00176   //  Note that r % Inf is r, and that r % 0 and Inf % r are undefined.
00177   inline vnl_rational& operator%=(vnl_rational const& r) {
00178     assert(r.numerator() != 0);
00179     if (den_ == r.denominator()) num_ %= r.numerator();
00180     else { long c = vnl_rational::gcd(den_,r.denominator()); if (c==0) c=1;
00181            num_ *= r.denominator()/c;
00182            num_ %= (den_/c)*r.numerator();
00183            den_ *= r.denominator()/c; }
00184     normalize(); return *this;
00185   }
00186   inline vnl_rational& operator%=(long r){assert(r);num_%=den_*r;normalize();return *this;}
00187 
00188   //: Pre-increment (++r).  No-op when +-Inf.
00189   inline vnl_rational& operator++() { num_ += den_; return *this; }
00190   //: Pre-decrement (--r).  No-op when +-Inf.
00191   inline vnl_rational& operator--() { num_ -= den_; return *this; }
00192   //: Post-increment (r++).  No-op when +-Inf.
00193   inline vnl_rational operator++(int){vnl_rational b=*this;num_+=den_;return b;}
00194   //: Post-decrement (r--).  No-op when +-Inf.
00195   inline vnl_rational operator--(int){vnl_rational b=*this;num_-=den_;return b;}
00196 
00197   inline bool operator<(vnl_rational const& rhs) const {
00198     if (den_ == rhs.denominator())   // If same denominator
00199       return num_ < rhs.numerator(); // includes the case -Inf < +Inf
00200     // note that denominator is always >= 0:
00201     else
00202       return num_ * rhs.denominator() < den_ * rhs.numerator();
00203   }
00204   inline bool operator>(vnl_rational const& r) const { return r < *this; }
00205   inline bool operator<=(vnl_rational const& r) const { return !operator>(r); }
00206   inline bool operator>=(vnl_rational const& r) const { return !operator<(r); }
00207   inline bool operator<(long r) const { return num_ < den_ * r; }
00208   inline bool operator>(long r) const { return num_ > den_ * r; }
00209   inline bool operator<=(long r) const { return !operator>(r); }
00210   inline bool operator>=(long r) const { return !operator<(r); }
00211   inline bool operator<(int r) const { return num_ < den_ * r; }
00212   inline bool operator>(int r) const { return num_ > den_ * r; }
00213   inline bool operator<=(int r) const { return !operator>(r); }
00214   inline bool operator>=(int r) const { return !operator<(r); }
00215   inline bool operator<(double r) const { return num_ < den_ * r; }
00216   inline bool operator>(double r) const { return num_ > den_ * r; }
00217   inline bool operator<=(double r) const { return !operator>(r); }
00218   inline bool operator>=(double r) const { return !operator<(r); }
00219 
00220   //: Converts rational value to integer by truncating towards zero.
00221   inline long truncate() const { assert(den_ != 0);  return num_/den_; }
00222   //: Converts rational value to integer by truncating towards negative infinity.
00223   inline long floor() const { long t = truncate();
00224     return num_<0L && (num_%den_) != 0 ? t-1 : t; }
00225   //: Converts rational value to integer by truncating towards positive infinity.
00226   inline long ceil() const { long t = truncate();
00227     return num_>0L && (num_%den_) != 0 ? t+1 : t; }
00228   //: Rounds rational to nearest integer.
00229   inline long round() const { long t = truncate();
00230     if (num_ < 0) return ((-num_)%den_) >= 0.5*den_ ? t-1 : t;
00231     else          return   (num_ %den_) >= 0.5*den_ ? t+1 : t;
00232   }
00233 
00234   // Implicit conversions
00235   inline operator short() {
00236     long t = truncate(); short r = (short)t;
00237     assert(r == t); // abort on underflow or overflow
00238     return r;
00239   }
00240   inline operator int() {
00241     long t = truncate(); int r = (int)t;
00242     assert(r == t); // abort on underflow or overflow
00243     return r;
00244   }
00245   inline operator long() const { return truncate(); }
00246   inline operator long() { return truncate(); }
00247   inline operator float() const { return ((float)num_)/((float)den_); }
00248   inline operator float() { return ((float)num_)/((float)den_); }
00249   inline operator double() const { return ((double)num_)/((double)den_); }
00250   inline operator double() { return ((double)num_)/((double)den_); }
00251 
00252   //: Calculate greatest common divisor of two integers.
00253   //  Used to simplify rational number.
00254   static inline long gcd (long l1, long l2) {
00255     while (l2!=0) { long t = l2; l2 = l1 % l2; l1 = t; }
00256     return l1<0 ? (-l1) : l1;
00257   }
00258 
00259  private:
00260   //: Private function to normalize numerator/denominator of rational number.
00261   //  If num_ and den_ are both nonzero, their gcd is made 1 and den_ made positive.
00262   //  Otherwise, the nonzero den_ is set to 1 or the nonzero num_ to +1 or -1.
00263   inline void normalize() {
00264     if (num_ == 0) { den_ = 1; return; } // zero
00265     if (den_ == 0) { num_ = (num_>0) ? 1 : -1; return; } // +-Inf
00266     if (num_ != 1 && num_ != -1 && den_ != 1) {
00267       long common = vnl_rational::gcd(num_, den_);
00268       if (common != 1) { num_ /= common; den_ /= common; }
00269     }
00270     // if negative, put sign in numerator:
00271     if (den_ < 0) { num_ *= -1; den_ *= -1; }
00272   }
00273 };
00274 
00275 //: formatted output
00276 // \relatesalso vnl_rational
00277 inline vcl_ostream& operator<<(vcl_ostream& s, vnl_rational const& r)
00278 {
00279   return s << r.numerator() << '/' << r.denominator();
00280 }
00281 
00282 //: simple input
00283 // \relatesalso vnl_rational
00284 inline vcl_istream& operator>>(vcl_istream& s, vnl_rational& r)
00285 {
00286   long n, d; s >> n >> d;
00287   r.set(n,d); return s;
00288 }
00289 
00290 //: Returns the sum of two rational numbers.
00291 // \relatesalso vnl_rational
00292 inline vnl_rational operator+(vnl_rational const& r1, vnl_rational const& r2)
00293 {
00294   vnl_rational result(r1); return result += r2;
00295 }
00296 
00297 inline vnl_rational operator+(vnl_rational const& r1, long r2)
00298 {
00299   vnl_rational result(r1); return result += r2;
00300 }
00301 
00302 inline vnl_rational operator+(vnl_rational const& r1, int r2)
00303 {
00304   vnl_rational result(r1); return result += (long)r2;
00305 }
00306 
00307 inline vnl_rational operator+(long r2, vnl_rational const& r1)
00308 {
00309   vnl_rational result(r1); return result += r2;
00310 }
00311 
00312 inline vnl_rational operator+(int r2, vnl_rational const& r1)
00313 {
00314   vnl_rational result(r1); return result += (long)r2;
00315 }
00316 
00317 //: Returns the difference of two rational numbers.
00318 // \relatesalso vnl_rational
00319 inline vnl_rational operator-(vnl_rational const& r1, vnl_rational const& r2)
00320 {
00321   vnl_rational result(r1); return result -= r2;
00322 }
00323 
00324 inline vnl_rational operator-(vnl_rational const& r1, long r2)
00325 {
00326   vnl_rational result(r1); return result -= r2;
00327 }
00328 
00329 inline vnl_rational operator-(vnl_rational const& r1, int r2)
00330 {
00331   vnl_rational result(r1); return result -= (long)r2;
00332 }
00333 
00334 inline vnl_rational operator-(long r2, vnl_rational const& r1)
00335 {
00336   vnl_rational result(-r1); return result += r2;
00337 }
00338 
00339 inline vnl_rational operator-(int r2, vnl_rational const& r1)
00340 {
00341   vnl_rational result(-r1); return result += (long)r2;
00342 }
00343 
00344 //: Returns the product of two rational numbers.
00345 // \relatesalso vnl_rational
00346 inline vnl_rational operator*(vnl_rational const& r1, vnl_rational const& r2)
00347 {
00348   vnl_rational result(r1); return result *= r2;
00349 }
00350 
00351 inline vnl_rational operator*(vnl_rational const& r1, long r2)
00352 {
00353   vnl_rational result(r1); return result *= r2;
00354 }
00355 
00356 inline vnl_rational operator*(vnl_rational const& r1, int r2)
00357 {
00358   vnl_rational result(r1); return result *= (long)r2;
00359 }
00360 
00361 inline vnl_rational operator*(long r2, vnl_rational const& r1)
00362 {
00363   vnl_rational result(r1); return result *= r2;
00364 }
00365 
00366 inline vnl_rational operator*(int r2, vnl_rational const& r1)
00367 {
00368   vnl_rational result(r1); return result *= (long)r2;
00369 }
00370 
00371 //: Returns the quotient of two rational numbers.
00372 // \relatesalso vnl_rational
00373 inline vnl_rational operator/(vnl_rational const& r1, vnl_rational const& r2)
00374 {
00375   vnl_rational result(r1); return result /= r2;
00376 }
00377 
00378 inline vnl_rational operator/(vnl_rational const& r1, long r2)
00379 {
00380   vnl_rational result(r1); return result /= r2;
00381 }
00382 
00383 inline vnl_rational operator/(vnl_rational const& r1, int r2)
00384 {
00385   vnl_rational result(r1); return result /= (long)r2;
00386 }
00387 
00388 inline vnl_rational operator/(long r1, vnl_rational const& r2)
00389 {
00390   vnl_rational result(r1); return result /= r2;
00391 }
00392 
00393 inline vnl_rational operator/(int r1, vnl_rational const& r2)
00394 {
00395   vnl_rational result((long)r1); return result /= r2;
00396 }
00397 
00398 //: Returns the remainder of r1 divided by r2.
00399 // \relatesalso vnl_rational
00400 inline vnl_rational operator%(vnl_rational const& r1, vnl_rational const& r2)
00401 {
00402   vnl_rational result(r1); return result %= r2;
00403 }
00404 
00405 inline vnl_rational operator%(vnl_rational const& r1, long r2)
00406 {
00407   vnl_rational result(r1); return result %= r2;
00408 }
00409 
00410 inline vnl_rational operator%(vnl_rational const& r1, int r2)
00411 {
00412   vnl_rational result(r1); return result %= (long)r2;
00413 }
00414 
00415 inline vnl_rational operator%(long r1, vnl_rational const& r2)
00416 {
00417   vnl_rational result(r1); return result %= r2;
00418 }
00419 
00420 inline vnl_rational operator%(int r1, vnl_rational const& r2)
00421 {
00422   vnl_rational result((long)r1); return result %= r2;
00423 }
00424 
00425 inline bool operator==(int  r1, vnl_rational const& r2) { return r2==r1; }
00426 inline bool operator==(long r1, vnl_rational const& r2) { return r2==r1; }
00427 inline bool operator!=(int  r1, vnl_rational const& r2) { return r2!=r1; }
00428 inline bool operator!=(long r1, vnl_rational const& r2) { return r2!=r1; }
00429 inline bool operator< (int  r1, vnl_rational const& r2) { return r2> r1; }
00430 inline bool operator< (long r1, vnl_rational const& r2) { return r2> r1; }
00431 inline bool operator> (int  r1, vnl_rational const& r2) { return r2< r1; }
00432 inline bool operator> (long r1, vnl_rational const& r2) { return r2< r1; }
00433 inline bool operator<=(int  r1, vnl_rational const& r2) { return r2>=r1; }
00434 inline bool operator<=(long r1, vnl_rational const& r2) { return r2>=r1; }
00435 inline bool operator>=(int  r1, vnl_rational const& r2) { return r2<=r1; }
00436 inline bool operator>=(long r1, vnl_rational const& r2) { return r2<=r1; }
00437 
00438 inline long truncate(vnl_rational const& r) { return r.truncate(); }
00439 inline long floor(vnl_rational const& r) { return r.floor(); }
00440 inline long ceil(vnl_rational const& r) { return r.ceil(); }
00441 inline long round(vnl_rational const& r) { return r.round(); }
00442 
00443 inline vnl_rational vnl_math_abs(vnl_rational const& x) { return x<0L ? -x : x; }
00444 inline vnl_rational vnl_math_squared_magnitude(vnl_rational const& x) { return x*x; }
00445 inline vnl_rational vnl_math_sqr(vnl_rational const& x) { return x*x; }
00446 inline bool vnl_math_isnan(vnl_rational const& ){return false;}
00447 inline bool vnl_math_isfinite(vnl_rational const& x){return x.denominator() != 0L;}
00448 
00449 #endif // vnl_rational_h_