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 static inline
  70 int gcd (int a, int b)
  71 {
  72   int t;
  73   while ((t = a % b))
  74     {
  75       a = b;
  76       b = t;
  77     }
  78   return b;
  79 }
  80
  81 void
  82 Rational::set_infinite (int s)
  83 {
  84   sign_ = ::sign (s) * 2;
  85 }
  86
  87 Rational
  88 Rational::operator - () const
  89 {
  90   Rational r (*this);
  91   r.negate ();
  92   return r;
  93 }
  94
  95 Rational
  96 Rational::div_rat (Rational div) const
  97 {
  98   Rational r (*this);
  99   r /= div;
 100   return r.trunc_rat ();
 101 }
 102
 103 Rational
 104 Rational::mod_rat (Rational div) const
 105 {
 106   Rational r (*this);
 107   r = (r / div - r.div_rat (div)) * div;
 108   return r;
 109 }
 110
 111 void
 112 Rational::normalize ()
 113 {
 114   if (!sign_)
 115     {
 116       den_ = 1;
 117       num_ = 0;
 118     }
 119   else if (!den_)
 120     {
 121       sign_ = 2;
 122       num_ = 1;
 123     }
 124   else if (!num_)
 125     {
 126       sign_ = 0;
 127       den_ = 1;
 128     }
 129   else
 130     {
 131       int g = gcd (num_, den_);
 132
 133       num_ /= g;
 134       den_ /= g;
 135     }
 136 }
 137 int
 138 Rational::sign () const
 139 {
 140   return ::sign (sign_);
 141 }
 142
 143 int
 144 Rational::compare (Rational const &r, Rational const &s)
 145 {
 146   if (r.sign_ < s.sign_)
 147     return -1;
 148   else if (r.sign_ > s.sign_)
 149     return 1;
 150   else if (r.is_infinity ())
 151     return 0;
 152   else if (r.sign_ == 0)
 153     return 0;
 154   return r.sign_ * ::sign (int (r.num_ * s.den_) - int (s.num_ * r.den_));
 155 }
 156
 157 int
 158 compare (Rational const &r, Rational const &s)
 159 {
 160   return Rational::compare (r, s);
 161 }
 162
 163 Rational &
 164 Rational::operator %= (Rational r)
 165 {
 166   *this = mod_rat (r);
 167   return *this;
 168 }
 169
 170 Rational &
 171 Rational::operator += (Rational r)
 172 {
 173   if (is_infinity ())
 174     ;
 175   else if (r.is_infinity ())
 176     *this = r;
 177   else
 178     {
 179       int lcm = (den_ / gcd (r.den_, den_)) * r.den_;
 180       int n = sign_ * num_ * (lcm / den_) + r.sign_ * r.num_ * (lcm / r.den_);
 181       int d = lcm;
 182       sign_ = ::sign (n) * ::sign (d);
 183       num_ = ::abs (n);
 184       den_ = ::abs (d);
 185       normalize ();
 186     }
 187   return *this;
 188 }
 189
 190 /*
 191   copied from libg++ 2.8.0
 192 */
 193 Rational::Rational (double x)
 194 {
 195   if (x != 0.0)
 196     {
 197       sign_ = ::sign (x);
 198       x *= sign_;
 199
 200       int expt;
 201       double mantissa = frexp (x, &expt);
 202
 203       const int FACT = 1 << 20;
 204
 205       /*
 206         Thanks to Afie for this too simple  idea.
 207
 208         do not blindly substitute by libg++ code, since that uses
 209         arbitrary-size integers.  The rationals would overflow too
 210         easily.
 211       */
 212
 213       num_ = (unsigned int) (mantissa * FACT);
 214       den_ = (unsigned int) FACT;
 215       normalize ();
 216       if (expt < 0)
 217         den_ <<= -expt;
 218       else
 219         num_ <<= expt;
 220       normalize ();
 221     }
 222   else
 223     {
 224       num_ = 0;
 225       den_ = 1;
 226       sign_ = 0;
 227       normalize ();
 228     }
 229 }
 230
 231 void
 232 Rational::invert ()
 233 {
 234   int r (num_);
 235   num_ = den_;
 236   den_ = r;
 237 }
 238
 239 Rational &
 240 Rational::operator *= (Rational r)
 241 {
 242   sign_ *= ::sign (r.sign_);
 243   if (r.is_infinity ())
 244     {
 245       sign_ = sign () * 2;
 246       goto exit_func;
 247     }
 248
 249   num_ *= r.num_;
 250   den_ *= r.den_;
 251
 252   normalize ();
 253  exit_func:
 254   return *this;
 255 }
 256
 257 Rational &
 258 Rational::operator /= (Rational r)
 259 {
 260   r.invert ();
 261   return (*this *= r);
 262 }
 263
 264 void
 265 Rational::negate ()
 266 {
 267   sign_ *= -1;
 268 }
 269
 270 Rational &
 271 Rational::operator -= (Rational r)
 272 {
 273   r.negate ();
 274   return (*this += r);
 275 }
 276
 277 string
 278 Rational::to_string () const
 279 {
 280   if (is_infinity ())
 281     {
 282       string s (sign_ > 0 ? "" : "-");
 283       return string (s + "infinity");
 284     }
 285
 286   string s = ::to_string (num ());
 287   if (den () != 1 && num ())
 288     s += "/" + ::to_string (den ());
 289   return s;
 290 }
 291
 292 int
 293 Rational::to_int () const
 294 {
 295   return (int) num () / den ();
 296 }
 297
 298 int
 299 sign (Rational r)
 300 {
 301   return r.sign ();
 302 }
 303
 304 bool
 305 Rational::is_infinity () const
 306 {
 307   return sign_ == 2 || sign_ == -2;
 308 }