rational.hpp

Go to the documentation of this file.

/*
  RATIONAL.hpp  -  support for precise rational arithmetics
 
   Copyright (C)
     2022,            Hermann Vosseler <Ichthyostega@web.de>
 
 **Lumiera** 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. See the file COPYING for further details.
 
*/
 
 
#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 (ullong num)
{
  return util::Rat{num};
}
 
 
#endif /*LIB_RATIONAL_H*/