flower/rational.cc

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