core/vnl/vnl_rational.cxx
Go to the documentation of this file.
00001 // This is core/vnl/vnl_rational.cxx
00002 #include "vnl_rational.h"
00003 //:
00004 // \file
00005 
00006 #include <vnl/vnl_numeric_traits.h> // for vnl_numeric_traits<long>::maxval
00007 #include <vcl_cassert.h>
00008 
00009 //: Creates a rational from a double.
00010 //  This is done by computing the continued fraction approximation for d.
00011 vnl_rational::vnl_rational(double d)
00012 {
00013   bool sign = d<0;
00014   if (sign) d = -d;
00015 
00016   // Continued fraction approximation of abs(d): recursively determined
00017   long den=0L, num=1L, prev_den=1L, prev_num=0L;
00018 
00019   while (d*num < 1e9 && d*den < 1e9) {
00020     long a = (long)d; // integral part of d
00021     d -= a; // certainly >= 0
00022     long temp = num; num = a*num + prev_num; prev_num = temp;
00023          temp = den; den = a*den + prev_den; prev_den = temp;
00024     if (d < 1e-6) break;
00025     d = 1/d;
00026   }
00027   num_ = num; den_ = den;
00028   if (sign) num_ = -num_;
00029   // no need to normalize() since prev_num and prev_den have guaranteed a gcd=1
00030 }
00031 
00032 //: Multiply/assign: replace lhs by lhs * rhs
00033 //  Note that 0 * Inf and Inf * 0 are undefined.
00034 //  Also note that there could be integer overflow during this calculation!
00035 //  In that case, an approximate result will be returned.
00036 vnl_rational& vnl_rational::operator*=(vnl_rational const& r)
00037 {
00038   assert(num_!=0 || den_ != 0); // 0 * Inf is undefined
00039   long a = vnl_rational::gcd(r.numerator(),den_),
00040        b = vnl_rational::gcd(r.denominator(),num_);
00041   num_ /= b; den_ /= a;
00042   a = r.numerator()/a; b = r.denominator()/b;
00043   // find out whether overflow would occur; in that case, return approximate result
00044   double n = double(a) * double(num_),
00045          d = double(b) * double(den_);
00046   if (n < vnl_numeric_traits<long>::maxval && d < vnl_numeric_traits<long>::maxval)
00047   { num_ *= a; den_ *= b; normalize(); return *this; }
00048   else
00049     return *this = vnl_rational(n/d);
00050 }
00051 
00052 //: Multiply/assign: replace lhs by lhs * rhs
00053 //  Note that there could be integer overflow during this calculation!
00054 //  In that case, an approximate result will be returned.
00055 vnl_rational& vnl_rational::operator*=(long r)
00056 {
00057   long a = vnl_rational::gcd(r,den_);
00058   den_ /= a; r /= a;
00059   // find out whether overflow would occur; in that case, return approximate result
00060   double n = double(r) * double(num_);
00061   if (n < vnl_numeric_traits<long>::maxval)
00062   { num_ *= r; normalize(); return *this; }
00063   else
00064     return *this = vnl_rational(n/double(den_));
00065 }
00066 
00067 //: Divide/assign: replace lhs by lhs / rhs
00068 //  Note that 0 / 0 and Inf / Inf are undefined.
00069 //  Also note that there could be integer overflow during this calculation!
00070 //  In that case, an approximate result will be returned.
00071 vnl_rational& vnl_rational::operator/=(vnl_rational const& r)
00072 {
00073   assert(num_!=0 || den_ != 0); // 0/0, Inf/Inf undefined
00074   long a = vnl_rational::gcd(r.numerator(),num_),
00075        b = vnl_rational::gcd(r.denominator(),den_);
00076   num_ /= a; den_ /= b;
00077   a = r.numerator()/a; b = r.denominator()/b;
00078   // find out whether overflow would occur; in that case, return approximate result
00079   double n = double(b) * double(num_),
00080          d = double(a) * double(den_);
00081   if (n < vnl_numeric_traits<long>::maxval && d < vnl_numeric_traits<long>::maxval)
00082   { num_ *= b; den_ *= a; normalize(); return *this; }
00083   else
00084     return *this = vnl_rational(n/d);
00085 }
00086 
00087 //: Divide/assign: replace lhs by lhs / rhs
00088 //  Note that 0 / 0 is undefined.
00089 //  Also note that there could be integer overflow during this calculation!
00090 //  In that case, an approximate result will be returned.
00091 vnl_rational& vnl_rational::operator/=(long r)
00092 {
00093   assert(num_!=0 || r != 0); // 0/0 undefined
00094   long a = vnl_rational::gcd(r,num_);
00095   num_ /= a; r /= a;
00096   // find out whether overflow would occur; in that case, return approximate result
00097   double d = double(r) * double(den_);
00098   if (d < vnl_numeric_traits<long>::maxval)
00099   { den_ *= r; normalize(); return *this; }
00100   else
00101     return *this = vnl_rational(double(num_)/d);
00102 }