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$
