2
$\begingroup$

The Problem

I am working on a problem that boils down to finding the closest representation of an arbitrary number ($x$) in the form:

$$x = A\times\frac{N}{D}$$

Where $A$ is a 32-bit integer, and $N$ and $D$ are 16-bit integers.

Typical inputs are real numbers between 1 and 1 billion.


My Attempts

I first attempted to use $N$ and $D$ to represent the fractional part of the number, but I can only correct for errors less than one part in 65 thousand, not very useful.

I then looked at naively testing all the possible combinations of $N$ and $D$ to find one with the lowest remainder, even when skipping like fractions (i.e., 2/4, 3/6, etc.), and skipping values of $N$ that would result in a fraction that is too large ($N \leq \frac{x}{D}$). I still need to search more than a billion combinations.

Some test cases I have found challenging; clearly, it's pretty easy to come up with these.

  • $1000\times \frac{255}{256}$
  • $ 123456 \times \frac{123}{456}$

Backround

This problem comes up as I am trying to configure a PLL to convert arbitrary input frequencies to a uniform output frequency with the lowest error possible.

Question

Is there a much more efficient way to solve this problem? Does this kind of problem have a name that I can look up?

$\endgroup$
3
  • 1
    $\begingroup$ How is the real number given as input in your problem? An input to an algorithm should be a finite string. $\endgroup$
    – Vincenzo
    Jul 1 at 21:28
  • $\begingroup$ @Vincenzo Well, in my application, it would be a double-precision float. Here I use the term "Algorithm" in its colloquial sense. I also reviewed the tag description to see if this question was a good fit, "An algorithm is a sequence of well-defined steps that defines an abstract solution to a problem. Use this tag when your issue is related to design and analysis of algorithms". Based on that, I thought the tag was a good fit. $\endgroup$
    – Mark Omo
    Jul 1 at 21:35
  • $\begingroup$ I am looking for an algorithm that is better (ideally substantially) than brute force for finding the optimal set of integers $A$, $N$, $D$ that gives the lowest error. $\endgroup$
    – Mark Omo
    Jul 1 at 21:49
3
$\begingroup$

Generally speaking, this problem is called diophantine approximation.

However, your inputs are not just real numbers but double-precision floating point (aka double) numbers, which are in fact rational. Any double number has the form:

$$ (-1)^{sign} ( 1 . b_{51} b_{50} ... b_{0}) \times 2^{e - 1023} $$ where $b$ encodes the mantissa (52 bits) and $e'=e-1023$ encodes the exponent (11 bits). You mention that your typical inputs are at least 1, so your exponent $e'$ is always nonnegative. The mantissa is always an integer multiple of $2^{-52}$, and will be an integer multiple of $2^{-\ell}$ whenever the least significant nonzero bit of the mantissa is the $\ell$-th one.
The bottom line is that to achieve an exact representation you could do the following:

  • look at the binary representation of your floating point number $x$, isolating the mantissa bit sequence $b$ and the exponent $e'=e-1023$;
  • let $\ell$ be the index of the least significant nonzero bit of the mantissa $b$;
  • your number $x$ is an integer multiple of $1/D = 2^{e' - \ell}$, and no fraction with a smaller denominator than $D$ can represent it exactly;
  • divide $x$ by $1/D = 2^{e'-\ell}$ to find $N \times A$.

In the worst case, this may however not guarantee that $N$ and $A$ satisfy your size requirement, since $\ell- e'$ could be up to $52$, while you only have $32 + 16 = 48$ bits available for $N \times A$. A practical solution (but perhaps non-optimal) would be to round the least significant $52-48 = 4$ bits in the mantissa of $x$.

If you insist on finding the best approximation given a constraint on the largest denominator $D$, we need to turn to diophantine approximation techniques and continued fractions.

A fraction $A/B$ is a best rational approximation to real number $x$ if $|A - Bx| \le |C - Dx|$ for all integers $C$ and $D$ where $0 < D \le B$.

Then it is known that the best rational approximation of $x$ is a so-called convergent to $x$, which you can recover from the continued fraction expansion of $x$: $$ x = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \ldots}} $$ The $m$-th convergent $C_m$ is defined by the fraction $A_m/B_m$ where $A_k$ and $B_k$ are defined recursively by

$A_{-1}=1$, $A_0=a_0$, $A_{k} = a_{k} A_{k-1} + A_{k-2}$

$B_{-1}=0$, $B_0=1$, $B_{k} = a_{k} B_{k-1} + B_{k-2}$.

In other words, for each $k=0,1,2,\ldots$:

  • generate $a_k$ by setting $a_k = \lfloor x_k \rfloor$, $x_{k+1} = 1/(x_k - \lfloor x_k \rfloor)$, where $x_0 = x$;
  • compute $A_k$ and $B_k$ with the formulas above
  • the $k$-th convergent is $A_k/B_k$.

This will give you a sequence of closer and closer approximations of $x$, each of which is a best rational approximation. In your case, you want to stop when $B_k$ does not fit anymore into $16$ bits.

In your examples, we get the following:

First example.

$x = 1000 \times 255 / 256 \approx 996.09375$

$A_0 / B_0 = 996 / 1 \approx 996.0$

$A_1 / B_1 = 9961 / 10 \approx 996.1$

$A_2 / B_2 = 10957 / 11 \approx 996.0909090909091$

$A_3 / B_3 = 20918 / 21 \approx 996.0952380952381$

$A_4 / B_4 = 31875 / 32 \approx 996.09375$

Second example.

$x = 123456 \times 123 / 456 \approx 33300.63157894737$

$A_0 / B_0 = 33300 / 1 \approx 33300.0$

$A_1 / B_1 = 33301 / 1 \approx 33301.0$

$A_2 / B_2 = 66601 / 2 \approx 33300.5$

$A_3 / B_3 = 99902 / 3 \approx 33300.666666666664$

$A_4 / B_4 = 266405 / 8 \approx 33300.625$

$A_5 / B_5 = 632712 / 19 \approx 33300.63157894737$

$\endgroup$
1
  • $\begingroup$ Thanks! This seems to work well for my needs in practice. However, for around 0.08% of inputs, it does not find an optimal solution. For example, with the input $795082588 \times (1763/65225)$, the result is $988571985/46$, with an error of $1.5 \times 10^{-14}$. You can see my implementation here. $\endgroup$
    – Mark Omo
    Jul 3 at 0:39
1
$\begingroup$

The "best" rational approximations to a real number are given by it's continuous fraction, check that out.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.