rational.hpp

Go to the documentation of this file.

/*
  RATIONAL.hpp  -  support for precise rational arithmetics

  Copyright (C)         Lumiera.org
    2022,               Hermann Vosseler <Ichthyostega@web.de>

  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 of
  the License, 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., 675 Mass Ave, Cambridge, MA 02139, USA.

*/


#ifndef LIB_RATIONAL_H
#define LIB_RATIONAL_H


#include "lib/integral.hpp"
#include "lib/util-quant.hpp"

#include <cmath>
#include <limits>
#include <stdint.h>
#include <boost/rational.hpp>


namespace util {
  
  using Rat = boost::rational<int64_t>;
  using boost::rational_cast;
  using std::abs;
  
  
  
  inline bool
  can_represent_Product (int64_t a, int64_t b)
  {
    return ilog2(abs(a))+1
         + ilog2(abs(b))+1
         < 63;
  }
  
  inline bool
  can_represent_Product (Rat a, Rat b)
  {
    return can_represent_Product(a.numerator(), b.numerator())
       and can_represent_Product(a.denominator(), b.denominator());
  }
  
  inline bool
  can_represent_Sum (Rat a, Rat b)
  {
    return can_represent_Product(a.numerator(), b.denominator())
       and can_represent_Product(b.numerator(), a.denominator());
  }
  
  
  inline int64_t
  reQuant (int64_t num, int64_t den, int64_t u)
  {
    u = 0!=u? u:1;
    auto [d,r] = util::iDiv (num, den);
    // round to smallest integer fraction, to shake off "number dust" 
    f128 const ROUND_ULP = 1 + 1/(f128(std::numeric_limits<int64_t>::max()) * 2);
    
    // construct approximation quantised to 1/u
    f128 frac = f128(r) / den;
    int64_t res = d*u + int64_t(frac*u * ROUND_ULP);
    ENSURE (abs (f128(res)/u - rational_cast<f128>(Rat{num,den})) <= 1.0/abs(u)
           ,"Requantisation error exceeded num=%li / den=%li -> res=%li / quant=%li"
           ,                                   num,       den,       res,        u);
    return res;
  }
  
  inline Rat
  reQuant (Rat src, int64_t u)
  {
    return Rat{reQuant (src.numerator(), src.denominator(), u), u};
  }
} // namespace util



inline util::Rat
operator""_r (unsigned long long num)
{
  return util::Rat{num};
}


#endif /*LIB_RATIONAL_H*/