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 (int n, int d)
  65 {
  66   sign_ = ::sign (n) * ::sign (d);
  67   num_ = ::abs (n);
  68   den_ = ::abs (d);
  69   normalize ();
  70 }
  71
  72 Rational::Rational (int n)
  73 {
  74   sign_ = ::sign (n);
  75   num_ = ::abs (n);
  76   den_ = 1;
  77 }
  78
  79
  80 void
  81 Rational::set_infinite (int s)
  82 {
  83   sign_ = ::sign (s) * 2;
  84   num_ = 1;
  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
 112 /*
 113   copy & paste from scm_gcd (GUILE).
 114  */
 115 static int
 116 gcd (long u, long v)
 117 {
 118   long result = 0;
 119   if (u == 0)
 120     result = v;
 121   else if (v == 0)
 122     result = u;
 123   else
 124     {
 125       long k = 1;
 126       long t;
 127       /* Determine a common factor 2^k */
 128       while (!(1 & (u | v)))
 129         {
 130           k <<= 1;
 131           u >>= 1;
 132           v >>= 1;
 133         }
 134       /* Now, any factor 2^n can be eliminated */
 135       if (u & 1)
 136         t = -v;
 137       else
 138         {
 139           t = u;
 140         b3:
 141           t = t >> 1;
 142         }
 143       if (!(1 & t))
 144         goto b3;
 145       if (t > 0)
 146         u = t;
 147       else
 148         v = -t;
 149       t = u - v;
 150       if (t != 0)
 151         goto b3;
 152       result = u * k;
 153     }
 154
 155   return result;
 156 }
 157
 158
 159 void
 160 Rational::normalize ()
 161 {
 162   if (!sign_)
 163     {
 164       den_ = 1;
 165       num_ = 0;
 166     }
 167   else if (!den_)
 168     {
 169       sign_ = 2;
 170       num_ = 1;
 171     }
 172   else if (!num_)
 173     {
 174       sign_ = 0;
 175       den_ = 1;
 176     }
 177   else
 178     {
 179       int g = gcd (num_, den_);
 180
 181       num_ /= g;
 182       den_ /= g;
 183     }
 184 }
 185 int
 186 Rational::sign () const
 187 {
 188   return ::sign (sign_);
 189 }
 190
 191 int
 192 Rational::compare (Rational const &r, Rational const &s)
 193 {
 194   if (r.sign_ < s.sign_)
 195     return -1;
 196   else if (r.sign_ > s.sign_)
 197     return 1;
 198   else if (r.is_infinity ())
 199     return 0;
 200   else if (r.sign_ == 0)
 201     return 0;
 202   return r.sign_ * ::sign (int (r.num_ * s.den_) - int (s.num_ * r.den_));
 203 }
 204
 205 int
 206 compare (Rational const &r, Rational const &s)
 207 {
 208   return Rational::compare (r, s);
 209 }
 210
 211 Rational &
 212 Rational::operator %= (Rational r)
 213 {
 214   *this = mod_rat (r);
 215   return *this;
 216 }
 217
 218 Rational &
 219 Rational::operator += (Rational r)
 220 {
 221   if (is_infinity ())
 222     ;
 223   else if (r.is_infinity ())
 224     *this = r;
 225   else
 226     {
 227       int lcm = (den_ / gcd (r.den_, den_)) * r.den_;
 228       int n = sign_ * num_ * (lcm / den_) + r.sign_ * r.num_ * (lcm / r.den_);
 229       int d = lcm;
 230       sign_ = ::sign (n) * ::sign (d);
 231       num_ = ::abs (n);
 232       den_ = ::abs (d);
 233       normalize ();
 234     }
 235   return *this;
 236 }
 237
 238 /*
 239   copied from libg++ 2.8.0
 240 */
 241 Rational::Rational (double x)
 242 {
 243   if (x != 0.0)
 244     {
 245       sign_ = ::sign (x);
 246       x *= sign_;
 247
 248       int expt;
 249       double mantissa = frexp (x, &expt);
 250
 251       const int FACT = 1 << 20;
 252
 253       /*
 254         Thanks to Afie for this too simple  idea.
 255
 256         do not blindly substitute by libg++ code, since that uses
 257         arbitrary-size integers.  The rationals would overflow too
 258         easily.
 259       */
 260
 261       num_ = (unsigned int) (mantissa * FACT);
 262       den_ = (unsigned int) FACT;
 263       normalize ();
 264       if (expt < 0)
 265         den_ <<= -expt;
 266       else
 267         num_ <<= expt;
 268       normalize ();
 269     }
 270   else
 271     {
 272       num_ = 0;
 273       den_ = 1;
 274       sign_ = 0;
 275       normalize ();
 276     }
 277 }
 278
 279 void
 280 Rational::invert ()
 281 {
 282   int r (num_);
 283   num_ = den_;
 284   den_ = r;
 285 }
 286
 287 Rational &
 288 Rational::operator *= (Rational r)
 289 {
 290   sign_ *= ::sign (r.sign_);
 291   if (r.is_infinity ())
 292     {
 293       sign_ = sign () * 2;
 294       goto exit_func;
 295     }
 296
 297   num_ *= r.num_;
 298   den_ *= r.den_;
 299
 300   normalize ();
 301  exit_func:
 302   return *this;
 303 }
 304
 305 Rational &
 306 Rational::operator /= (Rational r)
 307 {
 308   r.invert ();
 309   return (*this *= r);
 310 }
 311
 312 void
 313 Rational::negate ()
 314 {
 315   sign_ *= -1;
 316 }
 317
 318 Rational &
 319 Rational::operator -= (Rational r)
 320 {
 321   r.negate ();
 322   return (*this += r);
 323 }
 324
 325 string
 326 Rational::to_string () const
 327 {
 328   if (is_infinity ())
 329     {
 330       string s (sign_ > 0 ? "" : "-");
 331       return string (s + "infinity");
 332     }
 333
 334   string s = ::to_string (num ());
 335   if (den () != 1 && num ())
 336     s += "/" + ::to_string (den ());
 337   return s;
 338 }
 339
 340 int
 341 Rational::to_int () const
 342 {
 343   return (int) num () / den ();
 344 }
 345
 346 int
 347 sign (Rational r)
 348 {
 349   return r.sign ();
 350 }
 351
 352 bool
 353 Rational::is_infinity () const
 354 {
 355   return sign_ == 2 || sign_ == -2;
 356 }