Diophantine approximation with additional constraints

Question

I am trying to compute best rational approximations to various transcendental numbers $c$, subject to the following constraints: $$\frac {i j} {2^k} = c + \epsilon, \space\space2^n \le i, j \lt 2^{n+1}$$

For sufficiently small $n$, say $n \lt 40$, optimal values of $i,j,k$ which minimize $\vert\epsilon\vert$ can be found by exhaustive search. For example, for $n=23$ and $c=\pi$, the best rational approximation under my constraints is provided by $i=14120171, j=15656321, k=46$ such that $$\frac {14120171 \times 15656321}{2^{46}} = \pi + 3.1974({10^{-14}})$$

How would I go about finding optimal rational approximations for larger $n$, say $n=52$, for which exhaustive search is not feasible?

For large $n$ I am currently using heuristic searches based on methods like simulated annealing, but I am wondering whether there is a mathematical approach rather than an engineering approach I could utilize. I am aware that one can use the convergents of the continued fraction expansion of a transcendental number to determine best rational approximations, but I do not know how to then extend this approach to rational approximations that satisfy my constraints.

The practical relevance of my constraints is that for an appropriate choice of $n$, $a=\frac {i}{2^{\left \lceil \frac {k}{2} \right \rceil}}$ and $b=\frac {j}{2^{\left \lfloor \frac {k}{2} \right \rfloor}}$ are exactly representable as IEEE-754 floating-point numbers.

Modern floating-point processors frequently offer a fused multiply-add operation (FMA), which computes $ab+c$ such that the exact product $ab$ is used in the addition, and only the final sum is rounded. This operation is also exposed as a function in programming languages, for example the standard math functions fma() and fmaf() in ISO C99.

In conjunction with the desired rational approximations, this allows a more accurate computation of expressions like $\pi-x$ which can be programmed as fmaf(1.68325555f, 1.86637890f, -x) instead of simply using 3.14159265f - x.

Answer 1

I recently had to revisit this issue for an implementation of the primary branch of the real-valued Lambert $W$ function, $W_{0}(\cdot)$, for which I needed a highly accurate approximation to $e^{-1}$. For IEEE-754 double-precision computation, $n=52$. Since $2^{-2} \lt e^{-1} \lt 2^{-1}$, $k= 106$, so I was looking for $2^{52} \le i, j \lt 2^{53}$, such that $$\frac{i j}{2^{106}} = e^{-1} +\epsilon$$

Let $\lfloor x\rceil$ denote the nearest integer to $x$. $n = \lfloor 2^{k}e^{-1} \rceil = 29845926042406685857117349218881$. Factor $n, n-1, n+1, n-2, n+2, \ldots$, and check whether the factors found can be combined into $i, j$ that satisfy the requirements. Luckily, $n+1 = 29845926042406685857117349218882$ factors as follows: $2 \times 23 \times 71567107 \times 81105469 \times 111779876380849$. By suitable grouping of the factors, $i = 71567107 \times 81105469 = 5804483778208183$, and the remaining factors are then combined into $j = 2 \times 23 \times 111779876380849 = 5141874313519054$. This gives

$$ \frac{5804483778208183 \times 5141874313519054}{2^{106}} = e^{-1} + 8.2867 \cdot 10^{-33}$$

I tried this integer factorization approach with a few constants $c$ and noticed that a neighbor of $n$ with suitable factorization can be found among a fairly limited number of candidates. While this approach still requires a search, it can be concluded much faster than a brute-force search (which is now feasible, typically finding a solution within 250 to 350 hours on one core of a fast workstation).