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 <numeric>
#include <ostream>
 
 
using std::gcd;
 
struct Rational {
  Index numer;
  Index denom;
 
  constexpr Rational(const Index n = 0, const Index d = 1) noexcept
  : numer(d ? n / gcd(n, d) : 0), denom(d / gcd(n, d)) {}
  
  explicit Rational(const String& s);
  
  void simplify_in_place() noexcept;
 
  [[nodiscard]] constexpr bool isUndefined() const noexcept { return (denom == 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*numer) % denom);
  }
 
  [[nodiscard]] constexpr Index toIndex(int n=1) const noexcept {
    return (n*numer) / denom;;
  }
 
  [[nodiscard]] constexpr Numeric toNumeric() const noexcept {
    return Numeric(numer) / Numeric(denom);
  }
  
  [[nodiscard]] constexpr int toInt(int n=1) const noexcept { return int(toIndex(n)); }
  
  constexpr Rational& operator+=(const Rational& a) noexcept {
    numer = numer * a.denom + a.numer * denom;
    denom *= a.denom;
    return *this;
  }
 
  constexpr Rational& operator+=(const Index& a) noexcept {
    numer += denom * a;
    return *this;
  }
  
  constexpr Rational& operator+=(const int& a) noexcept {
    numer += denom * a;
    return *this;
  }
 
  constexpr Rational& operator-=(const Rational& a) noexcept {
    numer = numer * a.denom - a.numer * denom;
    denom *= a.denom;
    return *this;
  }
 
  constexpr Rational& operator-=(const Index& a) noexcept {
    numer -= denom * a;
    return *this;
  }
  
  constexpr Rational& operator-=(const int& a) noexcept {
    numer -= denom * a;
    return *this;
  }
 
  constexpr Rational& operator/=(const Rational& a) noexcept {
    numer *= a.denom;
    denom *= a.numer;
    return *this;
  }
 
  constexpr Rational& operator/=(const Index& a) noexcept {
    denom *= a;
    return *this;
  }
  
  constexpr Rational& operator/=(const int& a) noexcept {
    denom *= a;
    return *this;
  }
 
  constexpr Rational& operator*=(const Rational& a) noexcept {
    numer *= a.numer;
    denom *= a.denom;
    return *this;
  }
 
  constexpr Rational& operator*=(const Index& a) noexcept {
    numer *= a;
    return *this;
  }
  
  constexpr Rational& operator*=(const int& a) noexcept {
    numer *= a;
    return *this;
  }
 
  constexpr Rational operator++(int) const noexcept{
    return {numer + denom, denom};
  }
 
  constexpr Rational operator--(int) const noexcept {
    return {numer - denom, denom};
  }
 
  constexpr Rational& operator++() noexcept {
    numer += denom;
    return *this;
  }
 
  constexpr Rational& operator--() noexcept {
    numer -= denom;
    return *this;
  }
 
  explicit constexpr operator bool() const noexcept {
    return isDefined() and bool(numer);
  }
 
  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 >> numer >> denom;
    return bif;
  }
  
  bofstream& write(bofstream& bof) const {
    bof << numer << denom;
    return bof;
  }
  
  constexpr Rational& fixSign() noexcept {
    if (denom < 0) {
      numer = -numer;
      denom = -denom;
    }
    return *this;
  }
 
  friend std::ostream& operator<<(std::ostream& os, const Rational& a);
  
  friend std::istream& operator>>(std::istream& is, Rational& a);
};  // Rational;
 
constexpr Rational reduce_by_gcd(Rational a) noexcept {
  const Index div = gcd(a.numer, a.denom);
  if (div)
    return {a.numer / 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
  return signchange ? Rational{-nom, denom} : 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 {-a.numer, 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.numer + b.numer, a.denom)
             : Rational(a.numer * b.denom + b.numer * a.denom,
                        a.denom * b.denom);
}
 
constexpr Rational operator+(const Rational a, Index b) noexcept {
  return {a.numer + b * a.denom, a.denom};
}
 
constexpr Rational operator+(const Rational a, int b) noexcept {
  return {a.numer + 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.numer - b.numer, a.denom)
             : Rational(a.numer * b.denom - b.numer * a.denom,
                        a.denom * b.denom);
}
 
constexpr Rational operator-(const Rational a, Index b) noexcept {
  return {a.numer - b * a.denom, a.denom};
}
 
constexpr Rational operator-(const Rational a, int b) noexcept {
  return {a.numer - b * a.denom, a.denom};
}
 
constexpr Rational operator-(Index b, const Rational a) noexcept {
  return {-a.numer + b * a.denom, a.denom};
}
 
constexpr Rational operator-(int b, const Rational a) noexcept {
  return {-a.numer + b * a.denom, a.denom};
}
 
constexpr Rational operator/(const Rational a, const Rational b) noexcept {
  return {a.numer * b.denom, a.denom * b.numer};
}
 
constexpr Rational operator/(const Rational a, Index b) noexcept {
  return {a.numer, a.denom * b};
}
 
constexpr Rational operator/(const Rational a, int b) noexcept {
  return {a.numer, a.denom * b};
}
 
constexpr Rational operator/(Index b, const Rational a) noexcept {
  return {a.denom * b, a.numer};
}
 
constexpr Rational operator/(int b, const Rational a) noexcept {
  return {a.denom * b, a.numer};
}
 
constexpr Rational operator*(const Rational a, const Rational b) noexcept {
  return {a.numer * b.numer, a.denom * b.denom};
}
 
constexpr Rational operator*(const Rational a, Index b) noexcept {
  return {a.numer * b, a.denom};
}
 
constexpr Rational operator*(const Rational a, int b) noexcept {
  return {a.numer * 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.numer % b.numer, a.denom)
             : Rational((a.numer * b.denom) % (a.denom * b.numer),
                        a.denom * b.denom);
}
 
constexpr Rational operator%(const Rational a, Index b) noexcept {
  return {a.numer % (a.denom * b), a.denom};
}
 
constexpr Rational operator%(const Rational a, int b) noexcept {
  return {a.numer % (a.denom * b), a.denom};
}
 
constexpr Rational operator%(Index b, const Rational a) noexcept {
  return {(b * a.denom) % a.numer, a.denom};
}
 
constexpr Rational operator%(int b, const Rational a) noexcept {
  return {(b * a.denom) % a.numer, a.denom};
}
 
constexpr bool operator==(const Rational a, const Rational b) noexcept {
  return a.isDefined() and b.isDefined() and
         a.numer * b.denom == a.denom * b.numer;
}
 
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.numer * b.denom < a.denom * b.numer;
}
 
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.numer 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());
}
 
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