rational.h

Go to the documentation of this file.

/* Copyright (C) 2012 
Richard Larsson <ric.larsson@gmail.com>
 
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
 
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA. */
 
#ifndef rational_h
#define rational_h
 
#include "array.h"
#include "bifstream.h"
#include "bofstream.h"
#include "math_funcs.h"
#include "matpack.h"
#include <ostream>
 
 
constexpr Index gcd(Index a, Index b) noexcept {
  if (b == 0)
    return a;
  return gcd(b, a%b);
}
 
 
class Rational {
 public:
  constexpr Rational(const Index nom = 0, const Index denom = 1) noexcept
  : mnom(denom ? nom : 0), mdenom(denom) {
    const auto div = gcd(nom, denom);
    if (div) {
      mnom /= div;
      mdenom /= div;
    }
  }
  
  explicit Rational(const String& s);
 
  [[nodiscard]] constexpr Index Nom() const noexcept { return mnom; }
 
  [[nodiscard]] constexpr Index Denom() const noexcept { return mdenom; }
  
  constexpr Index& Nom() noexcept { return mnom; }
  
  constexpr Index& Denom() noexcept { return mdenom; }
  
  constexpr void Nom(Index x) noexcept { mnom = x; }
  
  constexpr void Denom(Index x) noexcept { mdenom = x; }
  
  void simplify_in_place() noexcept;
 
  [[nodiscard]] constexpr bool isUndefined() const noexcept { return (mdenom == 0); }
 
  [[nodiscard]] constexpr bool isDefined() const noexcept { return not isUndefined(); }
 
  [[nodiscard]] constexpr bool isIndex(int n=1) const noexcept {
    return isDefined() and not bool((n*mnom) % mdenom);
  }
 
  [[nodiscard]] constexpr Index toIndex(int n=1) const noexcept {
    return (n*mnom) / mdenom;;
  }
 
  [[nodiscard]] constexpr Numeric toNumeric() const noexcept {
    return Numeric(mnom) / Numeric(mdenom);
  }
  
  [[nodiscard]] constexpr int toInt(int n=1) const noexcept { return int(toIndex(n)); }
  
  constexpr Rational& operator+=(const Rational& a) noexcept {
    mnom = mnom * a.Denom() + a.Nom() * mdenom;
    mdenom *= a.Denom();
    return *this;
  }
 
  constexpr Rational& operator+=(const Index& a) noexcept {
    mnom += mdenom * a;
    return *this;
  }
  
  constexpr Rational& operator+=(const int& a) noexcept {
    mnom += mdenom * a;
    return *this;
  }
 
  constexpr Rational& operator-=(const Rational& a) noexcept {
    mnom = mnom * a.Denom() - a.Nom() * mdenom;
    mdenom *= a.Denom();
    return *this;
  }
 
  constexpr Rational& operator-=(const Index& a) noexcept {
    mnom -= mdenom * a;
    return *this;
  }
  
  constexpr Rational& operator-=(const int& a) noexcept {
    mnom -= mdenom * a;
    return *this;
  }
 
  constexpr Rational& operator/=(const Rational& a) noexcept {
    mnom *= a.Denom();
    mdenom *= a.Nom();
    return *this;
  }
 
  constexpr Rational& operator/=(const Index& a) noexcept {
    mdenom *= a;
    return *this;
  }
  
  constexpr Rational& operator/=(const int& a) noexcept {
    mdenom *= a;
    return *this;
  }
 
  constexpr Rational& operator*=(const Rational& a) noexcept {
    mnom *= a.Nom();
    mdenom *= a.Denom();
    return *this;
  }
 
  constexpr Rational& operator*=(const Index& a) noexcept {
    mnom *= a;
    return *this;
  }
  
  constexpr Rational& operator*=(const int& a) noexcept {
    mnom *= a;
    return *this;
  }
 
  constexpr Rational operator++(int) const noexcept{
    return Rational(mnom + mdenom, mdenom);
  }
 
  constexpr Rational operator--(int) const noexcept {
    return Rational(mnom - mdenom, mdenom);
  }
 
  constexpr Rational& operator++() noexcept {
    mnom += mdenom;
    return *this;
  }
 
  constexpr Rational& operator--() noexcept {
    mnom -= mdenom;
    return *this;
  }
 
  explicit constexpr operator bool() const noexcept {
    return isDefined() and bool(mnom);
  }
 
  explicit constexpr operator Numeric() const noexcept { return toNumeric(); }
 
  explicit constexpr operator Index() const noexcept { return toIndex(); }
 
  explicit constexpr operator int() const noexcept { return toInt(); }
  
  bifstream& read(bifstream& bif) {
    bif >> mnom >> mdenom;
    return bif;
  }
  
  bofstream& write(bofstream& bof) const {
    bof << mnom << mdenom;
    return bof;
  }
  
  constexpr Rational& fixSign() noexcept {
    if (mdenom < 0) {
      mnom = -mnom;
      mdenom = -mdenom;
    }
    return *this;
  }
 
 private:
  // Rational is supposed to be used rationally ( mnom / mdenom )
  Index mnom;
  Index mdenom;
};  // Rational;
 
constexpr Rational reduce_by_gcd(Rational a) noexcept {
  const Index div = gcd(a.Nom(), a.Denom());
  if (div)
    return Rational(a.Nom() / div, a.Denom() / div);
  return a;
}
 
constexpr Rational numeric2rational(Numeric x, size_t maxdec=4) noexcept {
  Index nom=0, denom=1;
  
  // Keep track of sign independently
  const bool signchange = x < 0;
  x = signchange ? -x : x;
  
  // Add numbers by keeping the floor
  size_t i=0;
  do {
    const auto xi=Index(x);
    nom += xi;
    x = 10 * (x - Numeric(xi));
    nom *= 10;
    denom *= 10;
    i++;
  } while (i<=maxdec);
  
  // Fix possible rounding error
  if (x >= 5)
    nom += 10;
  
  // Change sign or not
  if (signchange)
    return Rational(-nom, denom);
  return Rational(nom, denom);
}
 
 
// An undefined rational to be used everywhere for a undefined rationals
#define RATIONAL_UNDEFINED Rational(0, 0)
 
constexpr Rational operator-(const Rational a) noexcept {
  return Rational(-a.Nom(), a.Denom());
}
 
constexpr Rational operator+(const Rational a) noexcept { return a; }
 
constexpr Rational operator+(const Rational a, const Rational b) noexcept {
  return (a.Denom() == b.Denom())
             ? Rational(a.Nom() + b.Nom(), a.Denom())
             : Rational(a.Nom() * b.Denom() + b.Nom() * a.Denom(),
                        a.Denom() * b.Denom());
}
 
constexpr Rational operator+(const Rational a, Index b) noexcept {
  return Rational(a.Nom() + b * a.Denom(), a.Denom());
}
 
constexpr Rational operator+(const Rational a, int b) noexcept {
  return Rational(a.Nom() + b * a.Denom(), a.Denom());
}
 
constexpr Rational operator+(Index b, const Rational a) noexcept { return operator+(a, b); }
 
constexpr Rational operator+(int b, const Rational a) noexcept { return operator+(a, b); }
 
constexpr Rational operator-(const Rational a, const Rational b) noexcept {
  return (a.Denom() == b.Denom())
             ? Rational(a.Nom() - b.Nom(), a.Denom())
             : Rational(a.Nom() * b.Denom() - b.Nom() * a.Denom(),
                        a.Denom() * b.Denom());
}
 
constexpr Rational operator-(const Rational a, Index b) noexcept {
  return Rational(a.Nom() - b * a.Denom(), a.Denom());
}
 
constexpr Rational operator-(const Rational a, int b) noexcept {
  return Rational(a.Nom() - b * a.Denom(), a.Denom());
}
 
constexpr Rational operator-(Index b, const Rational a) noexcept {
  return Rational(-a.Nom() + b * a.Denom(), a.Denom());
}
 
constexpr Rational operator-(int b, const Rational a) noexcept {
  return Rational(-a.Nom() + b * a.Denom(), a.Denom());
}
 
constexpr Rational operator/(const Rational a, const Rational b) noexcept {
  return Rational(a.Nom() * b.Denom(), a.Denom() * b.Nom());
}
 
constexpr Rational operator/(const Rational a, Index b) noexcept {
  return Rational(a.Nom(), a.Denom() * b);
}
 
constexpr Rational operator/(const Rational a, int b) noexcept {
  return Rational(a.Nom(), a.Denom() * b);
}
 
constexpr Rational operator/(Index b, const Rational a) noexcept {
  return Rational(a.Denom() * b, a.Nom());
}
 
constexpr Rational operator/(int b, const Rational a) noexcept {
  return Rational(a.Denom() * b, a.Nom());
}
 
constexpr Rational operator*(const Rational a, const Rational b) noexcept {
  return Rational(a.Nom() * b.Nom(), a.Denom() * b.Denom());
}
 
constexpr Rational operator*(const Rational a, Index b) noexcept {
  return Rational(a.Nom() * b, a.Denom());
}
 
constexpr Rational operator*(const Rational a, int b) noexcept {
  return Rational(a.Nom() * b, a.Denom());
}
 
constexpr Rational operator*(Index b, const Rational a) noexcept { return operator*(a, b); }
 
constexpr Rational operator*(int b, const Rational a) noexcept { return operator*(a, b); }
 
constexpr Rational operator%(const Rational a, const Rational b) noexcept {
  return (a.Denom() == b.Denom())
             ? Rational(a.Nom() % b.Nom(), a.Denom())
             : Rational((a.Nom() * b.Denom()) % (a.Denom() * b.Nom()),
                        a.Denom() * b.Denom());
}
 
constexpr Rational operator%(const Rational a, Index b) noexcept {
  return Rational(a.Nom() % (a.Denom() * b), a.Denom());
}
 
constexpr Rational operator%(const Rational a, int b) noexcept {
  return Rational(a.Nom() % (a.Denom() * b), a.Denom());
}
 
constexpr Rational operator%(Index b, const Rational a) noexcept {
  return Rational((b * a.Denom()) % a.Nom(), a.Denom());
}
 
constexpr Rational operator%(int b, const Rational a) noexcept {
  return Rational((b * a.Denom()) % a.Nom(), a.Denom());
}
 
constexpr bool operator==(const Rational a, const Rational b) noexcept {
  return a.isDefined() and b.isDefined() and
         a.Nom() * b.Denom() == a.Denom() * b.Nom();
}
 
constexpr bool operator!=(const Rational a, const Rational b) noexcept {
  return not (a.isDefined() and b.isDefined() and operator==(a, b));
}
 
constexpr bool operator<(const Rational a, const Rational b) noexcept {
  return a.isDefined() and b.isDefined() and
         a.Nom() * b.Denom() < a.Denom() * b.Nom();
}
 
constexpr bool operator>(const Rational a, const Rational b) noexcept { return operator<(b, a) and a.isDefined() and b.isDefined(); }
 
constexpr bool operator<=(const Rational a, const Rational b) noexcept {
  return not operator>(a, b) and a.isDefined() and b.isDefined();
}
 
constexpr bool operator>=(const Rational a, const Rational b) noexcept {
  return not operator<(a, b) and a.isDefined() and b.isDefined();
}
 
constexpr bool operator!(const Rational a) noexcept { return a.Nom() and a.isDefined(); }
 
inline Numeric sqrt(const Rational r) { return std::sqrt(r.toNumeric()); }
 
inline Numeric pow(const Rational base, Numeric exp) {
  return std::pow(base.toNumeric(), exp);
}
 
inline Numeric pow(Numeric base, const Rational exp) {
  return std::pow(base, exp.toNumeric());
}
 
inline Numeric pow(const Rational base, const Rational exp) {
  return pow(base, exp.toNumeric());
}
 
std::ostream& operator<<(std::ostream& os, const Rational& a);
 
std::istream& operator>>(std::istream& is, Rational& a);
 
constexpr bool operator<(const Index a, const Rational b) noexcept {
  return Rational(a, 1) < b and b.isDefined();
}
 
constexpr bool operator<(const int a, const Rational b) noexcept {
  return Rational(a, 1) < b and b.isDefined();
}
 
constexpr bool operator<(const Rational a, const Index b) noexcept {
  return a < Rational(b, 1) and a.isDefined();
}
 
constexpr bool operator<(const Rational a, const int b) noexcept {
  return a < Rational(b, 1) and a.isDefined();
}
 
constexpr bool operator>(const Index a, const Rational b) noexcept {
  return Rational(a, 1) > b and b.isDefined();
}
 
constexpr bool operator>(const int a, const Rational b) noexcept {
  return Rational(a, 1) > b and b.isDefined();
}
 
constexpr bool operator>(const Rational a, const Index b) noexcept {
  return a > Rational(b, 1) and a.isDefined();
}
 
constexpr bool operator>(const Rational a, const int b) noexcept {
  return a > Rational(b, 1) and a.isDefined();
}
 
constexpr bool operator==(const Rational a, const Index b) noexcept {
  return a == Rational(b, 1) and a.isDefined();
}
 
constexpr bool operator==(const Rational a, const int b) noexcept {
  return a == Rational(b, 1) and a.isDefined();
}
 
constexpr bool operator!=(const Rational a, const Index b) noexcept {
  return not(a == b) and a.isDefined();
}
 
constexpr bool operator!=(const Rational a, const int b) noexcept {
  return not(a == b) and a.isDefined();
}
 
constexpr Rational abs(const Rational a) noexcept {
  return a < 0 ? -a : a;
}
 
constexpr Rational max(const Rational a, const Rational b) noexcept {
  return a < b ? b : a;
}  // Let other operators find out if this is allowed instead
 
constexpr Rational min(const Rational a, const Rational b) noexcept {
  return a < b ? a : b;
}  // Let other operators find out if this is allowed instead
 
using ArrayOfRational = Array<Rational>;
 
constexpr Rational operator ""_2(unsigned long long int n) noexcept {
  return Rational(n, 2);
};
 
constexpr bool iseven(const Rational r) noexcept {
  return 0 == (r % 2);
}
 
constexpr Numeric operator*(Rational y, Numeric x) noexcept {return y.toNumeric() * x;}
constexpr Numeric operator*(Numeric x, Rational y) noexcept {return x * y.toNumeric();}
 
constexpr Numeric operator/(Rational y, Numeric x) noexcept {return y.toNumeric() / x;}
constexpr Numeric operator/(Numeric x, Rational y) noexcept {return x / y.toNumeric();}
 
constexpr Numeric operator+(Rational y, Numeric x) noexcept {return y.toNumeric() + x;}
constexpr Numeric operator+(Numeric x, Rational y) noexcept {return x + y.toNumeric();}
 
constexpr Numeric operator-(Rational y, Numeric x) noexcept {return y.toNumeric() - x;}
constexpr Numeric operator-(Numeric x, Rational y) noexcept {return x - y.toNumeric();}
 
#endif  // rational_h