BestRationalApproximation.h

Go to the documentation of this file.

/* This file is part of gbLAB.
 *
 * gbLAB is distributed without any warranty under the MIT License.
 */
 
 
#ifndef gbLAB_BESTRATIONALAPPROXIMATION_H_
#define gbLAB_BESTRATIONALAPPROXIMATION_H_
 
#include <math.h>
 
// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
 
namespace oILAB {
 
class BestRationalApproximation {
public:
  // typedef int64_t LongIntType;
  typedef long long int LongIntType;
 
private:
  static std::pair<LongIntType, LongIntType> rat_approx(double f,
                                                        LongIntType md) {
    if (fabs(f) < 1.0 / md) {
      return std::pair<LongIntType, LongIntType>(0, 1);
    }
 
    /*  a: continued fraction coefficients. */
    LongIntType a, h[3] = {0, 1, 0}, k[3] = {1, 0, 0};
    LongIntType x, d, n = 1;
    long int i, neg = 0;
 
    if (md <= 1) { // *denom = 1; *num = (LongIntType) f; return;
      return std::pair<LongIntType, LongIntType>((LongIntType)f, 1);
    }
 
    if (f < 0.0) {
      neg = 1;
      f = -f;
    }
 
    while (f != floor(f)) {
      n <<= 1;
      f *= 2;
    }
 
    d = f;
 
    /* continued fraction and check denominator each step */
    for (i = 0; i < 64; i++) {
 
      a = n ? d / n : 0;
      if (i && !a)
        break;
 
      x = d;
      d = n;
 
      n = x % n;
 
      x = a;
      if (k[1] * a + k[0] >= md) {
        x = (md - k[0]) / k[1];
        if (x * 2 >= a || k[1] >= md)
          i = 65;
        else
          break;
      }
 
      h[2] = x * h[1] + h[0];
      h[0] = h[1];
      h[1] = h[2];
      k[2] = x * k[1] + k[0];
      k[0] = k[1];
      k[1] = k[2];
    }
    return std::pair<LongIntType, LongIntType>(neg ? -h[1] : h[1], k[1]);
  }
 
  const std::pair<LongIntType, LongIntType> intPair;
 
public:
  const LongIntType &num;
  const LongIntType &den;
 
  BestRationalApproximation(const double &f, const LongIntType &md)
      : /* init list */ intPair(
            isfinite(f) ? rat_approx(f, md)
                        : std::pair<LongIntType, LongIntType>(1, 0)),
        /* init list */ num(intPair.first),
        /* init list */ den(intPair.second) {}
        
    };
 
    } // namespace oILAB
#endif