Problem 318: 2011 Nines

Project Euler Problem 318 Solution

EulerSolve provides an optimized solution for Project Euler Problem 318, 2011 Nines, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every pair of positive integers $$1\le p \lt q,\qquad p+q\le 2011,$$ consider $$\alpha=\sqrt p+\sqrt q.$$ We look at the even powers $\alpha^{2n}$, and $C(p,q,n)$ is the number of consecutive 9s at the beginning of the fractional part of $\alpha^{2n}$. Let $N(p,q)$ be the smallest $n$ such that $$C(p,q,n)\ge 2011.$$ The goal is to compute $$\sum_{p+q\le 2011} N(p,q),$$ but only for those pairs where the fractional part of $\alpha^{2n}$ approaches $1$. Mathematical Approach 1) Introduce the conjugate partner. Set $$\beta=\sqrt q-\sqrt p.$$ Then $$\alpha^{2n}+\beta^{2n}$$ is always an integer. The reason is that in the binomial expansions of $(\sqrt q+\sqrt p)^{2n}$ and $(\sqrt q-\sqrt p)^{2n}$, all odd-radical terms cancel and only integer terms remain. So if we define $$A_n=\alpha^{2n}+\beta^{2n}\in\mathbb Z,$$ then $$\alpha^{2n}=A_n-\beta^{2n}.$$ 2) Exactly when does the fractional part approach $1$? Because $p \lt q$, we have $\beta \gt 0$. The fractional part of $\alpha^{2n}$ approaches $1$ exactly when $\beta^{2n}\to0$, i.e. $$0\lt \beta \lt 1.$$ If $\beta\ge1$, then $\beta^{2n}$ does not decay to $0$, so the fractional part cannot converge to $1$. Thus the admissible pairs are precisely those with $$\sqrt q-\sqrt p \lt 1.$$ 3) Fractional part as $1-\beta^{2n}$....

Detailed mathematical approach

Problem Summary

For positive integers $p\lt q$ put $\alpha=\sqrt p+\sqrt q$ and look at the fractional part of the even powers $\alpha^{2n}$. For $p=2$, $q=3$ the statement lists eight of them, from $(\sqrt2+\sqrt3)^2=9.8989794\ldots$ to $(\sqrt2+\sqrt3)^{16}$, whose fractional part is $0.99999998915\ldots$: they begin with more and more nines, approaching $1$. Only pairs with that property count. For such a pair, $C(p,q,n)$ is the number of consecutive nines at the beginning of the fractional part of $\alpha^{2n}$, and $N(p,q)$ is the least $n$ with $C(p,q,n)\ge 2011$. We must find ^[1]

$$\sum_{p\lt q,\ p+q\le 2011}N(p,q).$$

The whole solution in one sentence: the conjugate $\beta=\sqrt q-\sqrt p$ makes $\alpha^{2n}+\beta^{2n}$ a whole number, so for $\beta\lt1$ the fractional part of $\alpha^{2n}$ is exactly $1-\beta^{2n}$, a logarithm counts the nines, and the $41\,987$ qualifying pairs are added up one by one.

Mathematical Approach

1. The conjugate and a hidden integer

Nothing can be done with $\alpha^{2n}$ alone: for $p=2$, $q=3$ it is irrational, and we need its digits. The trick is a second number, invisible in the problem. Put

$$\alpha=\sqrt q+\sqrt p,\qquad \beta=\sqrt q-\sqrt p.$$

Expand both $2n$-th powers with the binomial theorem: $\alpha^{2n}$ adds every term $\binom{2n}{k}p^{k/2}q^{(2n-k)/2}$, while $\beta^{2n}$ adds it with the sign $(-1)^k$. The odd-$k$ terms cancel in the sum, and for $k=2j$ the powers $p^{j}q^{n-j}$ are whole numbers, so

$$\begin{aligned}A_n&=\alpha^{2n}+\beta^{2n}\\&=2\sum_{j=0}^{n}\binom{2n}{2j}p^{j}q^{n-j}\end{aligned}$$

is a positive integer for every $n\ge0$. These numbers also obey a two-term recurrence: $\alpha^2$ and $\beta^2$ add up to $2(p+q)$ and multiply to $(q-p)^2$, so they are the roots of

$$x^2-2(p+q)x+(q-p)^2=0,$$

and hence their power sums satisfy

$$A_n=2(p+q)A_{n-1}-(q-p)^2A_{n-2},$$

with $A_0=2$ and $A_1=2(p+q)$. For $p=2$, $q=3$ this is $A_n=10A_{n-1}-A_{n-2}$ and gives ^[2]

$$\begin{aligned}&2,\,10,\,98,\,970,\,9602,\\&95050,\,940898,\,9313930,\,92198402.\end{aligned}$$

Compare them with the powers in the statement, $9.898\ldots$, $97.989\ldots$, $969.998\ldots$: each integer part is one less than $A_n$.

2. The fractional part in closed form

Suppose $0\lt\beta\lt1$. Then $0\lt\beta^{2n}\lt1$ for every $n\ge1$, and since $\alpha^{2n}=A_n-\beta^{2n}$ with $A_n$ a whole number, $\alpha^{2n}$ sits just below $A_n$. Therefore

$$\begin{aligned}\lfloor\alpha^{2n}\rfloor&=A_n-1,\\\{\alpha^{2n}\}&=1-\beta^{2n},\end{aligned}$$

where $\{x\}$ is the fractional part. For $p=2$, $q=3$ we have $\beta=\sqrt3-\sqrt2=0.3178\ldots$ and $\beta^2=5-2\sqrt6=0.101020\ldots$, so $\beta^4=0.0102051\ldots$ and $\{\alpha^{4}\}=0.9897948\ldots$, which is what the statement prints. The irrational power has become a tiny explicit error term, and everything else follows.

3. Exactly which pairs qualify

The statement's condition is $\{\alpha^{2n}\}\to1$. We claim this happens precisely when $\beta\lt1$.

One direction is Section 2: if $0\lt\beta\lt1$ then $\{\alpha^{2n}\}=1-\beta^{2n}\to1$. For the converse, suppose $\beta\ge1$ and write $u=\alpha^2$, $v=\beta^2$, the two roots above, so $u\gt v\ge1$.

If $pq$ is a perfect square, $\sqrt{pq}$ and hence $v=p+q-2\sqrt{pq}$ are whole numbers, so $\alpha^{2n}=A_n-v^{n}$ is an integer for every $n$: the fractional part is $0$, never near $1$. Otherwise $\sqrt{pq}$ is irrational, ^[3] and the odd-numbered binomial terms, which cancelled in the sum, survive in the difference:

$$u^n-v^n=2\sqrt{pq}\,D_n,\qquad D_n\in\mathbb Z_{\gt0}.$$

With $u^n+v^n=A_n$ this shows that $u^n$ and $v^n$ are irrational for every $n\ge1$.

Now assume, for contradiction, that $\{\alpha^{2n}\}\to1$. Because $\alpha^{2n}=A_n-v^n$, this says $v^n$ comes ever closer to its nearest integer $c_n$; put $\eta_n=v^n-c_n\to0$. As $v$ is a root of the quadratic, $v^{n+1}=2(p+q)v^{n}-(q-p)^2v^{n-1}$, and subtracting the same combination of the integers $c_n$ gives

$$\begin{aligned}&c_{n+1}-2(p+q)c_n+(q-p)^2c_{n-1}\\&\quad=-\eta_{n+1}+2(p+q)\eta_n\\&\qquad-(q-p)^2\eta_{n-1}.\end{aligned}$$

The left side is a whole number and the right side tends to $0$, so both vanish for $n$ beyond some $n_0$. From there the errors obey the same recurrence as the powers, whose roots $u\ne v$ are distinct, so ^[4]

$$\eta_n=Cu^{n}+C'v^{n}\qquad(n\ge n_0).$$

As $u\gt v\ge1$, a nonzero $C$ would make $|\eta_n|$ grow without bound and a nonzero $C'$ would keep it at least $|C'|$; since $\eta_n\to0$, both vanish. Then $v^n=c_n$ for large $n$, contradicting the irrationality just proved, so $\beta\ge1$ is impossible. In other words $v$ would have to be a Pisot number, which it is not, because its conjugate $u$ exceeds $1$. ^[5]

The test needs no square roots. As $q\ge p+1$, the inequality $\sqrt q\lt1+\sqrt p$ is equivalent to $q\lt p+1+2\sqrt p$, that is to

$$(q-p-1)^2\lt 4p.$$

Equality means $\beta=1$, which happens for consecutive squares $(p,q)=(a^2,(a+1)^2)$: $(1,4)$, $(4,9)$, $(9,16)$, and so on. The borderline is tight. For $q-p-1=61$, compare $61^2=3721$ with $4p$: $(930,992)$ has $4p=3720$ and is excluded, while $(931,993)$ has $4p=3724$ and is kept, with $\beta=0.99960\ldots$ just under the limit.

4. Counting the nines

A fractional part $x\in(0,1)$ begins with at least $K$ nines exactly when $x\ge1-10^{-K}$, the number $0.99\ldots9$ with $K$ nines. By Section 2, $x=1-\beta^{2n}$, so

$$C(p,q,n)\ge K\iff \beta^{2n}\le10^{-K}.$$

Take logarithms to base $10$ and write

$$\lambda=-\log_{10}\beta^{2}=-2\log_{10}\beta\gt0.$$

Then $\beta^{2n}\le10^{-K}$ becomes $K\le n\lambda$. As this holds for every $K$ up to $n\lambda$, the count itself is

$$C(p,q,n)=\lfloor n\lambda\rfloor .$$

The count is thus a staircase of slope $\lambda$: automatically non-decreasing, the pattern the statement points out. For $p=2$, $q=3$, $\lambda=0.9955904\ldots$, so $\lfloor0.99559n\rfloor$ equals $0,1,2,3,4,5,6,7$ for $n=1,\dots,8$ — exactly the eight lines of the statement. The slope is just below $1$, so the apparent rule "one nine more per step" must eventually slip, and it does: the first repetition is $C(2,3,226)=C(2,3,227)=225$, as $1/(1-\lambda)=226.779185\ldots$ predicts.

5. The minimal exponent

With $K=2011$, the least $n$ with $C(p,q,n)\ge K$ is the least $n$ with $n\lambda\ge K$, so

$$N(p,q)=\lfloor K/\lambda\rfloor+1 .$$

The ceiling $\lceil K/\lambda\rceil$ gives the same value, because $K/\lambda$ is never a whole number: $K=n\lambda$ would make $\beta^{2n}=10^{-K}$ rational, while $\beta^{2n}$ is irrational for every qualifying pair. (Section 3 needs $pq$ not to be a perfect square, automatic here: $p=da^2$ and $q=db^2$ would force $\beta=(b-a)\sqrt d\ge1$.) The two forms agree in exact arithmetic; in floating point either fails if the quotient rounds across an integer, which Section 8 examines.

6. A complete small example

Take the bound $p+q\le5$ and ask for one nine, $K=1$. The pairs with $p\lt q$ are $(1,2)$, $(1,3)$, $(1,4)$ and $(2,3)$. The pair $(1,4)$ fails the test, $(4-1-1)^2=4=4p$, so $\beta=2-1=1$ and its powers are integers. For the others:

$(1,2)$: $\beta^2=3-2\sqrt2=0.1715728\ldots$, $\lambda=0.7655513\ldots$, $1/\lambda=1.3062\ldots$, so $N=2$;
$(1,3)$: $\beta^2=4-2\sqrt3=0.5358983\ldots$, $\lambda=0.2709175\ldots$, $1/\lambda=3.6911\ldots$, so $N=4$;
$(2,3)$: $\beta^2=5-2\sqrt6=0.1010205\ldots$, $\lambda=0.9955904\ldots$, $1/\lambda=1.0044\ldots$, so $N=2$.

The sum is $2+4+2=8$, the second checkpoint of the C++ program. Its first reads the statement through the same formula: $1/\lambda$, $2/\lambda$, $3/\lambda$ for $(2,3)$ are $1.0044\ldots$, $2.0088\ldots$, $3.013287\ldots$, so the first power with one, two or three nines is $n=2,3,4$, matching $97.98\ldots$, $969.998\ldots$, $9601.9998\ldots$. At $K=2011$ that pair gives $K/\lambda=2019.9069\ldots$ and $N(2,3)=2020$.

7. The shape of the sum

For fixed $p$, the condition $(q-p-1)^2\lt4p$ says $q-p-1\lt2\sqrt p$, that is

$$p\lt q\le p+\lceil 2\sqrt p\rceil .$$

The ceiling is exact in both cases: for a non-square $p$ the number $2\sqrt p$ is irrational, and for $p=a^2$ it equals $2a$, the last admissible value, since $q=(a+1)^2$ gives equality. So $p=1$ admits only $q=2,3$; $p=100$ the $20$ values $q=101,\dots,120$; $p=101$ one more, up to $q=122$. For most $p$ only about $2\sqrt p$ of the $2011-2p$ permitted values of $q$ survive, so qualifying pairs are scarce: $41\,987$ of the $1005^2=1\,010\,025$ pairs, about one in $24$.

Their terms are also wildly unequal. A pair whose $\beta$ is far below $1$ has a large $\lambda$ and a small $N$: $(1005,1006)$ has $\beta=0.01576\ldots$, $\lambda=3.6044\ldots$, $N=558$. A pair whose $\beta$ nearly reaches $1$ has a minute $\lambda$ and an enormous $N$: the largest term is $N(931,993)=5\,933\,985$, with $\lambda=0.000338895\ldots$. Multiplying the number of pairs by that largest term bounds the total by $41\,987\cdot5\,933\,985\lt2.5\cdot10^{11}$, which is what makes a $64$-bit accumulator the safe choice.

8. How accurately must $\lambda$ be known?

Everything rests on one floor function: $\lfloor K/\lambda\rfloor$ is wrong whenever the computed $K/\lambda$ lands on the other side of an integer. The first danger, the admissibility test, is removed entirely by the integer form $(q-p-1)^2\lt4p$; floating-point $\beta$ would also be safe here, since over all $1\,010\,025$ pairs the smallest nonzero $|\beta-1|$ is $0.00013\ldots$, at $(930,992)$.

The second danger is the floor itself. Here $K/\lambda$ reaches almost $6\cdot10^{6}$, so tiny relative errors matter, and they are largest where it hurts: subtracting two square roots leaves $\beta$ with an absolute error of up to $5\cdot10^{-15}$, and when $\beta$ is close to $1$ the value $\lambda=-2\log_{10}\beta$ is tiny, so that error becomes a large relative error in $\lambda$ and in the quotient. Measured in binary64 here against a $50$-digit computation, the largest absolute deviation of $K/\lambda$ is $2.4\cdot10^{-5}$, at $(813,871)$ — a figure that depends on the floating-point library. That is not negligible: the closest any $K/\lambda$ comes to an integer is $2.7\cdot10^{-5}$, at $(625,635)$, where $K/\lambda=1435.000027089\ldots$ and $N=1436$, not $1435$. But the two extremes fall on different pairs, and pair by pair the margin wins: the tightest ratio of distance-to-integer against rounding error is about $4\,000$, at $(871,931)$. Double precision is adequate here by a verified margin, not by a proof: Baker's theorem bounds nonzero linear forms in logarithms such as $K\ln10+2n\ln\beta$ away from zero, ^[6] but no such bound is evaluated here.

How the Code Works

The C++ program ^[7] transcribes Sections 3 to 5 directly. A struct Options holds max_sum and target_nines, both $2011$ by default, and a flag for the self-checks; parse_arguments and parse_int_after_prefix recognise --max-sum=, --target-nines= and --skip-checkpoints, reject unknown arguments and non-digit values, and require $\mathtt{max\_sum}\ge3$, $\mathtt{target\_nines}\ge1$ (their decimal accumulation into a signed int is unchecked, so an absurd value is not cleanly rejected). Those adjustable parameters make the derivation testable.

minimal_n(p, q, target_nines) computes $\beta$ as sqrt(q) - sqrt(p), then beta2 and lambda = -log10(beta2), the $\lambda$ of Section 4. It starts from the truncation n = (u64)(target / lambda), lifts $n$ to $1$ if that is zero, then walks: the first loop raises $n$ while $n\lambda$ is below the target, the second lowers it while $n-1$ would already suffice. In exact arithmetic that returns the least $n$ with $n\lambda\ge K$, namely $\lfloor K/\lambda\rfloor+1$, and the guess is within a step of it; the loops compare products formed from the same computed $\lambda$, so they repair the truncation, not the logarithm.

solve(max_sum, target_nines) runs the double loop over $p$ and over $q=p+1$ while $p+q\le\mathtt{max\_sum}$, recomputes $\beta$, skips the pair when beta >= 1.0, and otherwise adds minimal_n to a uint64_t total. run_checkpoints verifies $N(2,3)=2,3,4$ for one, two and three nines and that the small sum of Section 6 equals $8$, returning false on failure; main parses, runs it unless --skip-checkpoints was given, exits with code $2$ if it fails, and otherwise prints the sum.

Python and Java carry only the two computing functions: Python minimal_n and solve(max_sum=2011, target_nines=2011), whose inner bound range(p + 1, max_sum - p + 1) is the same set of $q$; Java minimalN and a parameterless solve() with both bounds fixed at $2011$. Neither parses arguments nor runs the checkpoints, and both return the total as a string. Their float and double are IEEE binary64, ^[8] while the precision of C++ long double is implementation-defined — here the same binary64, wider on some platforms; extra width certifies nothing on its own, so all three were compared against a high-precision computation. The accumulator is $64$-bit in C++ and Java, justified by the bound of Section 7, and unbounded in Python.

Two details deserve a warning rather than imitation. At the boundary $\beta=1$ the entries are consecutive squares, whose square roots are exact in binary64, so beta >= 1.0 decides those pairs exactly; away from it the test is right only because $|\beta-1|$ never drops below $0.00013$. And both comparisons in minimal_n carry an absolute slack of 1e-15: at $\mathtt{target}=2011$ that is far below the binary64 spacing there, about $2.3\cdot10^{-13}$, so in Python and Java target - 1e-15 rounds back to target and does nothing. A slack $\varepsilon$ there moves the quotient threshold by only $\varepsilon/\lambda$, at most $3\cdot10^{-12}$ here, and removing it leaves all $41\,987$ exponents unchanged.

Complexity and Verification

The double loop visits every pair with $p\lt q$ and $p+q\le M$: for each $p$ with $2p\lt M$ there are $M-2p$ values of $q$, hence $\lfloor (M-1)^2/4\rfloor$ pairs in all, which is $1\,010\,025$ for $M=2011$. Every candidate costs two square roots for the filter, and each of the $41\,987$ survivors two more, a logarithm, a division and a bounded walk: $O(M^2)$ time, $O(1)$ memory. The scan can be cut to the range of Section 7, $q\le\min(M-p,\ p+\lceil2\sqrt p\rceil)$, which visits exactly the qualifying pairs and costs $O(M^{3/2})$; the ceiling matters, since $p+\lfloor2\sqrt p\rfloor+1$ would also admit the $31$ boundary pairs whose $p$ is a square. That shortened scan gives the same total. On the machine used here all three programs finish in a fraction of a second.

The C++ program checks itself first, unless --skip-checkpoints is given: $N(2,3)$ must be $2$, $3$, $4$ for one, two and three nines, and the sum for $p+q\le5$ with one nine must be $8$; Python and Java carry no such checks. For this article the following were verified by computation.

The eight powers of $\sqrt2+\sqrt3$ in the statement, to $80$ digits, with their nine counts $0,1,\dots,7$ and the integers $A_n$; the binomial formulas for $A_n$, $D_n$ and the recurrence agree for all $p\le8$, $q\le14$, $n\le8$.
The identity $C(p,q,n)=\lfloor n\lambda\rfloor$ by exact integers. Writing $\beta^{2n}=x_n-y_n\sqrt{pq}$, with $x_0=1$, $y_0=0$ and the simultaneous update $x\mapsto(p+q)x+2pqy$, $y\mapsto2x+(p+q)y$, the test $\beta^{2n}\le10^{-K}$ reads $10^Kx_n-1\le10^Ky_n\sqrt{pq}$; both sides are non-negative, so it is equivalent to $(10^Kx_n-1)^2\le10^{2K}y_n^2pq$. The staircase holds for all $1\,025$ qualifying triples with $p\le9$, $q\le19$, $n\le25$.
$N(p,q)=\lfloor K/\lambda\rfloor+1$ against that test for all $434$ qualifying triples with $p\le12$, $q\le25$, $K\le7$, and the $397$ with $N\le45$ against the printed digits of $\alpha^{2N}$ and $\alpha^{2N-2}$.
The criterion $(q-p-1)^2\lt4p$ against high-precision $\sqrt q-\sqrt p\lt1$ for all pairs with $q\le4000$, and against the programs' float test for all $1\,010\,025$ target pairs: no disagreement.
Every $N(p,q)$ of the target recomputed at $50$ digits: all $41\,987$ match the double-precision values, and the total is the same at $120$ digits. Section 8's margins come from this comparison.
Sums for smaller parameters, by formula and by the exact test: $8$, $14$, $20$ for one, two, three nines at $p+q\le5$; $40$ and $75$ at $p+q\le10$; $149$ and $671$ for one and five nines at $p+q\le20$. At the target, C++, Python, Java and the shortened scan all print the value of the $50$-digit computation.

Footnotes and References

The sources below supply the background facts; the closed form of the fractional part, the characterisation of the qualifying pairs, the staircase and the precision margins are derived above.

Project Euler 318 — 2011 Nines. The official statement lists eight even powers of the sum of the square roots of 2 and 3, notes that such a fractional part tends to 1, defines C(p, q, n) as the number of leading nines of that fractional part and N(p, q) as the least n with C at least 2011, and asks for the sum of N(p, q) over the positive integer pairs p less than q with p + q at most 2011 whose fractional parts tend to 1.
OEIS A087799. Defined by a(n) = 10a(n−1) − a(n−2) with a(0) = 2 and a(1) = 10, terms 2, 10, 98, 970, 9602, 95050, 940898, 9313930, closed form the sum of the n-th powers of 5 + sqrt(24) and 5 − sqrt(24) — the squares of the sum and difference of the square roots of 2 and 3, so this is the sequence of Section 1 for (2, 3).
Wikipedia — Quadratic irrational number. The article states that the square root of a natural number which is not a perfect square is irrational; Section 3 uses this for the square root of pq.
Wikipedia — Linear recurrence with constant coefficients. When the characteristic roots are distinct, every solution of a homogeneous linear recurrence is a fixed combination of their n-th powers; Section 3 applies this to the rounding errors.
Wikipedia — Pisot–Vijayaraghavan number. A Pisot number is a real algebraic integer greater than 1 whose other conjugates are smaller than 1 in absolute value; the distance from its powers to the nearest integer then tends to 0 exponentially. The article also records the converse for algebraic numbers, the dichotomy proved in Section 3.
Wikipedia — Baker's theorem. A nonzero linear form in logarithms of algebraic numbers with algebraic coefficients is bounded away from zero by an effectively computable negative power of the coefficients' height; Section 8 appeals only to the existence of such a bound.
C++ — Euler318.cpp. The linked immutable revision contains parse_arguments, minimal_n, solve and run_checkpoints as described above, and is the implementation record for this article. The high-precision comparisons and the shortened scan are extra validation, not packaged tests.
Wikipedia — Double-precision floating-point format. The format has a 53-bit significand and a maximum relative rounding error of 2 to the power −53, about 1.11 times 10 to the power −16, the figures used in Section 8.

Problem 318 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;

struct Options {
    int max_sum = 2011;
    int target_nines = 2011;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--max-sum=", options.max_sum) ||
            parse_int_after_prefix(arg, "--target-nines=", options.target_nines)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.max_sum >= 3 && options.target_nines >= 1;
}

u64 minimal_n(const int p, const int q, const int target_nines) {
    const long double beta = std::sqrt(static_cast<long double>(q)) -
                             std::sqrt(static_cast<long double>(p));
    const long double beta2 = beta * beta;
    const long double lambda = -std::log10(beta2);
    const long double target = static_cast<long double>(target_nines);

    u64 n = static_cast<u64>(target / lambda);
    if (n == 0ULL) {
        n = 1ULL;
    }
    while (static_cast<long double>(n) * lambda < target - 1e-15L) {
        ++n;
    }
    while (n > 1ULL && static_cast<long double>(n - 1ULL) * lambda >= target - 1e-15L) {
        --n;
    }
    return n;
}

u64 solve(const int max_sum, const int target_nines) {
    u64 total = 0ULL;
    for (int p = 1; p < max_sum; ++p) {
        for (int q = p + 1; p + q <= max_sum; ++q) {
            const long double beta = std::sqrt(static_cast<long double>(q)) -
                                     std::sqrt(static_cast<long double>(p));
            if (beta >= 1.0L) {
                continue;
            }
            total += minimal_n(p, q, target_nines);
        }
    }
    return total;
}

bool run_checkpoints() {
    if (minimal_n(2, 3, 1) != 2ULL || minimal_n(2, 3, 2) != 3ULL || minimal_n(2, 3, 3) != 4ULL) {
        std::cerr << "Checkpoint failed for sqrt(2)+sqrt(3) nines progression" << '\n';
        return false;
    }
    if (solve(5, 1) != 8ULL) {
        std::cerr << "Checkpoint failed for small max-sum=5 target=1" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.max_sum, options.target_nines) << '\n';
    return 0;
}

Python

import math

def minimal_n(p, q, target_nines):
    beta = math.sqrt(q) - math.sqrt(p)
    beta2 = beta * beta
    lam = -math.log10(beta2)
    target = float(target_nines)
    
    n = int(target / lam)
    if n == 0:
        n = 1
    while n * lam < target - 1e-15:
        n += 1
    while n > 1 and (n - 1) * lam >= target - 1e-15:
        n -= 1
    return n

def solve(max_sum=2011, target_nines=2011):
    total = 0
    for p in range(1, max_sum):
        for q in range(p + 1, max_sum - p + 1):
            beta = math.sqrt(q) - math.sqrt(p)
            if beta >= 1.0:
                continue
            total += minimal_n(p, q, target_nines)
    return str(total)

if __name__ == '__main__':
    print(solve())

Java

public class Euler318 {
    static long minimalN(int p, int q, int targetNines) {
        double beta = Math.sqrt(q) - Math.sqrt(p);
        double beta2 = beta * beta;
        double lambda = -Math.log10(beta2);
        double target = (double) targetNines;

        long n = (long) (target / lambda);
        if (n == 0) {
            n = 1;
        }
        while (n * lambda < target - 1e-15) {
            n++;
        }
        while (n > 1 && (n - 1) * lambda >= target - 1e-15) {
            n--;
        }
        return n;
    }

    public static String solve() {
        int maxSum = 2011;
        int targetNines = 2011;
        long total = 0;

        for (int p = 1; p < maxSum; ++p) {
            for (int q = p + 1; p + q <= maxSum; ++q) {
                double beta = Math.sqrt(q) - Math.sqrt(p);
                if (beta >= 1.0) {
                    continue;
                }
                total += minimalN(p, q, targetNines);
            }
        }
        return String.valueOf(total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}