# Diophantine approximation with additional constraints

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.

• I just asked a similar question over on the CS Stack Exchange and was excited to find your question! However, I was disappointed by the lack of answers... Commented Jul 1, 2021 at 23:40
• @MarkOmo You may be interested that there is an answer after almost 9 years. Commented Sep 23, 2023 at 2:14

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).

• Wow, it took almost 9 years before an answer was suitable enough. I myself have experienced something similar. I guess the problem is sometimes questions get buried. Commented Sep 23, 2023 at 2:11
• @Tito Piezas III This was a case of a solution staring me right in the face at the time I asked the question, but me failing to connect the dots. In my notes, I referred to the computation of $i, j$ from $c$ as a factorization! And before I worked on the implementation of Lambert $W$, where I needed to split a new constant, I did not revisit this issue. I felt so stupid when I realized that I had overlooked in 2014 what looks so obvious now. Commented Sep 23, 2023 at 2:51