flower/rational.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 1997--2010 Han-Wen Nienhuys <hanwen@xs4all.nl>
   5
   6   LilyPond is free software: you can redistribute it and/or modify
   7   it under the terms of the GNU General Public License as published by
   8   the Free Software Foundation, either version 3 of the License, or
   9   (at your option) any later version.
  10
  11   LilyPond is distributed in the hope that it will be useful,
  12   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14   GNU General Public License for more details.
  15
  16   You should have received a copy of the GNU General Public License
  17   along with LilyPond.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19
  20 #include "rational.hh"
  21
  22 #include <cmath>
  23 #include <cassert>
  24 #include <cstdlib>
  25 using namespace std;
  26
  27 #include "string-convert.hh"
  28 #include "libc-extension.hh"
  29
  30 double
  31 Rational::to_double () const
  32 {
  33   if (sign_ == -1 || sign_ == 1 || sign_ == 0)
  34     return ((double)sign_) * num_ / den_;
  35   if (sign_ == -2)
  36     return -HUGE_VAL;
  37   else if (sign_ == 2)
  38     return HUGE_VAL;
  39   else
  40     assert (false);
  41
  42   return 0.0;
  43 }
  44
  45
  46 #ifdef STREAM_SUPPORT
  47 ostream &
  48 operator << (ostream &o, Rational r)
  49 {
  50   o << r.string ();
  51   return o;
  52 }
  53 #endif
  54
  55 Rational
  56 Rational::abs () const
  57 {
  58   return Rational (num_, den_);
  59 }
  60
  61 Rational
  62 Rational::trunc_rat () const
  63 {
  64   if (is_infinity())
  65     return *this;
  66   return Rational ((num_ - (num_ % den_)) * sign_, den_);
  67 }
  68
  69 Rational::Rational ()
  70 {
  71   sign_ = 0;
  72   num_ = den_ = 1;
  73 }
  74
  75 Rational::Rational (I64 n, I64 d)
  76 {
  77   sign_ = ::sign (n) * ::sign (d);
  78   num_ = ::abs (n);
  79   den_ = ::abs (d);
  80   normalize ();
  81 }
  82
  83 Rational::Rational (I64 n)
  84 {
  85   sign_ = ::sign (n);
  86   num_ = ::abs (n);
  87   den_ = 1;
  88 }
  89
  90 Rational::Rational (U64 n)
  91 {
  92   sign_ = 1;
  93   num_ = n;
  94   den_ = 1;
  95 }
  96
  97 Rational::Rational (int n)
  98 {
  99   sign_ = ::sign (n);
 100   num_ = ::abs (n);
 101   den_ = 1;
 102 }
 103
 104
 105 void
 106 Rational::set_infinite (int s)
 107 {
 108   sign_ = ::sign (s) * 2;
 109   num_ = 1;
 110 }
 111
 112 Rational
 113 Rational::operator - () const
 114 {
 115   Rational r (*this);
 116   r.negate ();
 117   return r;
 118 }
 119
 120 Rational
 121 Rational::div_rat (Rational div) const
 122 {
 123   Rational r (*this);
 124   r /= div;
 125   return r.trunc_rat ();
 126 }
 127
 128 Rational
 129 Rational::mod_rat (Rational div) const
 130 {
 131   Rational r (*this);
 132   r = (r / div - r.div_rat (div)) * div;
 133   return r;
 134 }
 135
 136
 137 /*
 138   copy & paste from scm_gcd (GUILE).
 139  */
 140 static I64
 141 gcd (I64 u, I64 v)
 142 {
 143   I64 result = 0;
 144   if (u == 0)
 145     result = v;
 146   else if (v == 0)
 147     result = u;
 148   else
 149     {
 150       I64 k = 1;
 151       I64 t;
 152       /* Determine a common factor 2^k */
 153       while (!(1 & (u | v)))
 154         {
 155           k <<= 1;
 156           u >>= 1;
 157           v >>= 1;
 158         }
 159       /* Now, any factor 2^n can be eliminated */
 160       if (u & 1)
 161         t = -v;
 162       else
 163         {
 164           t = u;
 165         b3:
 166           t = t >> 1;
 167         }
 168       if (!(1 & t))
 169         goto b3;
 170       if (t > 0)
 171         u = t;
 172       else
 173         v = -t;
 174       t = u - v;
 175       if (t != 0)
 176         goto b3;
 177       result = u * k;
 178     }
 179
 180   return result;
 181 }
 182
 183
 184 void
 185 Rational::normalize ()
 186 {
 187   if (!sign_)
 188     {
 189       den_ = 1;
 190       num_ = 0;
 191     }
 192   else if (!den_)
 193     {
 194       sign_ = 2;
 195       num_ = 1;
 196     }
 197   else if (!num_)
 198     {
 199       sign_ = 0;
 200       den_ = 1;
 201     }
 202   else
 203     {
 204       I64 g = gcd (num_, den_);
 205
 206       num_ /= g;
 207       den_ /= g;
 208     }
 209 }
 210 int
 211 Rational::sign () const
 212 {
 213   return ::sign (sign_);
 214 }
 215
 216 int
 217 Rational::compare (Rational const &r, Rational const &s)
 218 {
 219   if (r.sign_ < s.sign_)
 220     return -1;
 221   else if (r.sign_ > s.sign_)
 222     return 1;
 223   else if (r.is_infinity ())
 224     return 0;
 225   else if (r.sign_ == 0)
 226     return 0;
 227   return r.sign_ * ::sign ((I64) (r.num_ * s.den_) - (I64) (s.num_ * r.den_));
 228 }
 229
 230 int
 231 compare (Rational const &r, Rational const &s)
 232 {
 233   return Rational::compare (r, s);
 234 }
 235
 236 Rational &
 237 Rational::operator %= (Rational r)
 238 {
 239   *this = mod_rat (r);
 240   return *this;
 241 }
 242
 243 Rational &
 244 Rational::operator += (Rational r)
 245 {
 246   if (is_infinity ())
 247     ;
 248   else if (r.is_infinity ())
 249     *this = r;
 250   else
 251     {
 252       I64 lcm = (den_ / gcd (r.den_, den_)) * r.den_;
 253       I64 n = sign_ * num_ * (lcm / den_) + r.sign_ * r.num_ * (lcm / r.den_);
 254       I64 d = lcm;
 255       sign_ = ::sign (n) * ::sign (d);
 256       num_ = ::abs (n);
 257       den_ = ::abs (d);
 258       normalize ();
 259     }
 260   return *this;
 261 }
 262
 263 /*
 264   copied from libg++ 2.8.0
 265 */
 266 Rational::Rational (double x)
 267 {
 268   if (x != 0.0)
 269     {
 270       sign_ = ::sign (x);
 271       x *= sign_;
 272
 273       int expt;
 274       double mantissa = frexp (x, &expt);
 275
 276       const int FACT = 1 << 20;
 277
 278       /*
 279         Thanks to Afie for this too simple  idea.
 280
 281         do not blindly substitute by libg++ code, since that uses
 282         arbitrary-size integers.  The rationals would overflow too
 283         easily.
 284       */
 285
 286       num_ = (U64) (mantissa * FACT);
 287       den_ = (U64) FACT;
 288       normalize ();
 289       if (expt < 0)
 290         den_ <<= -expt;
 291       else
 292         num_ <<= expt;
 293       normalize ();
 294     }
 295   else
 296     {
 297       num_ = 0;
 298       den_ = 1;
 299       sign_ = 0;
 300       normalize ();
 301     }
 302 }
 303
 304 void
 305 Rational::invert ()
 306 {
 307   I64 r (num_);
 308   num_ = den_;
 309   den_ = r;
 310 }
 311
 312 Rational &
 313 Rational::operator *= (Rational r)
 314 {
 315   sign_ *= ::sign (r.sign_);
 316   if (r.is_infinity ())
 317     {
 318       sign_ = sign () * 2;
 319       goto exit_func;
 320     }
 321
 322   num_ *= r.num_;
 323   den_ *= r.den_;
 324
 325   normalize ();
 326  exit_func:
 327   return *this;
 328 }
 329
 330 Rational &
 331 Rational::operator /= (Rational r)
 332 {
 333   r.invert ();
 334   return (*this *= r);
 335 }
 336
 337 void
 338 Rational::negate ()
 339 {
 340   sign_ *= -1;
 341 }
 342
 343 Rational &
 344 Rational::operator -= (Rational r)
 345 {
 346   r.negate ();
 347   return (*this += r);
 348 }
 349
 350 string
 351 Rational::to_string () const
 352 {
 353   if (is_infinity ())
 354     {
 355       string s (sign_ > 0 ? "" : "-");
 356       return string (s + "infinity");
 357     }
 358
 359   string s = ::to_string (num ());
 360   if (den () != 1 && num ())
 361     s += "/" + ::to_string (den ());
 362   return s;
 363 }
 364
 365 int
 366 Rational::to_int () const
 367 {
 368   return (int) num () / den ();
 369 }
 370
 371 int
 372 sign (Rational r)
 373 {
 374   return r.sign ();
 375 }
 376
 377 bool
 378 Rational::is_infinity () const
 379 {
 380   return sign_ == 2 || sign_ == -2;
 381 }