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