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