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