Detailed mathematical approach

Problem Summary

For each non-square integer \(d \lt 100\), we approximate \(\pi\) by numbers of the form

$a+b\sqrt d,\qquad |a|,|b|\le 10^{13}.$

Among all such pairs \((a,b)\), we choose the one that minimizes

$\left|\pi-\left(a+b\sqrt d\right)\right|.$

The quantity called the integral part is the coefficient \(a\), not the floor of the real number \(a+b\sqrt d\). The required output is therefore

$\sum_{d\in S} |a_d|,$

where \(S\) is the set of non-square integers from \(2\) to \(99\), and \(a_d\) is the integer coefficient belonging to the best approximation for that specific \(d\).

Mathematical Approach

Write \(\alpha=\sqrt d\) and \(\tau=\pi-3\). The central observation is that once \(b\) is fixed, the optimal \(a\) is forced to be the nearest integer to \(\pi-b\alpha\). So the real difficulty is to find the right \(b\) without scanning all \(2\cdot 10^{13}+1\) possibilities.

Step 1: Reduce the problem to a nearest-integer search

For fixed \(b\), define

$x_b=\pi-b\alpha.$

Then minimizing \(\left|\pi-(a+b\alpha)\right|\) over all admissible \(a\) is the same as minimizing

$|x_b-a|$

over integers \(a\) with \(|a|\le 10^{13}\). Therefore the best \(a\) is the nearest integer to \(x_b\), and the two-variable problem becomes

$\min_{|b|\le 10^{13}} \left\| \pi-b\alpha \right\|_{\mathbb Z},$

where \(\|y\|_{\mathbb Z}\) means the distance from \(y\) to the nearest integer.

There is also a strong bound on \(b\). The trivial candidate \(a=3\), \(b=0\) has error \(\tau=\pi-3\lt 1\), so any true optimum must have error below \(1\). If \(|b|\alpha \gt 10^{13}+\pi+1\), then for every admissible \(a\),

$\left|a+b\alpha-\pi\right|\ge |b|\alpha-|a|-\pi \gt 1,$

which cannot be optimal. The implementation therefore uses the safe bound

$|b|\le B,\qquad B=\min\!\left(10^{13},\left\lfloor\frac{10^{13}+\pi+2}{\alpha}\right\rfloor\right).$

Step 2: Use convergents of \(\sqrt d\)

Because \(\alpha=\sqrt d\) is a quadratic irrational, its continued fraction is periodic, and its convergents \(p/q\) give exceptionally good rational approximations. For each convergent we write

$q\alpha=p+\delta,$

so \(\delta=q\alpha-p\) is very small.

This matters because changing \(b\) by \(q\) changes \(b\alpha\) by almost the integer \(p\):

$ (b+q)\alpha=b\alpha+p+\delta.$

So within the arithmetic progression \(b=b_0+tq\), the fractional part of \(b\alpha\) drifts only by the tiny amount \(t\delta\). That turns a huge search over \(b\) into a very small local search around each convergent.

Step 3: Split the search into residue classes modulo \(q\)

Since every convergent satisfies \(\gcd(p,q)=1\), the inverse of \(p\) modulo \(q\) exists. For each residue \(j\in\{0,1,\dots,q-1\}\), choose \(b_0\) so that

$b_0p\equiv j \pmod q.$

Then \(b_0p=j+kq\) for some integer \(k\). If we now write

$b=b_0+tq,$

we get

$\begin{aligned} b\alpha &=b_0\alpha+tq\alpha \\ &=\frac{b_0(p+\delta)}{q}+t(p+\delta) \\ &=k+tp+\frac{j}{q}+\frac{b_0\delta}{q}+t\delta. \end{aligned}$

Since \(\pi=3+\tau\), it follows that

$\pi-b\alpha=\left(3-k-tp\right)+\left(\tau-\frac{j}{q}-\frac{b_0\delta}{q}-t\delta\right).$

The first parenthesis is an integer, so the quality of the approximation is controlled entirely by the small residual

$\operatorname{diff}(t)=\frac{j}{q}+\frac{b_0\delta}{q}-\tau+t\delta.$

We want this quantity to be as close to \(0\) as possible.

Step 4: Locate the best \(t\) from a linear estimate

For fixed \(q\), \(p\), and residue class \(j\), the expression \(\operatorname{diff}(t)\) is linear in \(t\). Therefore the best \(t\) is near the real minimizer

$t\approx -\frac{\frac{j}{q}+\frac{b_0\delta}{q}-\tau}{\delta}.$

Only residue classes with \(j/q\) near \(\tau\) can compete, because \(|b_0\delta/q|\le |\delta|\) is tiny for a useful convergent. The implementations therefore inspect only a small fixed window of \(j\)-values around \(\tau q\), and for each such class they inspect only a small fixed window of \(t\)-values around the estimate above.

The admissible \(t\)-range is also clipped by \(|b|\le B\). Its endpoints are checked explicitly so that a best solution on the boundary is never missed.

Step 5: Recover \(a\) and compare candidates safely

Once a candidate \(b\) is known, we compute

$x_b=\pi-b\alpha$

and take the nearest integer \(a_0\) to \(x_b\). Because the computation uses finite-precision arithmetic, the implementations test the three integers

$a_0-1,\qquad a_0,\qquad a_0+1$

instead of trusting one rounding step blindly. This guarantees that the true nearest admissible integer is captured even if \(x_b\) lies very close to a half-integer.

The best pair is the one with smallest absolute error. If two candidates are numerically tied, the implementations keep the one with smaller \(|a|\), and then smaller \(|b|\), to make the output deterministic.

Worked Example: \(d=2\), \(N=10\)

Here \(\alpha=\sqrt 2\approx 1.414213562\) and \(\tau=\pi-3\approx 0.141592654\).

Take the convergent

$\frac{p}{q}=\frac{3}{2},$

for which

$\delta=2\sqrt 2-3\approx -0.171572875.$

Since \(\tau q\approx 0.283\), the most relevant residue is \(j=0\). Because \(p\equiv 1\pmod 2\), the corresponding base residue is \(b_0\equiv 0\pmod 2\), so we take \(b_0=0\).

Then

$\operatorname{diff}(t)=-\tau+t\delta.$

The linear estimate gives

$t\approx -\frac{-\tau}{\delta}\approx -0.825,$

so \(t=-1\) is the natural nearby candidate. This yields

$b=b_0+tq=-2.$

Now compute

$\pi-b\sqrt 2=\pi+2\sqrt 2\approx 5.969019778,$

whose nearest integer is \(a=6\). Hence the approximation is

$6-2\sqrt 2\approx 3.171572875,$

with error about

$|\,\pi-(6-2\sqrt 2)\,|\approx 0.029980221.$

This matches the hard-coded sanity check used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical search. For each non-square \(d\lt 100\), they compute \(\sqrt d\) with high precision, derive the safe bound \(B\), and generate continued-fraction convergents until the denominator exceeds that bound.

For every convergent, the implementation forms the small error term \(\delta=q\sqrt d-p\), looks only at residue classes \(j\) near \((\pi-3)q\), converts each class into a base value \(b_0\) by modular inversion, and derives the permitted interval of \(t\)-values from the constraint \(|b|\le B\).

Instead of scanning every possible \(t\), the implementation evaluates only a short interval around the linear estimate for the best \(t\), plus the interval endpoints. Each resulting \(b\) is converted into a nearby integer coefficient \(a\), and the error is checked explicitly.

The C++ version distributes the independent \(d\)-values across worker threads. The Python version performs the same search directly with decimal arithmetic. The Java version delegates execution to the same compiled numeric core, so its output agrees with the C++ behavior.

Complexity Analysis

For fixed \(d\), the denominators of the convergents of a quadratic irrational grow exponentially, so only \(O(\log B)\) convergents are needed before the denominator passes the search bound. The inspected windows in \(j\) and \(t\) are fixed-size constants, so each convergent contributes only \(O(1)\) candidate evaluations.

Therefore the practical running time per \(d\) is \(O(\log B)\), which here is effectively \(O(\log N)\), and the memory usage is \(O(\log B)\) if the convergents are stored explicitly. Since there are only finitely many non-square values \(d\lt 100\), the total cost is this per-\(d\) work multiplied by a small constant, with straightforward parallelization in C++.

Footnotes and References

1. Problem page: Project Euler 591
2. Continued fractions: Wikipedia - Continued fraction
3. Convergents of continued fractions: Wikipedia - Convergent
4. Quadratic irrationals: Wikipedia - Quadratic irrational
5. Diophantine approximation: Wikipedia - Diophantine approximation