Rational Reconstruction
Since I started using exact precision arithmetic several people have asked me the following question.
''Given a floating point number, how do I figure out what number the computer really means?''
Most computations on computers are performed using floating point numbers. While fast for computation, they are not exact representations. If a program outputs the number .33333, it likely represents 1/3, but what if it outputs 0.42857?
It turns out there is a technique to recover the exact value of a rational number given an accurate approximation and bound on its denominator. We will explain how to recover a rational number and provide a simple code to do the calculation.
Basic Theory
More formally, we can state the above problem as follows; given a number r, find the rational number p/q with denominator at most B which is the closest to r. To solve this problem we can use continued fraction approximations. The continued fraction representation of r is defined by
The full continued fraction expansion of r gives its exact representation, and number defined by the first i terms of the expansion is called the ith convergent of r. The convergents provide increasingly good approximations of r, and will solve our problem. If we are given an accurate enough floating point approximation r of a rational number x, and know that the denominator of x is at most B, then x will appear as the last convergent of r having denominator at most B. The continued fraction convergents can be computed by running the Extended Euclidean Algorithm.
This technique of finding a low denominator rational approximation of a number has several names. It is referred to as diophantine approximation or rational roundoff. There is a directly related problem which is most often referred to as rational number reconstruction where instead of looking for a low denominator rational approximation of a number, we are looking for a low denominator rational number p/q satisfying qn= p mod M for given n, M. The method for solving this modular version is highly related and also uses the Extended Euclidean Algorithm.
For more on continued fractions there are several books but the wikipedia entry is a good place to start. For more on diophantine approximation I recommend Chapter 6 of Theory of Linear and Integer Programming by A. Schrijver, and for rational number reconstruction I recommend Section 5.10 of Modern Computer Algebra by J. von zur Gathen and J. Gerhard.
Reconstruction Code
|
Here the source code, it is written in C for Linux/UNIX. [Download] It requires the GNU Multiple Precision (GMP) library for multiple precision arithmetic. The text to the right demonstrates the command line interface. The code includes functions and instructions to solve both the diophantine approximation and rational number reconstruction problems. |
$ ./reconstruct 0.42857 Input is: 0.42857 float rep. is: 4.285700E-01 Approximation is: 3/7 $ |
Related Links
GNU Multiple Precision Arithmetic LibraryGMP is a high performance multiprecision arithmetic library. Our code relies on this library.
NTL - A Library for doing Number Theory
A high-performance, portable C++ library providing data structures and algorithms.