# How could I implement, in C++, Mathematica's capability of rationalizing machine reals?

If I have a variable x bound to a machine precision real in Mathematica, I can use

 y = FromDigits[RealDigits[x]]


then y is bound to a completely equivalent rational number which has infinite precision.

What are the basic principles behind

 FromDigits


and

 RealDigits


How could I implement a similar conversion in C++ to coerce a double in a rational, or at least obtain the long integer denominator and numerator of a double?

-
You may find the wiki page on best rational approximation helpful. –  Sasha Nov 26 '13 at 6:16
This question appears to be off-topic because it is not directly related to Mathematica. –  Sjoerd C. de Vries Nov 26 '13 at 6:40
Only Mathematica can obtain the desired conversion result for my prolbem, I also asked similar question here: stackoverflow.com/questions/19808261/… , however, Matlab and C++ performs similar. So I forward similar question here to see whether there is possible C++ solution from Mathematica users. –  LCFactorization Nov 26 '13 at 6:52
The question is not about how to do something in Mathematica, but how to do something that you can do in Mathematica in another programming language. That is clearly off-topic here. –  Sjoerd C. de Vries Nov 26 '13 at 7:15

## migrated from mathematica.stackexchange.comNov 26 '13 at 9:21

This question came from our site for users of Mathematica.

On most current machines a double-precision real number will be represented in IEEE 754's binary64 format.

  63 62 x x x x x x 52 51 x x x x x x x x x x x x x x x x x x x x x 0
| s |      exp        |                mantissa                      |


s=0 implies positive, s=1 implies negative. exp is an 11 bit biased signed integer in the range [-1023, 1024] (take the unsigned number from the bits and subtract 1023). mantissa stores the bottom 52 bits of the fraction.

Numbers with exponents in the range [-1022, 1023] are normalized (bit 52 of the mantissa is implicitly "1".) The number represented by a particular bit pattern is $-1^s (1 + (\mathrm{mantissa} / 2^{52})) 2^\mathrm{exp}$.

Numbers with exp=1024 are special: if the mantissa is 0 then they represent + or - Infinity (from overflow or divide by 0, for example). If the mantissa is non-zero then the represent NaN (for example sqrt(-1)).

Numbers with exp=-1023 are special: if the mantissa is 0 then the represent + or - 0. If the mantissa is non-zero then the represent denormalized numbers: $-1^s (\mathrm{mantissa}/2^{1075})$.

In C++ on most current machines the following code will probably do most of what you want (I adapted this from some binary32 code without testing, so test first.) We will use uint64_t from the cstdint header for the bit-level representation, and a union for extracting the bits. The following code does not deal with zero, denormals, NAN or positive or negative Infinity (HUGE_VAL).

#include <cstdint>    // for uint64_t

typedef union
{
double   value;
uint64_t bits;
} ieee754_binary64_union;

#define IEEE754_MANTISSA_BITS 52
#define IEEE754_EXPONENT_BITS 11

#define IEEE754_HIDDEN_BIT    (((uint64_t)1) << IEEE754_MANTISSA_BITS)
#define IEEE754_EXPONENT_MASK ((1 << IEEE754_EXPONENT_BITS) - 1)
#define IEEE754_EXPONENT_BIAS ((1 << (IEEE754_EXPONENT_BITS - 1)) - 1)

static inline
double
ieee754_bits2float(uint64_t u)
{
ieee754_binary64_union fiu;
fiu.bits = u;
return fiu.value;
}

static inline
uint64_t
ieee754_float2bits(double f)
{
ieee754_binary64_union fiu;
fiu.value = f;
return fiu.bits;
}

static inline
int
ieee754_get_sign(double f)
{
return ieee754_float2bits(f) >> 63;
}

static inline
int
ieee754_get_exponent(double f)
{
uint64_t bits = ieee754_float2bits(f);
return (((bits >> IEEE754_MANTISSA_BITS) &
IEEE754_EXPONENT_BIAS);
}

static inline
uint64_t
ieee754_get_mantissa(double f)
{
uint64_t bits = ieee754_float2bits(f);
return ((bits & IEEE754_MANTISSA_MASK) + IEEE754_HIDDEN_BIT);
}

static inline
double
ieee754_make_double(int      sign,     // 1 (neg) or 0 (pos)
int      exponent, // range -1023 -> 1024
uint64_t mantissa) // with or without hidden bit
{
uint64_t collected_bits =
((uint64_t)sign << (IEEE754_MANTISSA_BITS + IEEE754_EXPONENT_BITS))     |
((uint64_t)((exponent + IEEE754_EXPONENT_BIAS) & IEEE754_EXPONENT_MASK) <<
IEEE754_MANTISSA_BITS)                                      |
return ieee754_bits2float(collected_bits);
}

-
The continued fraction algorithm is easy enough to implement. The first step is to compute the continued fraction of the input $x = [c_0;c_1,\ldots]$. You start with $x_0 = x$, and use the recurrence $c_i = \lfloor x_i \rfloor$, $x_{i+1} = 1/(x_i - c_i)$. You stop when $x_i - c_i$ is "small enough". The next step is to compute the convergent of the continued fraction. You can find the appropriate formulas on the web. In practice, both steps are performed at the same time.
The resulting rational number is not actually equivalent to the floating point number - there's probably some small error. In fact, the rational approximation may be more accurate in some cases, and off the mark in others. However, floating point numbers are secretly rational numbers - they are of the form $\pm 2^x M$, where $x$ is the integer exponent (positive, zero or negative) and $M$, the mantissa, is a rational number between $0$ and $1$ (or $0$ and $2$). In fact, it's always of the form $A/2^B$, where $B$ is the length of the mantissa and $A$ is some non-negative integer. If you really care for it, you could decode a real number into this representation, and thus obtain a truly equivalent rational number.