Problem 255: Rounded Square Roots

Project Euler Problem 255 Solution

EulerSolve provides an optimized solution for Project Euler Problem 255, Rounded Square Roots, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every $d$-digit integer $n$, the problem applies a rounded square-root iteration starting from a digit-based initial guess $x_0$. We must count how many iteration steps are needed until the estimate stops changing, then average that count over all $d$-digit numbers. If we write the current estimate as $x$, the next estimate is defined by $$c(n,x)=\left\lceil\frac{n}{x}\right\rceil,\qquad y(n,x)=\left\lfloor\frac{x+c(n,x)}{2}\right\rfloor.$$ The process stops when $y(n,x)=x$. The brute-force approach would repeat this for every $n$ separately, but for $d=14$ there are $9\cdot 10^{13}$ inputs, so the code instead groups numbers into intervals that share the same next iterate. Mathematical Approach 1. The Rounded Newton Step The map $$x\longmapsto \left\lfloor\frac{x+\left\lceil n/x\right\rceil}{2}\right\rfloor$$ is an integer version of Heron's method / Newton's method for $\sqrt{n}$. The real Newton step would use $n/x$; here that quotient is first rounded upward by the ceiling function and then averaged with the current guess. For a fixed $n$, define the iteration count recursively by $$\tau(n,x)= \begin{cases} 1,& y(n,x)=x,\\ 1+\tau(n,y(n,x)),& y(n,x)\ne x. \end{cases}$$ The leading $1$ means that the current application of the update rule always counts as one step. 2. When Does the Iteration Stabilize?...

Detailed mathematical approach

Problem Summary

For a positive integer $n$ with $d$ decimal digits, the procedure of the problem starts from $x_0=2\cdot10^{(d-1)/2}$ when $d$ is odd and from $x_0=7\cdot10^{(d-2)/2}$ when $d$ is even, and repeats

$$x_{k+1}=\left\lfloor\frac{x_k+\lceil n/x_k\rceil}{2}\right\rfloor$$

until $x_{k+1}=x_k$. It stops at the rounded square root of $n$, the integer nearest to $\sqrt n$. Every use of the formula counts as one iteration, the last one included: for $n=4321$ the values are $x_0=70$, $x_1=66$ and $x_2=66$, two iterations, and indeed $\sqrt{4321}=65.734\ldots$. The $90\,000$ five-digit numbers need $3.2102888889$ iterations on average, rounded to $10$ decimal places. We must find the average number of iterations over all fourteen-digit numbers, $10^{13}\le n\lt10^{14}$, rounded to $10$ decimal places. ^[1] The procedure is an integer form of Heron's method, which improves a guess $x$ for $\sqrt n$ by averaging $x$ and $n/x$. ^[2]

The whole solution in one sentence: for a fixed current value $x$ the next value depends on $n$ only through the rounded-up quotient $\lceil n/x\rceil$, so the numbers that have taken the same values so far form blocks of consecutive integers, and we count the steps of whole blocks instead of single numbers.

Mathematical Approach

We first find out where the procedure stops and why it needs so few steps. Then we group the numbers into blocks that take the same steps and count the iterations block by block.

1. One step of the procedure

For a fixed $n$ we write one step as

$$s_n(x)=\left\lfloor\frac{x+\lceil n/x\rceil}{2}\right\rfloor,$$

and $\tau(n)$ for the number of iterations, the smallest $k\ge1$ with $x_k=x_{k-1}$. In the example of the statement, $4321/70=61.73\ldots$ is rounded up to $62$ and $\lfloor(70+62)/2\rfloor=66$; then $4321/66=65.47\ldots$ is rounded up to $66$, and $\lfloor(66+66)/2\rfloor=66$ ends the run with $\tau(4321)=2$. Three two-digit numbers, which all start from $x_0=7$:

$n=10$: $7\to4\to3\to3$, three iterations;
$n=24$: $7\to5\to5$, two iterations;
$n=99$: $7\to11\to10\to10$, three iterations; the first step goes up, because $7$ is below $\sqrt{99}\approx9.95$.

Without the rounding, the step $x\mapsto(x+n/x)/2$ is Newton's method for the equation $x^2-n=0$. ^[2] The usual integer square root iteration uses $\lfloor n/x\rfloor$ instead and aims at $\lfloor\sqrt n\rfloor$; it can alternate between two neighbouring integers, for $n=8$ between $3$ and $2$, and needs a different stopping test. ^[3] With the ceiling, as the next two sections show, every run settles on a single value.

2. Where the procedure stops

Claim: $s_n(x)=x$ holds exactly when $x(x-1)\lt n\le x(x+1)$, and this happens for exactly one $x$, the integer $r$ nearest to $\sqrt n$.

Proof. Put $c=\lceil n/x\rceil$. By the definition of the floor, $s_n(x)=x$ means $x\le(x+c)/2\lt x+1$, that is $x\le c\le x+1$, since $c$ is an integer. By the definition of the ceiling, $c\ge x$ means $n/x\gt x-1$, and $c\le x+1$ means $n/x\le x+1$; together, $x^2-x\lt n\le x^2+x$. ^[4] Now $(x-\tfrac12)^2=x^2-x+\tfrac14$ and $(x+\tfrac12)^2=x^2+x+\tfrac14$, and no integer lies strictly between $x^2\pm x$ and $x^2\pm x+\tfrac14$. For an integer $n$ the condition is therefore the same as $(x-\tfrac12)^2\lt n\lt(x+\tfrac12)^2$, that is $x-\tfrac12\lt\sqrt n\lt x+\tfrac12$, which says that $x$ is the integer nearest to $\sqrt n$. There is never a tie, because $(x+\tfrac12)^2$ is not an integer.

So the procedure stops at $r$ for exactly the $2r$ numbers from $r^2-r+1$ to $r^2+r$, which we call the stable band of $r$. The band of $66$ runs from $4291$ to $4422$ and contains $4321$; the band of $7$ runs from $43$ to $56$.

3. Every run slides down to the rounded root

Two facts describe every run. First, $s_n(x)\ge r$ for every $x\ge1$. The inequality of arithmetic and geometric means gives $x+n/x\ge2\sqrt n$, ^[5] and $2\sqrt n\gt2r-1$ because $\sqrt n\gt r-\tfrac12$. So the integer $x+\lceil n/x\rceil$ is greater than $2r-1$, hence at least $2r$, and halving and rounding down gives at least $r$. Second, $s_n(x)\lt x$ whenever $x\gt r$. Indeed $x\ge r+1$ gives $x(x-1)\ge(r+1)r\ge n$ by Section 2, so $n/x\le x-1$, $\lceil n/x\rceil\le x-1$ and $s_n(x)\le\lfloor(2x-1)/2\rfloor=x-1$.

Hence $x_1\ge r$ whatever $x_0$ is, and from then on the values decrease strictly until they reach $r$, where they stay. The procedure always stops, and always at the rounded square root. If $x_0=r$, then $\tau(n)=1$; otherwise $\tau(n)=j+1$, where $j\ge1$ is the first index with $x_j=r$, because one more iteration is needed to see that $r$ does not change. The run of $99$ shows both facts: it rises from $7$ to $11$, above $r=10$, and then slides down to $10$.

4. Why so few steps

For the exact Heron step a little algebra gives

$$\frac{x+n/x}{2}-\sqrt n=\frac{(x-\sqrt n)^2}{2x}:$$

the error is squared and divided by $2x$, so once it is small compared with $\sqrt n$, the number of correct digits roughly doubles with every step. ^[2] The integer step adds one rounding error: $x+\lceil n/x\rceil$ lies between $x+n/x$ and $x+n/x+1$, the upper end excluded, and rounding its half down lowers that half by at most $\tfrac12$. With $e=x-\sqrt n$ the new error therefore satisfies

$$\frac{e^2}{2x}-\frac12\le s_n(x)-\sqrt n\lt\frac{e^2}{2x}+\frac12.$$

The square root of a fourteen-digit number lies between $\sqrt{10^{13}}\approx3162278$ and $10^7$, so the ratio of $x_0=7\cdot10^6$ to the true root is greater than $0.7$ and at most about $2.21$. For the smallest such number, $n=10^{13}$ with $\sqrt n\approx3162277.66$, the run is

$$\begin{aligned}&7000000,\ 4214286,\ 3293584,\\&3164895,\ 3162279,\\&3162278,\ 3162278,\end{aligned}$$

with errors of about $3837722$, $1052008$, $131306$, $2617$, $1.34$ and $0.34$, six iterations in all; for instance $2617^2/(2\cdot3164895)\approx1.08$, and the next value lies $1.34$ above $\sqrt n$. The largest count among the fourteen-digit numbers is $7$. An example is $n=10000005311561$, which equals $r(r+1)-1$ for $r=3162278$: its square root lies just below $r+\tfrac12$, so the rounding error matters until the end, and the run finishes with $3162280$, $3162279$, $r$, $r$.

5. The first step for all two-digit numbers at once

Instead of running the numbers one by one, we follow all $90$ two-digit numbers through their first step from $x_0=7$. As $n$ runs from $10$ to $99$, the rounded-up quotient $\lceil n/7\rceil$ runs from $2$ to $15$, and the next value $\lfloor(7+\lceil n/7\rceil)/2\rfloor$ runs from $4$ to $11$:

$10\ldots14$ go to $4$, $15\ldots28$ to $5$ and $29\ldots42$ to $6$;
$43\ldots56$ go to $7$; this is the stable band of $7$, so these $14$ numbers are finished after one iteration;
$57\ldots70$ go to $8$, $71\ldots84$ to $9$, $85\ldots98$ to $10$, and $99$ alone goes to $11$.

Apart from the two ends, every group has exactly $14=2\cdot7$ numbers.

6. Numbers with the same next value form a block

Fix the current value $x$. By the definition of the ceiling, $\lceil n/x\rceil=c$ holds exactly for $(c-1)x\lt n\le cx$, a block of $x$ consecutive integers. ^[4] The next value $\lfloor(x+c)/2\rfloor$ equals $y$ exactly when $x+c$ is $2y$ or $2y+1$, that is, for the two ceilings $c=2y-x$ and $c=2y-x+1$. Joining their blocks, the positive integers $n$ whose next value is $y$ are those with

$$(2y-x-1)\,x\lt n\le(2y-x+1)\,x,$$

a range of $2x$ consecutive integers, except that for even $x$ and the smallest next value $y=x/2$ only its upper half is positive. For $x=7$ and $y=5$ this gives $14\lt n\le28$, the group $15\ldots28$ above. For $y=x$ it gives $x(x-1)\lt n\le x(x+1)$, the stable band of Section 2, as it must. The groups for $y,y+1,y+2,\ldots$ follow one another without gaps, so cutting any block of numbers that share the value $x$ along these groups leaves blocks again. Starting from the block of all $d$-digit numbers, which share $x_0$, it follows step by step that the numbers with the same values $x_0,x_1,\ldots,x_k$ always form a block of consecutive integers.

7. Counting iterations block by block

For a block $L\le n\le H$ of numbers that all have the current value $x$, let $T(L,H,x)$ be the total number of iterations that these numbers still need, counted from the value $x$. Each of them performs the next step, and those whose next value $y$ differs from $x$ continue from $y$. Grouping them by $y$ as in Section 6 gives

$$\begin{aligned}&T(L,H,x)=(H-L+1)\\&\qquad+\sum_{y\ne x}T(L_y,H_y,y),\end{aligned}$$

where $L_y\le n\le H_y$ is the part of the block with next value $y$; the group $y=x$ adds nothing more, because its numbers stop. The recursion ends, because every number reaches its stable band after finitely many steps (Section 3). A number with $\tau$ iterations takes part in $\tau$ nested calls, so the nesting depth equals the largest iteration count. The requested average is

$$\frac{T(10^{13},\,10^{14}-1,\,7\cdot10^6)}{9\cdot10^{13}}.$$

8. Worked example: all two-digit numbers

The first step of Section 5 contributes $90$, one iteration for every number. Each group with a next value other than $7$ is then cut again by the rule of Section 6:

$57\ldots70$ at $8$ lies inside the stable band $57\ldots72$ of $8$ and contributes $14$;
$10\ldots14$ at $4$: $13$ and $14$ stay at $4$, while $10$, $11$ and $12$ move to $3$ and stop there, so $5+3=8$;
$15\ldots28$ at $5$: $15\ldots20$ move to $4$ and stop, $21\ldots28$ stay, so $14+6=20$;
$29\ldots42$ at $6$: only $29$ and $30$ move, to $5$, and stop, so $14+2=16$;
$71\ldots84$ at $9$: only $71$ and $72$ move, to $8$, so $14+2=16$;
$85\ldots98$ at $10$: $85\ldots90$ move to $9$, so $14+6=20$;
$99$ at $11$ moves to $10$ and stops there, so $1+1=2$.

The total is $90+14+8+20$${}+16+16+20+2=186$, and two-digit numbers need $186/90=31/15\approx2.0667$ iterations on average. Counting number by number agrees: $14$ numbers need one iteration, $56$ need two and $20$ need three, and $14+2\cdot56+3\cdot20=186$. For the five-digit numbers the same recursion, started from $x_0=200$, visits $925$ blocks and finds the total $288\,926$, and $288\,926/90\,000=3.21028888\ldots$ is the value of the statement.

9. How many blocks there are

At the top level every group except the first and the last has exactly $2x_0$ numbers, so the $9\cdot10^{13}$ fourteen-digit numbers fall into about $9\cdot10^{13}/(1.4\cdot10^7)\approx6.4$ million groups, exactly $6\,428\,572$, one of them the stable band of $x_0$. Below the top level a block never splits into more than three groups. Such a block is one group of its parent, so it has at most $2x$ numbers, where $x$ is the value of the parent, and its own value $y=\lfloor(x+c)/2\rfloor$ is at least $x/2$, because $c\ge1$. Across the block $n/y$ therefore changes by less than $2x/y\le4$, the ceiling $\lceil n/y\rceil$ takes at most five consecutive values, and since every next value collects two consecutive ceilings (Section 6), at most three next values occur. With at most seven nested calls, the number of calls is thus at most a fixed multiple of the number of top-level groups. Most blocks are in fact not cut at all: the program makes $57\,941\,869$ calls, about nine for every top-level group.

10. Rounding the average

The programs divide the total $T$ by $9\cdot10^{13}$ in floating point, so we must make sure that rounding to $10$ decimals cannot go wrong. No number needs more than $7$ iterations, so $T\le7\cdot9\cdot10^{13}\lt2^{53}$, and $T$ and $9\cdot10^{13}$ are integers that a double holds exactly. ^[6] In the default round-to-nearest mode an IEEE division returns the double nearest to the exact quotient, ^[7] and below $8$ neighbouring doubles are at most $2^{-50}\approx8.9\cdot10^{-16}$ apart, ^[6] so all three programs obtain the average with an error below $8.9\cdot10^{-16}$, even the C++ program, which rounds first to long double and then to double. The value rounded to ten decimals changes only at the halfway points $(2m+1)/(2\cdot10^{10})$, and

$$\begin{aligned}&\frac{T}{9\cdot10^{13}}-\frac{2m+1}{2\cdot10^{10}}\\&\qquad=\frac{T-4500\,(2m+1)}{9\cdot10^{13}}.\end{aligned}$$

So the exact average is either exactly halfway, when $T$ is an odd multiple of $4500$, or at least $1/(9\cdot10^{13})\approx1.1\cdot10^{-14}$ away from every halfway point, more than ten times the floating-point error. A separate check shows that $T$ is not an odd multiple of $4500$, so the computed doubles round to the same ten decimals as the exact average.

How the Code Works

The three programs implement Section 7 directly. The C++ helper ceil_div(a, b) computes $\lceil a/b\rceil$ as $\lfloor(a+b-1)/b\rfloor$. ^[4] The recursive function total_iterations(lo, hi, x) returns the total $T$ of Section 7 for the block from lo to hi: $0$ for an empty block, and otherwise the block length hi - lo + 1 for the current step plus the recursive totals of the groups. It computes the ceilings cmin and cmax of lo/x and hi/x and from them the smallest and largest next values ymin and ymax. For each y in between, cl = 2y - x and ch = cl + 1 are the two ceilings of Section 6 that lead to y; clipped to the range from cmin to cmax they become from and to, and the block from nl = (from - 1)·x + 1 to nr = to·x, clipped to the range from lo to hi, is the group of y. The group y == x is skipped, because its numbers are finished; every other group gets a recursive call. The test from > to never succeeds, because every y between ymin and ymax is reached by a ceiling in the range.

average_iterations_for_digits(d) builds the block from $10^{d-1}$ to $10^d-1$ and the starting value of the statement, calls total_iterations once and divides by the count $9\cdot10^{d-1}$. main first asserts that the five-digit average lies strictly between $3.2102888888$ and $3.2102888890$, the checkpoint of the statement, and then prints the fourteen-digit average with std::fixed and std::setprecision(10). The block bounds are unsigned 64-bit integers, ample since to·x stays below $10^{14}+x$; the difference 2y - x is signed, and the sum is an unsigned __int128, a GCC extension that Clang also accepts, ^[8] although every total stays below $2^{53}$. The division is done in long double, which with GCC on x86 typically has the 80-bit extended format and on 64-bit ARM under macOS the same format as double; ^[9] by Section 10 this does not affect the printed digits. The assert disappears if the program is compiled with NDEBUG.

The Python and Java versions follow the C++ program line by line; from becomes fr in Python, where it is a keyword, and fromVal in Java. Python nests the helpers inside solve and uses // on its unbounded integers throughout the recursion; only the final average uses /, which yields a float, ^[10] and it is formatted with f'{ans:.10f}'. Java uses long throughout and String.format(Locale.US, "%.10f", ans), so the decimal separator is a point in every locale. Neither repeats the five-digit check. The recursion is at most seven calls deep, far below any stack limit.

Complexity and Verification

A call costs a fixed number of operations plus one pass of its loop per group, so the time is proportional to the calls plus the passes: $57\,941\,869$ and $87\,275\,460$ for fourteen digits. The top-level groups, about $9\cdot10^{d-1}/(2x_0)$ of them, grow like the square root of the count of $d$-digit numbers; below them the calls branch at most threefold per level (Section 9), and since the error shrinks quadratically, the number of levels grows only very slowly with $d$; here it is seven. From $12$ to $14$ digits the range grows a hundredfold, but the calls only about elevenfold, from $5\,177\,248$ to $57\,941\,869$. Memory is the recursion stack of at most seven frames. In our runs on a desktop computer the C++ and Java programs finished in under a second, the Python program in about $40$ seconds.

The C++ program checks the five-digit average of the statement first; this checkpoint uses the starting rule for an odd number of digits, while the main computation uses the rule for an even number. For this article the following checks were run as well; they are not part of the published solutions:

A direct run of the procedure for every $d$-digit number with $d\le11$, $9\cdot10^{10}$ numbers for $d=11$ alone, gives exactly the totals of the published code.
For every $d\le14$ the published total_iterations agrees exactly with a separately written recursion that follows each ceiling $c$ on its own instead of the pairs of Section 6, and for every $d\le12$ the C++, Python and Java code give the same totals.
Direct runs of $76\,000\,003$ fourteen-digit numbers in $304$ windows, namely the first and the last $10^7$ numbers, the stable band of $x_0$ with $5\cdot10^6$ numbers on each side, the $2\,000\,002$ numbers around $10000005311561$ and $300$ random windows of $10^5$ numbers, agree window by window with total_iterations called on the window. All numbers in the blocks at the deepest level of the recursion need exactly $7$ iterations.
The claims of Sections 2 and 3 and the error bound of Section 4, checked in exact integer arithmetic, hold for every $n\le10^6$ and every start value $x\le3000$.
The C++ program, built with GCC and with Clang, and the Python and Java programs print the same ten decimals; the fourteen-digit total is below $2^{53}$ and not an odd multiple of $4500$, as Section 10 requires.

Footnotes and References

The references below supply the background facts used above. The stable bands, the shape of every run, the blocks of numbers with a common next value, the recursion and the rounding argument are derived in this article.

Project Euler 255 — Rounded Square Roots. The official statement defines the procedure with its starting values for odd and even digit counts, works the example n = 4321, gives the average 3.2102888889 for five-digit numbers and asks for the average number of iterations over all 14-digit numbers, rounded to 10 decimal places.
Wikipedia — Methods of computing square roots. Heron's method replaces an estimate x of √S by (x + S/x)/2 and is Newton's method for x² − S = 0. With x_n = √S(1 + ε_n) the relative errors satisfy ε_{n+1} = ε_n²/(2(1 + ε_n)), so the convergence is quadratic and the number of correct digits roughly doubles with each step.
Wikipedia — Integer square root. The integer form of Newton's method, x_{k+1} = ⌊(x_k + ⌊n/x_k⌋)/2⌋, computes ⌊√n⌋; its sequence eventually becomes constant or alternates between two consecutive integers, and the starting value decides how many steps are needed.
Wikipedia — Floor and ceiling functions. ⌊x⌋ = m exactly when m ≤ x < m + 1, ⌈x⌉ = n exactly when n − 1 < x ≤ n, and ⌈n/m⌉ = ⌊(n + m − 1)/m⌋ for a positive integer m.
Wikipedia — AM–GM inequality. For non-negative x and y, (x + y)/2 ≥ √(xy), with equality exactly when x = y; the proof expands (x − y)² ≥ 0.
Wikipedia — Double-precision floating-point format. A double has a 53-bit significand, all integers from −2^53 to 2^53 are exactly representable, and between 2^n and 2^(n+1) neighbouring doubles are 2^(n−52) apart.
D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys, 1991. In the IEEE standard each basic operation, division included, is computed exactly and then rounded; by default rounding means round to nearest, and three directed rounding modes are provided as well.
GCC manual — 128-bit Integers. GCC supports __int128 and unsigned __int128 on targets with an integer mode wide enough to hold 128 bits.
Wikipedia — long double. With GCC on x86, long double is the 80-bit extended format; with Microsoft Visual C++ and on 64-bit ARM under macOS and Windows it has the same format as double, which the article calls a synonym for double.
Python language reference — Binary arithmetic operations. Division of integers with / yields a float, while floor division with // yields an integer.

Problem 255 source code

C++

#include <cassert>
#include <cstdint>
#include <iomanip>
#include <iostream>

using u64 = std::uint64_t;
using i64 = std::int64_t;
using u128 = unsigned __int128;

static inline u64 ceil_div(u64 a, u64 b) {
    return (a + b - 1) / b;
}

static u128 total_iterations(u64 lo, u64 hi, u64 x) {
    if (lo > hi) return 0;

    u128 total = static_cast<u128>(hi - lo + 1);  // Current step always counts.

    const u64 cmin = ceil_div(lo, x);
    const u64 cmax = ceil_div(hi, x);
    const u64 ymin = (x + cmin) / 2;
    const u64 ymax = (x + cmax) / 2;

    for (u64 y = ymin; y <= ymax; ++y) {
        i64 cl = static_cast<i64>(2 * y) - static_cast<i64>(x);
        i64 ch = cl + 1;
        i64 from = std::max<i64>(cl, static_cast<i64>(cmin));
        i64 to = std::min<i64>(ch, static_cast<i64>(cmax));
        if (from > to) continue;

        u64 nl = static_cast<u64>(from - 1) * x + 1;
        u64 nr = static_cast<u64>(to) * x;
        if (nl < lo) nl = lo;
        if (nr > hi) nr = hi;

        if (y == x) continue;
        total += total_iterations(nl, nr, y);
    }
    return total;
}

static long double average_iterations_for_digits(int d) {
    u64 lo = 1;
    for (int i = 1; i < d; ++i) lo *= 10;
    u64 hi = lo * 10 - 1;

    u64 p10 = 1;
    int exp = (d % 2 == 1) ? ((d - 1) / 2) : ((d - 2) / 2);
    for (int i = 0; i < exp; ++i) p10 *= 10;
    u64 x0 = (d % 2 == 1) ? (2 * p10) : (7 * p10);
    u128 tot = total_iterations(lo, hi, x0);
    u64 cnt = hi - lo + 1;
    return static_cast<long double>(tot) / static_cast<long double>(cnt);
}

int main() {
    long double check = average_iterations_for_digits(5);
    assert(check > 3.2102888888L && check < 3.2102888890L);

    long double ans = average_iterations_for_digits(14);
    std::cout << std::fixed << std::setprecision(10) << static_cast<double>(ans) << '\n';
    return 0;
}

Python

def solve():
    def ceil_div(a, b):
        return (a + b - 1) // b

    def total_iterations(lo, hi, x):
        if lo > hi:
            return 0
        total = hi - lo + 1
        cmin = ceil_div(lo, x)
        cmax = ceil_div(hi, x)
        ymin = (x + cmin) // 2
        ymax = (x + cmax) // 2
        for y in range(ymin, ymax + 1):
            cl = 2 * y - x
            ch = cl + 1
            fr = max(cl, cmin)
            to = min(ch, cmax)
            if fr > to:
                continue
            nl = (fr - 1) * x + 1
            nr = to * x
            nl = max(nl, lo)
            nr = min(nr, hi)
            if y == x:
                continue
            total += total_iterations(nl, nr, y)
        return total

    def average_iterations_for_digits(d):
        lo = 10 ** (d - 1)
        hi = 10 ** d - 1
        exp = (d - 1) // 2 if d % 2 == 1 else (d - 2) // 2
        p10 = 10 ** exp
        x0 = 2 * p10 if d % 2 == 1 else 7 * p10
        tot = total_iterations(lo, hi, x0)
        cnt = hi - lo + 1
        return tot / cnt

    ans = average_iterations_for_digits(14)
    return f'{ans:.10f}'

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

Java

import java.util.Locale;

public class Euler255 {
    static long ceilDiv(long a, long b) {
        return (a + b - 1) / b;
    }

    static long totalIterations(long lo, long hi, long x) {
        if (lo > hi)
            return 0;

        long total = hi - lo + 1;

        long cmin = ceilDiv(lo, x);
        long cmax = ceilDiv(hi, x);
        long ymin = (x + cmin) / 2;
        long ymax = (x + cmax) / 2;

        for (long y = ymin; y <= ymax; ++y) {
            long cl = 2 * y - x;
            long ch = cl + 1;

            long fromVal = Math.max(cl, cmin);
            long toVal = Math.min(ch, cmax);
            if (fromVal > toVal)
                continue;

            long nl = (fromVal - 1) * x + 1;
            long nr = toVal * x;
            if (nl < lo)
                nl = lo;
            if (nr > hi)
                nr = hi;

            if (y == x)
                continue;
            total += totalIterations(nl, nr, y);
        }

        return total;
    }

    static double averageIterationsForDigits(int d) {
        long lo = 1;
        for (int i = 1; i < d; ++i)
            lo *= 10;
        long hi = lo * 10 - 1;

        long p10 = 1;
        int exp = (d % 2 == 1) ? ((d - 1) / 2) : ((d - 2) / 2);
        for (int i = 0; i < exp; ++i)
            p10 *= 10;

        long x0 = (d % 2 == 1) ? (2 * p10) : (7 * p10);
        long tot = totalIterations(lo, hi, x0);
        long cnt = hi - lo + 1;

        return (double) tot / (double) cnt;
    }

    public static String solve() {
        double ans = averageIterationsForDigits(14);
        return String.format(Locale.US, "%.10f", ans);
    }

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