// RationalApproximation.h
//---------------------------------------------------------------------------
#ifndef __RATIONALAPPROXIMATION_H
#define __RATIONALAPPROXIMATION_H

//---------------------------------------------------------------------------
class TRationalApproximation
{
  public:
    TRationalApproximation();

    // Modifiers
    void Calculate(double value);    
    void Calculate(double value, double tolerance);

    // Selectors
    inline double GetNumerator() const { return m_numerator; }
    inline double GetDenominator() const { return m_denominator; }

  protected:
    const double m_defaultTolerance;
    double m_numerator;
    double m_denominator;
};

//---------------------------------------------------------------------------
#endif

// RationalApproximation.cpp
//---------------------------------------------------------------------------
#include "RationalApproximation.h"
#include <math.h>
#include <algorithm>
#include <limits> 

//---------------------------------------------------------------------------

TRationalApproximation::TRationalApproximation()
  : m_defaultTolerance(1e-6)
  , m_numerator(0.0)
  , m_denominator(0.0)
{
}
//---------------------------------------------------------------------------

void TRationalApproximation::Calculate(double value)
{
  Calculate(value, m_defaultTolerance * fabs(value));
}
//---------------------------------------------------------------------------

void TRationalApproximation::Calculate(double value, double tolerance)
{
  double x = value,
         c11 = 1.0, c12 = 0.0,
         c21 = 0.0, c22 = 1.0,
         epsilonTolerance = std::numeric_limits<double>::epsilon()*fabs(value),
         maxTolerance = std::max(tolerance, epsilonTolerance);

  while(1)
  {
    double d = (x > 0) ? floor(x + 0.5) : -floor(-x + 0.5); // round(x)
    x -= d;

    double tmp_c12 = c12;
    double tmp_c22 = c22;

    c12 = c11;
    c22 = c21;
    c11 = c11 * d + tmp_c12;
    c21 = c21 * d + tmp_c22;

    if(x == 0 || fabs(c11/c21 - value) <=  maxTolerance)
    {
      break;
    }

    x = 1.0 / x;
  }

  m_numerator = (c21 > 0) ? c11 : -c11;
  m_denominator = fabs(c21);
}
//---------------------------------------------------------------------------

#include <iostream>
#include <vector>

#include "RationalApproximation.h"

int main(int argc, char* argv[])
{
  typedef std::vector<double> values_type;
  values_type vals;
  vals.push_back(M_PI);
  vals.push_back(0.1);
  vals.push_back(0.02);
  vals.push_back(0.003);
  vals.push_back(0.0004);
  vals.push_back(0.00005);
  vals.push_back(0.000006);
  vals.push_back(0.0000007);
  vals.push_back(0.00000008);
  vals.push_back(0.000000009);
  vals.push_back(247037.0 / 4.0);
  vals.push_back(1.999999999);
  vals.push_back(0.12345678901);
  
  TRationalApproximation rat;

  std::cout.unsetf(std::ios::skipws);
  std::cout.precision(std::numeric_limits<values_type::value_type>::digits10 + 1);

  const double tol = 1e-6;
  std::cout << "Tolerance = " << tol << std::endl;

  for(values_type::const_iterator it = vals.begin(); it != vals.end(); ++it)
  {
    rat.Calculate(*it, tol);

    std::cout << *it << " = ";
    std::cout << rat.GetNumerator() << " / " << rat.GetDenominator() << std::endl;
  }

  return 0;;
}