Approximating a decimal with a fraction (32-bit fixed point to two 23-bit numbers). Think binary, ease of computation.

Question

I have a 32-bit fixed-width decimal number between 0 and 1.0 (actually its guaranteed to be between 0.001 and 0.02, so loss of range is acceptable in the approximation). The binary representation is defined as most-significant-bit representing 1/2, bit after that 1/4, and so on to 32 bits.

I must represent this number as a ratio of two 23-bit numbers (actually one of them is 23 and one is 18 bits but let's simplify for the sake of discussion). It does not matter what the 2 numbers are as long as each of them fits in 23 bits.

Another way to think about this problem is I have a fraction with a 32-bit numerator and a known constant denominator of 2^32. I must approximate this fraction as another fraction with 23-bit numerator and 23-bit denominator, such that an error of no more than 1/(2^32) is acceptable.

Is there a reasonably efficient algorithm that would allow me to do this without doing a brute force guess-and-check? Some form of "bit twiddling" trick would be ideal as this must be done on a fairly small microcontroller. Barring that, the next best thing would be integer arithmetic.

I understand that since I am repacking 32 bits of information into a total 46 bits of information, theoretically it is possible. A brute force solution of guessing and checking until some pair of number works is pretty trivial. What I am looking for is some more efficient algorithm or at least a way of figuring out a good starting point for guess-and-check. Ideally, there would be a defined time bound for the execution of this algorithm.

For example, my input (A) is 0.005777551792562008. If we take the binary representation of this

0b 0000 0001 0111 1010 1010 0011 0011 1100

, it can be reinterpreted as the 32-bit unsigned integer 24,814,396.

This number can be represented as the ratio (N/D) $$ 0.005777551792562008 = \frac{24814396}{2^{32}} \approx \frac{41572}{7195435} = 0.005777552017355449 $$.

In this case, the error of the approximation is 2.248e-10 which is smaller than 1/(2^32)=2.328e-10 .

The situation described in this question is pretty specific because it's based on a real world scenario. However, please feel free to generalize, even if your generalization makes the answer not fit the specific situation I describe.

The point of this question is to get some ideas about where to start and hopefully leave something that is applicable in more situations. Computing a "close enough" fractional approximation of a decimal is a pretty general problem, but one I have not seen addressed specifically in the context of computers and binary representation.

Answer 1

Finding rational approximations to a number is a field known as Diophantine Approximation. In your case, you are trying to find the most accurate approximation within some limits on the values of the numerator and the denominator. It just so happens that you can find the best rational approximation to a number by using its Continued Fraction representation.

Algorithm

1. Take the integer part of the number and write it down.
2. Divide $1$ by the fractional part of the number.
3. Repeat steps $1$ and $2$ until the necessary precision is achieved.

I will give a quick example of finding the continued fraction for $\pi$.

$$\begin{align} \pi&=3+0.14159...\\ &=3+\frac1{7+0.06251...}\\ &=3+\frac1{7+\frac1{15+0.99659}}\\ \pi&\approx3+\frac1{7+\frac1{15+\frac11}}\\ &\approx\frac{355}{113} \end{align}$$

$\frac{355}{113}$ is the best rational approximation of $\pi$ with a denominator less than $113$. It is accurate up to fifteen digits. Generally, the error between the approximation and the number is $$\left|n-\frac pq\right|<\frac1{q^2}$$

This means that as long as you can get a denominator greater than $2^{16}=65536$, you’re guaranteed to get a result within the bounds you gave above.

Notes

Continued fractions are often written in the notation $$\pi\approx[3;7,15,1]= 3+\frac1{7+\frac1{15+\frac11}}$$ where the integer part is the first number and all the other numbers go in the denominator.

In your specific example of $0.005777551792562008$, the continued fraction of your approximation looks like

while the continued fraction approximation looks like

Notice how similar the two continued fractions are, with the first five terms being identical. However, the continued fraction approximation within the range you used is still more accurate, with your approximation having an error greater than the continued fraction error by $2.247823\cdot10^{-10}$.

Answer 2

Worked out example using continued fractions:

The fraction $24814396/2^{32}$ simplifies to $6203599/1073741824$, whose continued fraction is computed as follows:

$$\begin{align} 1073741824 &= \mathbf{173} \cdot 6203599 + 519197\\ 6203599 &= \mathbf{11} \cdot 519197 + 492432 \\ 519197 &= \mathbf{1} \cdot 492432 + 26765 \\ 492432 &= \mathbf{18} \cdot 26765 + 10662 \\ 26765 &= \mathbf{2} \cdot 10662 + 5441 \\ 10662 &= \mathbf{1} \cdot 5441 + 5221 \\ 5441 &= \mathbf{1} \cdot 5221 + 220 \\ 5221 &= \mathbf{23} \cdot 220 + 161 \\ 220 &= \mathbf{1} \cdot 161 + 59 \\ 161 &= \mathbf{2} \cdot 59 + 43 \\ 59 &= \mathbf{1} \cdot 43 + 16 \\ 43 &= \mathbf{2} \cdot 16 + 11 \\ 16 &= \mathbf{1} \cdot 11 + 5 \\ 11 &= \mathbf{2} \cdot 5 + 1 \\ 5 &= \mathbf{5} \cdot 1. \end{align}$$

So the (finite) continued fraction in this case is $[0;173,11,1,18,2,1,1,23,1,2,1,2,1,2,5]$.

Computing the convergents of this continued fraction (there is a very simple recursion for this computation) yields the sequence of approximations

173: 1/173 = 0.005780346820809248
 11: 11/1904 = 0.005777310924369748
  1: 12/2077 = 0.005777563793933558
 18: 227/39290 = 0.005777551539832018
  2: 466/80657 = 0.005777551855387629
  1: 693/119947 = 0.00577755175202381
  1: 1159/200604 = 0.005777551793583378
 23: 27350/4733839 = 0.005777551792530334
  1: 28509/4934443 = 0.005777551792573143
  2: 84368/14602725 = 0.005777551792559265
  1: 112877/19537168 = 0.00577755179256277
  2: 310122/53677061 = 0.005777551792561817
  1: 422999/73214229 = 0.005777551792562071
  2: 1156120/200105519 = 0.005777551792562004
  5: 6203599/1073741824 = 0.005777551792562008

We can see that the denominator 4934443 is small enough for a 23-bit integer but 14602725 is too large. So one answer would be to cut it off at $28509/4934443$ for an error of 1.113e-14. This is probably optimal because $4934443$ is quite close to the upper limit, but if the limit were just a little bit higher we could lower the coefficient $2$ on the next row to $1$ and get $55859/9668282$ which has a slightly better error of less than 1e-14.