Converting from floating point to rational

//
// If FLT_RADIX == 2 as expected, we think we can do floating point to rational
// conversions exactly from the exponent and mantissa, otherwise we'll fall back
// on continued fractions.
//
#if FLT_RADIX != 2 && !defined(RATIONAL_USE_CONTINUED_FRACTIONS)
  #define RATIONAL_USE_CONTINUED_FRACTIONS
#endif
#ifdef RATIONAL_USE_CONTINUED_FRACTIONS
  #ifndef RATIONAL_FLOATING_POINT_ACCURACY
    #define RATIONAL_FLOATING_POINT_ACCURACY 0.00001
  #endif
#endif

The idea is that, if we happen to be compiling for something like an IBM mainframe (where, IIRC, FLT_RADIX is 16), then we’ll stick with the good old continued fractions algorithm rather than trying to get my tricky code to work.  In this case, we’ll need to exit a loop when the value is “close enough”; and the …ACCURACY macro provides a default of ±0.001% of the passed value.  (I had no reason to pick that number; I just pulled it out from my rear end.)

Users can define either or both of those macros for themselves if they want to force the continued fractions algorithm or provide a different accuracy default.

Excerpts from the rational class:

class rational final
{
    bigint num, den;

  #ifdef RATIONAL_USE_CONTINUED_FRACTIONS
    void init(long double, double);
  #else
    void init(long double);
  #endif

    // ...

public:

Explicit construction uses SFINAE to make sure that these functions apply only to floating point arguments.

  #ifdef RATIONAL_USE_CONTINUED_FRACTIONS
    template<typename F>
    explicit rational(F val, double acc = RATIONAL_FLOATING_POINT_ACCURACY,
      typename std::enable_if<std::is_floating_point<F>::value,int>::type = 0)
      : num(), den()
    {
        init(long double(val), acc);
    }
  #else
    template<typename F>
    explicit rational(F val,
      typename std::enable_if<std::is_floating_point<F>::value,int>::type = 0)
      : num(), den()
    {
        init(long double(val));
    }
  #endif

Same idea for explicit assignment:

  #ifdef RATIONAL_USE_CONTINUED_FRACTIONS
    template<typename F>
    rational& assign(F val, double acc = RATIONAL_FLOATING_POINT_ACCURACY,
      typename std::enable_if<std::is_floating_point<F>::value,int>::type = 0)
    {
        init(long double(val), acc);
        return *this;
    }
  #else
    template<typename F>
    rational& assign(F val,
      typename std::enable_if<std::is_floating_point<F>::value,int>::type = 0)
    {
        init(long double(val));
        return *this;
    }
  #endif

    // ...
};

And in rational.cpp:

#ifdef RATIONAL_USE_CONTINUED_FRACTIONS

void rational::init(long double val, double accuracy)
{
    if (std::isnan(val) || std::isinf(val))
    {
        throw std::invalid_argument(
            "Converting non-finite floating point value to rational");
    }

    long double max_error = std::fabs(val * accuracy);

    bool neg = val < 0.0;
    if (neg)
    {
        val = -val;
    }

    //
    // continued fractions (Knuth, 4.5.3):
    //

    long double val0 = val; // current fraction term

    bigint int0 = bigint(val0);

    bigint num0 = int0,  // current approximation
           den0 = 1U;    // first guess:  int(val) / 1

    bigint num1 = 1U, num2; // previous and 2nd-previous numerator
    bigint den1 = 0U, den2; // previous and 2nd-previous denominator

    long double err0 = std::fabs(val0 - long double(int0));

    while (err0 > max_error)
    {
        val0 = 1.0 / (val0 - long double(int0));
        int0 = bigint(val0);

        num2 = num1;
        num1 = num0;
        num0 = int0 * num1 + num2;

        den2 = den1;
        den1 = den0;
        den0 = int0 * den1 + den2;

        long double err1 = err0;
        err0 = std::fabs(val - long double(num0) / long double(den0));

        if (err0 >= err1) // oops...no longer converging
        {
            throw std::logic_error(
                "Overflow converting floating point to rational");
        }
    }

    num = std::move(num0);
    den = std::move(den0);
    if (neg)
    {
        num.negate();
    }
}

#else

void rational::init(long double val)
{
    if (std::isnan(val) || std::isinf(val))
    {
        throw std::invalid_argument(
            "Converting non-finite floating point value to rational");
    }

    if (val == 0.0)
    {
        num.clear();
        den.set_to_one();
        return;
    }

    bool neg = val < 0.0;
    if (neg)
    {
        val = -val;
    }

    union
    {
        long double ldval;
        unsigned char ucval[sizeof(long double)];
    } mant;

    int exp;
    mant.ldval = std::frexp(val, &exp);

    //
    // Find the mantissa bits:
    //
    static constexpr unsigned mant_msb = // index for mant.ucval
        sizeof(long double) - LDBL_MANT_DIG / CHAR_BIT;
    static constexpr unsigned phantom_bit = LDBL_MANT_DIG % CHAR_BIT;

    //
    // Put the most significant 1 back in:
    //
    mant.ucval[mant_msb] |= 1 << phantom_bit;

    //
    // And mask off the higher bits just in case:
    //
    unsigned char mask = 1;
    for (unsigned bit = phantom_bit; bit > 0; --bit)
    {
        mask |= 1 << bit;
    }
    mant.ucval[mant_msb] &= mask;

    //
    // Make a fraction equal to just the raw mantissa bits:
    //
    num.clear();
    for (unsigned idx = mant_msb; idx < sizeof(long double); ++idx)
    {
        num <<= CHAR_BIT;
        num += mant.ucval[idx];
    }
    den.set_to_one();

    if (exp < 0)
    {
        den <<= -exp;
    }
    else if (exp > 0)
    {
        num <<= exp;
    }

    lowest_terms();
    if (neg)
    {
        num.negate();
    }
}

#endif // not RATIONAL_USE_CONTINUED_FRACTIONS

And just in case you’re wondering what the lowest_terms() member function is about:

rational& rational::lowest_terms()
{
    bigint gcdval(std::move(dbacc::gcd(num, den)));
    if (gcdval > 1)
    {
        num /= gcdval;
        den /= gcdval;
    }
    return *this;
}