# Computing all “suboptimal” rational approximations to $\pi/2$

I have an irrational number $\alpha$ ($\alpha=\frac\pi2$), and I would like to determine all integers $n\in[1,N]$ ($N=10^{16}$) that satisfy $$n \epsilon(n)^2 \leq \tau$$ where $\tau$ is a known real number ($\tau=78$), and $\epsilon(n)$ is the distance between $n$ and the closest multiple of $\alpha$: $$n = m \alpha + \epsilon(n), \qquad |\epsilon(n)| < \tfrac12\alpha.$$

I know that if I approximate $\alpha$ by a rational number $A/B$, then the problem can be solved with extended Euclid's algorithm applied to $$Bn-Am = B\epsilon,$$ where $B\epsilon$ is an integer, with an upper bound derived from $n\epsilon^2<\tau$, $B\epsilon<B\sqrt{\tau}$. For each value of $B\epsilon$ I can compute corresponding solutions $n\in[1,N]$.

The problem is that if $A/B$ is a sufficiently accurate rational approximation to $\alpha$ (i.e., $B\sim N$ ensures that $n\epsilon^2$ is approximated accurately), then the bound $B\epsilon < B\sqrt\tau$ (e.g., $B\sqrt\tau=10^{17}$) leads me to consider prohibitively many different $n$.

The standard theory of rational approximations tells me how to compute the best rational approximations with continued fractions, but I need to compute all rational approximations $n/m$ that are good enough in the above sense. If I write $$\alpha = \frac nm-\frac \epsilon m = \frac nm+\delta,$$ the equivalent problem is to compute all $(n,m)$ with $$\big|\alpha - \frac nm\big| < \frac{\sqrt\tau}{m \sqrt{n}} \sim \frac{\sqrt{\tau/\alpha}}{m^{3/2}},$$ but the right-hand side here is quite different from $1/m^2$ (or $1/m^\mu$ where $\mu\geq2$ is the irrationality measure of $\alpha$). Since $\tfrac32\leq\mu$, that tells me there are infinitely many solutions, but I don't know how to actually compute them.

Is there an efficient algorithm for computing $n$ with $n\epsilon^2<\tau$?

• How efficient do you expect this algorithm to be? There will be approximately $N^{1/2}$ solutions to $|m\alpha - n| < N^{-1/2}$ with $m,n < N$, so magically listing them all will already take $N^{1/2}$ steps. – Felipe Voloch Oct 11 '14 at 16:42
• Yes, it's possible to list all such $(m,n)$ pairs efficiently (they're basically the small linear combinations of consecutive best approximations); but [as I see Felipe Voloch commented while I was writing this] there will be about $\tau N^{1/2}$ of them, which for $(N,\tau) = (10^{16},78)$ is almost $10^{10}$ $-$ is that really what you want? – Noam D. Elkies Oct 11 '14 at 16:44
• @NoamD.Elkies $10^{10}$ is all right, I would consider that feasible; it's $10^{16}$ that's too large. – Kirill Oct 11 '14 at 16:59
• @FelipeVoloch How do you know how many solutions there should be? If there is an optimal algorithm that takes $O(\text{number of solutions})$, then I would like to know it. – Kirill Oct 11 '14 at 17:01
• Are you trying to evaluate numerically something like $\sum_{n=1}^N \sin^n n$?... – Noam D. Elkies Oct 12 '14 at 14:48

## 1 Answer

As noted in the comments (by Felipe Voloch and myself), one expects about $O_\tau(N^{1/2})$ solutions, because that's the area $|R|$ of the region $$R = \{ (m,n) \in {\bf R}^2 \colon 1 \leq n \leq N, \ (\alpha m - n)^2 < \tau/n \}$$ whose intersection with ${\bf Z}^2$ we want to list. (In the comment I got the dependence on $\tau$ wrong: it's $(\tau N)^{1/2}$, not $\tau N^{1/2}$, because the upper bound on $|\alpha m - n|$ is $\sqrt{\tau/n}$, not $\tau\,/\sqrt n$ as I misread; so for $(N,\tau) = (10^{16},78)$ we expect about $10^9$ pairs $(m,n)$, not $10^{10}$.)

To list $R \cap {\bf Z}^2$ in $\,\tilde{\!O}(|R|)$ time:

i) Partition $R$ into about $\log_2 N$ parts $R_1,R_2,\ldots R_k$ where $R_j = \{(m,n) \in R \colon 2^{-j} N \leq n < 2^{1-j} N$. The $R_j$ are all equivalent under affine transformations (except for $R_k$ [unless $N$ is exactly a power of $2$], but $R_k$ is negligible anyway). Then do the remaining steps for each $j$.

ii) Inscribe $R_j$ in an ellipse $E_j$ whose area is $O(|R_j|)$. For example $E_j$ might be $$\{ (x,y) \in {\bf R}^2 \colon \frac{4^{j+1}}{N^2}\Bigl(x - \frac32 2^{-j} N\Bigr)^2 + \frac {2^{j-1}\tau}{N} (\alpha x - y)^2 \leq 2 \}$$ if I did the algebra right.

iii) Choose an affine transformation $T_j$ that takes $E_j$ to the unit circle $D_1: X^2 + Y^2 \leq 1$.

iv) Apply $2$-dimensional lattice basis reduction (a.k.a.\ the Euclidean algorithm) to $T_j {\bf Z}^2 - T_j (0,0)$. Use the reduced basis to efficiently list $D_1 \cap T_j {\bf Z}^2$, and thus $E_j \cap {\bf Z}^2$.

v) Check whether each of the resulting integer vectors $(m,n)$ is actually in $R_j$, and output those that pass this test (which should be roughly the constant fraction $|R_j| \, / \, |E_j|$), QEF.