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