# Implement 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 using common machine operations (fixed-size integers and floating point operations) 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
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
This is a dubious migration, since it's borderline off-topic here, too. The confusingly named Computational Science Stack Exchange, which covers scientific computation, seems most appropriate. –  David Richerby Nov 26 '13 at 9:26

## 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);
}
-
Please get rid of the source code and replace it with ideas, pseudo code and arguments of correctness. See here and here for related meta discussions. –  Raphael Mar 24 at 14:13
@Raphael: Go look at the comments under Yuval's answer, they are the reason I answered at all. In my post the ideas are above the code, and the code adds value that helped the OP (and will be helpful to anyone who reaches this post via google.) Note that the OP specifically asked for the helper functions that use C unions to cast back and forth to the bit-wise representation of a float. I don't know how to write that in any language other than C. And I don't think that C's union and casting facilities are "obvious" to someone who originally asked a question on mathematica.stackexchange. –  Wandering Logic Mar 24 at 14:42
I would think it is possible to discuss arbitrary precision reals without tying oneself to a programming language. If not, this question is probably offtopic here and should have been migrated so Stack Overflow, not here. (Bit late for that.) –  Raphael Mar 24 at 16:10
Yes; late for yet another migration. This is why I despise the stupid balkanization of related topics on stack exchange. OP asked in Mathematica, it got migrated here. David Richerby said it should be migrated to Computational Science, you say it should be migrated to SO. Meanwhile, I actually answered the OPs question (which is about machine precision reals, not arbitrary precision), but get critiqued for being long winded, or having bad style, or whatever. What's more important, perfection, or helping the OP get an adequate answer to the question they actually asked? –  Wandering Logic Mar 24 at 17:14
The integrity of any site's scope is important for several reasons (maybe just not for you?). If you want, we can discuss this in Computer Science Chat. –  Raphael Mar 24 at 20:50
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.