flower/rational.cc

   1 /*
   2   rational.cc -- implement Rational
   3
   4   source file of the Flower Library
   5
   6   (c) 1997--2006 Han-Wen Nienhuys <hanwen@xs4all.nl>
   7 */
   8
   9 #include "rational.hh"
  10
  11 #include <cmath>
  12 #include <cassert>
  13 #include <cstdlib>
  14 using namespace std;
  15
  16 #include "string-convert.hh"
  17 #include "libc-extension.hh"
  18
  19 Rational::operator double () const
  20 {
  21   if (sign_ == -1 || sign_ == 1 || sign_ == 0)
  22     return ((double)sign_) * num_ / den_;
  23   if (sign_ == -2)
  24     return -HUGE_VAL;
  25   else if (sign_ == 2)
  26     return HUGE_VAL;
  27   else
  28     assert (false);
  29
  30   return 0.0;
  31 }
  32
  33
  34 #ifdef STREAM_SUPPORT
  35 ostream &
  36 operator << (ostream &o, Rational r)
  37 {
  38   o << r.string ();
  39   return o;
  40 }
  41 #endif
  42
  43 Rational
  44 Rational::abs () const
  45 {
  46   return Rational (num_, den_);
  47 }
  48
  49 Rational
  50 Rational::trunc_rat () const
  51 {
  52   return Rational (num_ - (num_ % den_), den_);
  53 }
  54
  55 Rational::Rational ()
  56 {
  57   sign_ = 0;
  58   num_ = den_ = 1;
  59 }
  60
  61 Rational::Rational (int n, int d)
  62 {
  63   sign_ = ::sign (n) * ::sign (d);
  64   num_ = ::abs (n);
  65   den_ = ::abs (d);
  66   normalize ();
  67 }
  68
  69 Rational::Rational (int n)
  70 {
  71   sign_ = ::sign (n);
  72   num_ = ::abs (n);
  73   den_ = 1;
  74 }
  75
  76
  77 /*
  78   We can actually do a little better. See Knuth 4.5.2
  79  */
  80 static inline
  81 int gcd (int a, int b)
  82 {
  83   int t;
  84   while ((t = a % b))
  85     {
  86       a = b;
  87       b = t;
  88     }
  89   return b;
  90 }
  91
  92 void
  93 Rational::set_infinite (int s)
  94 {
  95   sign_ = ::sign (s) * 2;
  96 }
  97
  98 Rational
  99 Rational::operator - () const
 100 {
 101   Rational r (*this);
 102   r.negate ();
 103   return r;
 104 }
 105
 106 Rational
 107 Rational::div_rat (Rational div) const
 108 {
 109   Rational r (*this);
 110   r /= div;
 111   return r.trunc_rat ();
 112 }
 113
 114 Rational
 115 Rational::mod_rat (Rational div) const
 116 {
 117   Rational r (*this);
 118   r = (r / div - r.div_rat (div)) * div;
 119   return r;
 120 }
 121
 122 void
 123 Rational::normalize ()
 124 {
 125   if (!sign_)
 126     {
 127       den_ = 1;
 128       num_ = 0;
 129     }
 130   else if (!den_)
 131     {
 132       sign_ = 2;
 133       num_ = 1;
 134     }
 135   else if (!num_)
 136     {
 137       sign_ = 0;
 138       den_ = 1;
 139     }
 140   else
 141     {
 142       int g = gcd (num_, den_);
 143
 144       num_ /= g;
 145       den_ /= g;
 146     }
 147 }
 148 int
 149 Rational::sign () const
 150 {
 151   return ::sign (sign_);
 152 }
 153
 154 int
 155 Rational::compare (Rational const &r, Rational const &s)
 156 {
 157   if (r.sign_ < s.sign_)
 158     return -1;
 159   else if (r.sign_ > s.sign_)
 160     return 1;
 161   else if (r.is_infinity ())
 162     return 0;
 163   else if (r.sign_ == 0)
 164     return 0;
 165   return r.sign_ * ::sign (int (r.num_ * s.den_) - int (s.num_ * r.den_));
 166 }
 167
 168 int
 169 compare (Rational const &r, Rational const &s)
 170 {
 171   return Rational::compare (r, s);
 172 }
 173
 174 Rational &
 175 Rational::operator %= (Rational r)
 176 {
 177   *this = mod_rat (r);
 178   return *this;
 179 }
 180
 181 Rational &
 182 Rational::operator += (Rational r)
 183 {
 184   if (is_infinity ())
 185     ;
 186   else if (r.is_infinity ())
 187     *this = r;
 188   else
 189     {
 190       int lcm = (den_ / gcd (r.den_, den_)) * r.den_;
 191       int n = sign_ * num_ * (lcm / den_) + r.sign_ * r.num_ * (lcm / r.den_);
 192       int d = lcm;
 193       sign_ = ::sign (n) * ::sign (d);
 194       num_ = ::abs (n);
 195       den_ = ::abs (d);
 196       normalize ();
 197     }
 198   return *this;
 199 }
 200
 201 /*
 202   copied from libg++ 2.8.0
 203 */
 204 Rational::Rational (double x)
 205 {
 206   if (x != 0.0)
 207     {
 208       sign_ = ::sign (x);
 209       x *= sign_;
 210
 211       int expt;
 212       double mantissa = frexp (x, &expt);
 213
 214       const int FACT = 1 << 20;
 215
 216       /*
 217         Thanks to Afie for this too simple  idea.
 218
 219         do not blindly substitute by libg++ code, since that uses
 220         arbitrary-size integers.  The rationals would overflow too
 221         easily.
 222       */
 223
 224       num_ = (unsigned int) (mantissa * FACT);
 225       den_ = (unsigned int) FACT;
 226       normalize ();
 227       if (expt < 0)
 228         den_ <<= -expt;
 229       else
 230         num_ <<= expt;
 231       normalize ();
 232     }
 233   else
 234     {
 235       num_ = 0;
 236       den_ = 1;
 237       sign_ = 0;
 238       normalize ();
 239     }
 240 }
 241
 242 void
 243 Rational::invert ()
 244 {
 245   int r (num_);
 246   num_ = den_;
 247   den_ = r;
 248 }
 249
 250 Rational &
 251 Rational::operator *= (Rational r)
 252 {
 253   sign_ *= ::sign (r.sign_);
 254   if (r.is_infinity ())
 255     {
 256       sign_ = sign () * 2;
 257       goto exit_func;
 258     }
 259
 260   num_ *= r.num_;
 261   den_ *= r.den_;
 262
 263   normalize ();
 264  exit_func:
 265   return *this;
 266 }
 267
 268 Rational &
 269 Rational::operator /= (Rational r)
 270 {
 271   r.invert ();
 272   return (*this *= r);
 273 }
 274
 275 void
 276 Rational::negate ()
 277 {
 278   sign_ *= -1;
 279 }
 280
 281 Rational &
 282 Rational::operator -= (Rational r)
 283 {
 284   r.negate ();
 285   return (*this += r);
 286 }
 287
 288 string
 289 Rational::to_string () const
 290 {
 291   if (is_infinity ())
 292     {
 293       string s (sign_ > 0 ? "" : "-");
 294       return string (s + "infinity");
 295     }
 296
 297   string s = ::to_string (num ());
 298   if (den () != 1 && num ())
 299     s += "/" + ::to_string (den ());
 300   return s;
 301 }
 302
 303 int
 304 Rational::to_int () const
 305 {
 306   return (int) num () / den ();
 307 }
 308
 309 int
 310 sign (Rational r)
 311 {
 312   return r.sign ();
 313 }
 314
 315 bool
 316 Rational::is_infinity () const
 317 {
 318   return sign_ == 2 || sign_ == -2;
 319 }