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