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