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 <cassert>
#include <ostream>
#include "array.h"
#include "bifstream.h"
#include "bofstream.h"
#include "math_funcs.h"
#include "matpack.h"


constexpr Index gcd(Index a, Index b) {
  if (b == 0)
    return a;
  else
    return gcd(b, a%b);
}


class Rational {
 public:
  explicit constexpr Rational(const Index nom = 0, const Index denom = 1)
  : mnom(denom ? nom : 0), mdenom(denom) {
    const auto div = gcd(nom, denom);
    if (div) {
      mnom /= div;
      mdenom /= div;
    }
  }
  
  explicit Rational(const String& s);

  constexpr Index Nom() const { return mnom; }

  constexpr Index Denom() const { return mdenom; }
  
  Index& Nom() { return mnom; }
  
  Index& Denom() { return mdenom; }
  
  void Nom(Index x) { mnom = x; }
  
  void Denom(Index x) { mdenom = x; }
  
  void simplify_in_place();

  constexpr bool isUndefined() const { return (mdenom == 0); }

  constexpr bool isDefined() const { return mdenom not_eq 0; }

  constexpr bool isIndex(int n=1) const {
    return isDefined() and not bool((n*mnom) % mdenom);
  }

  constexpr Index toIndex(int n=1) const {
    return isIndex(n) ? (n*mnom) / mdenom
                     : throw std::logic_error(
                           "Rational is not representative of an integer");
  }

  constexpr Numeric toNumeric() const {
    return Numeric(mnom) / Numeric(mdenom);
  }
  
  constexpr int toInt(int n=1) const { return int(toIndex(n)); }
  
  Rational& operator+=(const Rational& a) {
    mnom = mnom * a.Denom() + a.Nom() * mdenom;
    mdenom *= a.Denom();
    return *this;
  }

  Rational& operator+=(const Index& a) {
    mnom += mdenom * a;
    return *this;
  }

  Rational& operator-=(const Rational& a) {
    mnom = mnom * a.Denom() - a.Nom() * mdenom;
    mdenom *= a.Denom();
    return *this;
  }

  Rational& operator-=(const Index& a) {
    mnom -= mdenom * a;
    return *this;
  }

  Rational& operator/=(const Rational& a) {
    mnom *= a.Denom();
    mdenom *= a.Nom();
    return *this;
  }

  Rational& operator/=(const Index& a) {
    mdenom *= a;
    return *this;
  }

  Rational& operator*=(const Rational& a) {
    mnom *= a.Nom();
    mdenom *= a.Denom();
    return *this;
  }

  Rational& operator*=(const Index& a) {
    mnom *= a;
    return *this;
  }

  constexpr Rational operator++(int) const {
    return isDefined()
               ? Rational(mnom + mdenom, mdenom)
               : throw std::logic_error("Undefined Rational in operator++");
  }

  constexpr Rational operator--(int) const {
    return isDefined()
               ? Rational(mnom - mdenom, mdenom)
               : throw std::logic_error("Undefined Rational in operator--");
  }

  Rational& operator++() {
    mnom += mdenom;
    return *this;
  }

  Rational& operator--() {
    mnom -= mdenom;
    return *this;
  }

  explicit constexpr operator bool() const {
    return isDefined() and bool(mnom);
  }

  explicit constexpr operator Numeric() const { return toNumeric(); }

  explicit constexpr operator Index() const { return toIndex(); }

  explicit constexpr operator int() const { return toInt(); }
  
  bifstream& read(bifstream& bif) {
    bif >> mnom >> mdenom;
    return bif;
  }
  
  bofstream& write(bofstream& bof) const {
    bof << mnom << mdenom;
    return bof;
  }
  
  constexpr Rational& fixSign() {
    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(const Rational a) {
  const Index div = gcd(a.Nom(), a.Denom());
  if (div)
    return Rational(a.Nom() / div, a.Denom() / div);
  else
    return a;
}

constexpr Rational numeric2rational(Numeric x, size_t maxdec=4) {
  Index nom=0, denom=1;
  
  // Keep track of sign independently
  const bool signchange = x < 0;
  x = x < 0 ? -x : x;
  
  // Add numbers by keeping the floor
  size_t i=0;
  do {
    const Index 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);
  else
    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) {
  return Rational(-a.Nom(), a.Denom());
}

constexpr Rational operator+(const Rational a) { return a; }

constexpr Rational operator+(const Rational a, const Rational b) {
  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) {
  return Rational(a.Nom() + b * a.Denom(), a.Denom());
}

constexpr Rational operator+(Index b, const Rational a) { return operator+(a, b); }

constexpr Rational operator-(const Rational a, const Rational b) {
  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) {
  return Rational(a.Nom() - b * a.Denom(), a.Denom());
}

constexpr Rational operator-(Index b, const Rational a) {
  return Rational(-a.Nom() + b * a.Denom(), a.Denom());
}

constexpr Rational operator/(const Rational a, const Rational b) {
  return Rational(a.Nom() * b.Denom(), a.Denom() * b.Nom());
}

constexpr Rational operator/(const Rational a, Index b) {
  return Rational(a.Nom(), a.Denom() * b);
}

constexpr Rational operator/(Index b, const Rational a) {
  return Rational(a.Denom() * b, a.Nom());
}

constexpr Rational operator*(const Rational a, const Rational b) {
  return Rational(a.Nom() * b.Nom(), a.Denom() * b.Denom());
}

constexpr Rational operator*(const Rational a, Index b) {
  return Rational(a.Nom() * b, a.Denom());
}

constexpr Rational operator*(Index b, const Rational a) { return operator*(a, b); }

constexpr Rational operator%(const Rational a, const Rational b) {
  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) {
  return Rational(a.Nom() % (a.Denom() * b), a.Denom());
}

constexpr Rational operator%(Index b, const Rational a) {
  return Rational((b * a.Denom()) % a.Nom(), a.Denom());
}

constexpr bool operator==(const Rational a, const Rational b) {
  return a.isDefined() and b.isDefined() and
         a.Nom() * b.Denom() == a.Denom() * b.Nom();
}

constexpr bool operator!=(Rational a, Rational b) {
  return not operator==(a, b);
}

constexpr bool operator<(Rational a, Rational b) {
  return a.isDefined() and b.isDefined() and
         a.Nom() * b.Denom() < a.Denom() * b.Nom();
}

constexpr bool operator>(Rational a, Rational b) { return operator<(b, a); }

constexpr bool operator<=(Rational a, Rational b) {
  return not operator>(a, b);
}

constexpr bool operator>=(const Rational a, const Rational b) {
  return not operator<(a, b);
}

constexpr bool operator!(const Rational a) { return a.Nom() and a.isDefined(); }

inline Numeric fac(const Rational r) { return (::fac(r.toIndex())); }

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) {
  return Rational(a, 1) < b;
}

constexpr bool operator<(const Rational a, const Index b) {
  return a < Rational(b, 1);
}

constexpr bool operator>(const Index a, const Rational b) {
  return Rational(a, 1) > b;
}

constexpr bool operator>(const Rational a, const Index b) {
  return a > Rational(b, 1);
}

constexpr bool operator==(const Rational a, const Index b) {
  return a == Rational(b, 1);
}

constexpr bool operator!=(const Rational a, const Index b) {
  return not(a == b);
}

constexpr Rational abs(const Rational a) {
  return a < 0 ? -a : a;
}

constexpr Rational max(const Rational a, const Rational b) {
  return a < b ? b : a;
}  // Let other operators find out if this is allowed instead

typedef Array<Rational> ArrayOfRational;

constexpr Rational operator ""_2(unsigned long long int n) {
  return Rational(n, 2);
};

constexpr Rational operator ""_rat(unsigned long long int n) {
  return Rational(n, 1);
};

constexpr bool even(const Rational r) {
  if (r % 2)
    return false;
  else
    return true;
}

#endif  // rational_h