flower/rational.cc

   1 /*
   2   rational.cc -- implement Rational
   3
   4   source file of the Flower Library
   5
   6   (c) 1997--2007 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 double
  20 Rational::to_double () const
  21 {
  22   if (sign_ == -1 || sign_ == 1 || sign_ == 0)
  23     return ((double)sign_) * num_ / den_;
  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
  35 #ifdef STREAM_SUPPORT
  36 ostream &
  37 operator << (ostream &o, Rational r)
  38 {
  39   o << r.string ();
  40   return o;
  41 }
  42 #endif
  43
  44 Rational
  45 Rational::abs () const
  46 {
  47   return Rational (num_, den_);
  48 }
  49
  50 Rational
  51 Rational::trunc_rat () const
  52 {
  53   if (is_infinity())
  54     return *this;
  55   return Rational ((num_ - (num_ % den_)) * sign_, den_);
  56 }
  57
  58 Rational::Rational ()
  59 {
  60   sign_ = 0;
  61   num_ = den_ = 1;
  62 }
  63
  64 Rational::Rational (I64 n, I64 d)
  65 {
  66   sign_ = ::sign (n) * ::sign (d);
  67   num_ = ::abs (n);
  68   den_ = ::abs (d);
  69   normalize ();
  70 }
  71
  72 Rational::Rational (I64 n)
  73 {
  74   sign_ = ::sign (n);
  75   num_ = ::abs (n);
  76   den_ = 1;
  77 }
  78
  79 Rational::Rational (U64 n)
  80 {
  81   sign_ = 1;
  82   num_ = n;
  83   den_ = 1;
  84 }
  85
  86 Rational::Rational (int n)
  87 {
  88   sign_ = ::sign (n);
  89   num_ = ::abs (n);
  90   den_ = 1;
  91 }
  92
  93
  94 void
  95 Rational::set_infinite (int s)
  96 {
  97   sign_ = ::sign (s) * 2;
  98   num_ = 1;
  99 }
 100
 101 Rational
 102 Rational::operator - () const
 103 {
 104   Rational r (*this);
 105   r.negate ();
 106   return r;
 107 }
 108
 109 Rational
 110 Rational::div_rat (Rational div) const
 111 {
 112   Rational r (*this);
 113   r /= div;
 114   return r.trunc_rat ();
 115 }
 116
 117 Rational
 118 Rational::mod_rat (Rational div) const
 119 {
 120   Rational r (*this);
 121   r = (r / div - r.div_rat (div)) * div;
 122   return r;
 123 }
 124
 125
 126 /*
 127   copy & paste from scm_gcd (GUILE).
 128  */
 129 static I64
 130 gcd (I64 u, I64 v)
 131 {
 132   I64 result = 0;
 133   if (u == 0)
 134     result = v;
 135   else if (v == 0)
 136     result = u;
 137   else
 138     {
 139       I64 k = 1;
 140       I64 t;
 141       /* Determine a common factor 2^k */
 142       while (!(1 & (u | v)))
 143         {
 144           k <<= 1;
 145           u >>= 1;
 146           v >>= 1;
 147         }
 148       /* Now, any factor 2^n can be eliminated */
 149       if (u & 1)
 150         t = -v;
 151       else
 152         {
 153           t = u;
 154         b3:
 155           t = t >> 1;
 156         }
 157       if (!(1 & t))
 158         goto b3;
 159       if (t > 0)
 160         u = t;
 161       else
 162         v = -t;
 163       t = u - v;
 164       if (t != 0)
 165         goto b3;
 166       result = u * k;
 167     }
 168
 169   return result;
 170 }
 171
 172
 173 void
 174 Rational::normalize ()
 175 {
 176   if (!sign_)
 177     {
 178       den_ = 1;
 179       num_ = 0;
 180     }
 181   else if (!den_)
 182     {
 183       sign_ = 2;
 184       num_ = 1;
 185     }
 186   else if (!num_)
 187     {
 188       sign_ = 0;
 189       den_ = 1;
 190     }
 191   else
 192     {
 193       I64 g = gcd (num_, den_);
 194
 195       num_ /= g;
 196       den_ /= g;
 197     }
 198 }
 199 int
 200 Rational::sign () const
 201 {
 202   return ::sign (sign_);
 203 }
 204
 205 int
 206 Rational::compare (Rational const &r, Rational const &s)
 207 {
 208   if (r.sign_ < s.sign_)
 209     return -1;
 210   else if (r.sign_ > s.sign_)
 211     return 1;
 212   else if (r.is_infinity ())
 213     return 0;
 214   else if (r.sign_ == 0)
 215     return 0;
 216   return r.sign_ * ::sign ((I64) (r.num_ * s.den_) - (I64) (s.num_ * r.den_));
 217 }
 218
 219 int
 220 compare (Rational const &r, Rational const &s)
 221 {
 222   return Rational::compare (r, s);
 223 }
 224
 225 Rational &
 226 Rational::operator %= (Rational r)
 227 {
 228   *this = mod_rat (r);
 229   return *this;
 230 }
 231
 232 Rational &
 233 Rational::operator += (Rational r)
 234 {
 235   if (is_infinity ())
 236     ;
 237   else if (r.is_infinity ())
 238     *this = r;
 239   else
 240     {
 241       I64 lcm = (den_ / gcd (r.den_, den_)) * r.den_;
 242       I64 n = sign_ * num_ * (lcm / den_) + r.sign_ * r.num_ * (lcm / r.den_);
 243       I64 d = lcm;
 244       sign_ = ::sign (n) * ::sign (d);
 245       num_ = ::abs (n);
 246       den_ = ::abs (d);
 247       normalize ();
 248     }
 249   return *this;
 250 }
 251
 252 /*
 253   copied from libg++ 2.8.0
 254 */
 255 Rational::Rational (double x)
 256 {
 257   if (x != 0.0)
 258     {
 259       sign_ = ::sign (x);
 260       x *= sign_;
 261
 262       int expt;
 263       double mantissa = frexp (x, &expt);
 264
 265       const int FACT = 1 << 20;
 266
 267       /*
 268         Thanks to Afie for this too simple  idea.
 269
 270         do not blindly substitute by libg++ code, since that uses
 271         arbitrary-size integers.  The rationals would overflow too
 272         easily.
 273       */
 274
 275       num_ = (U64) (mantissa * FACT);
 276       den_ = (U64) FACT;
 277       normalize ();
 278       if (expt < 0)
 279         den_ <<= -expt;
 280       else
 281         num_ <<= expt;
 282       normalize ();
 283     }
 284   else
 285     {
 286       num_ = 0;
 287       den_ = 1;
 288       sign_ = 0;
 289       normalize ();
 290     }
 291 }
 292
 293 void
 294 Rational::invert ()
 295 {
 296   I64 r (num_);
 297   num_ = den_;
 298   den_ = r;
 299 }
 300
 301 Rational &
 302 Rational::operator *= (Rational r)
 303 {
 304   sign_ *= ::sign (r.sign_);
 305   if (r.is_infinity ())
 306     {
 307       sign_ = sign () * 2;
 308       goto exit_func;
 309     }
 310
 311   num_ *= r.num_;
 312   den_ *= r.den_;
 313
 314   normalize ();
 315  exit_func:
 316   return *this;
 317 }
 318
 319 Rational &
 320 Rational::operator /= (Rational r)
 321 {
 322   r.invert ();
 323   return (*this *= r);
 324 }
 325
 326 void
 327 Rational::negate ()
 328 {
 329   sign_ *= -1;
 330 }
 331
 332 Rational &
 333 Rational::operator -= (Rational r)
 334 {
 335   r.negate ();
 336   return (*this += r);
 337 }
 338
 339 string
 340 Rational::to_string () const
 341 {
 342   if (is_infinity ())
 343     {
 344       string s (sign_ > 0 ? "" : "-");
 345       return string (s + "infinity");
 346     }
 347
 348   string s = ::to_string (num ());
 349   if (den () != 1 && num ())
 350     s += "/" + ::to_string (den ());
 351   return s;
 352 }
 353
 354 int
 355 Rational::to_int () const
 356 {
 357   return (int) num () / den ();
 358 }
 359
 360 int
 361 sign (Rational r)
 362 {
 363   return r.sign ();
 364 }
 365
 366 bool
 367 Rational::is_infinity () const
 368 {
 369   return sign_ == 2 || sign_ == -2;
 370 }