Detailed mathematical approach

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?

The stopping condition is

$\left\lfloor\frac{x+\left\lceil n/x\right\rceil}{2}\right\rfloor=x.$

Because \(x+\lceil n/x\rceil\) is an integer, this is equivalent to

$2x\le x+\left\lceil\frac{n}{x}\right\rceil \le 2x+1,$

hence

$x\le \left\lceil\frac{n}{x}\right\rceil \le x+1.$

Converting this back to inequalities in \(n\) gives the exact stabilization band

$x(x-1) < n \le x(x+1).$

So all integers in

$[x(x-1)+1,\;x(x+1)]$

stop immediately when the current estimate is \(x\). The code exploits this implicitly: whenever the next iterate equals the current one, recursion stops on that subinterval.

3. Initial Guess from the Number of Digits

The initial estimate depends only on the digit length \(d\):

$x_0= \begin{cases} 2\cdot 10^{(d-1)/2},& d\text{ odd},\\ 7\cdot 10^{(d-2)/2},& d\text{ even}. \end{cases}$

This is natural because the square roots of \(d\)-digit numbers all lie in one decimal magnitude band.

If \(d=2k+1\) is odd, then

$10^k \le \sqrt{n} < \sqrt{10}\,10^k,$

so \(2\cdot 10^k\) is a reasonable central first guess.

If \(d=2k\) is even, then

$\sqrt{10}\,10^{k-1} \le \sqrt{n} < 10^k,$

and \(7\cdot 10^{k-1}\) lies comfortably inside that range.

4. Summing Iterations over an Interval

Instead of processing one number \(n\) at a time, the code processes an entire interval \([L,H]\) that currently shares the same estimate \(x\). Define

$T([L,H],x)=\sum_{n=L}^{H}\tau(n,x).$

Every number contributes at least one step for the current update, so

$T([L,H],x)=|[L,H]|+\sum_{n=L}^{H}\mathbf{1}_{y(n,x)\ne x}\,\tau(n,y(n,x)).$

If we group together all \(n\) that have the same next value \(y\), we obtain the interval recursion

$T([L,H],x)=|[L,H]|+\sum_{y\ne x}T(I_y,y),$

where \(I_y\subseteq[L,H]\) is the set of numbers whose next iterate is exactly \(y\).

5. How the Subintervals \(I_y\) Are Constructed

The quantity

$c(n,x)=\left\lceil\frac{n}{x}\right\rceil$

is piecewise constant in \(n\). For a fixed integer \(c\), the condition \(\lceil n/x\rceil=c\) is equivalent to

$ (c-1)x < n \le cx,$

so over integers the interval is

$[(c-1)x+1,\;cx].$

Within such a block the next iterate is constant, because

$y=\left\lfloor\frac{x+c}{2}\right\rfloor.$

Conversely, for a fixed next value \(y\), the floor identity

$\left\lfloor\frac{x+c}{2}\right\rfloor=y$

means that \(x+c\) can only be \(2y\) or \(2y+1\). Therefore the possible \(c\)-values are

$c=2y-x\qquad\text{or}\qquad c=2y-x+1.$

This is exactly why the code uses

$cl=2y-x,\qquad ch=cl+1.$

After intersecting these with the actual range \(c_{\min}\le c\le c_{\max}\), it converts them back to an integer subinterval in \(n\).

6. Worked Interval Example for Two-Digit Numbers

For \(d=2\), the initial guess is \(x_0=7\), and the full domain is \([10,99]\).

Here

$c_{\min}=\left\lceil\frac{10}{7}\right\rceil=2,\qquad c_{\max}=\left\lceil\frac{99}{7}\right\rceil=15.$

The interval partitions by next iterate as follows:

$[10,14]\to 4,\qquad [15,28]\to 5,\qquad [29,42]\to 6,$

$[43,56]\to 7,\qquad [57,70]\to 8,\qquad [71,84]\to 9,$

$[85,98]\to 10,\qquad [99,99]\to 11.$

The middle block \([43,56]\) is already stable because it is exactly the band

$7\cdot 6 < n \le 7\cdot 8,$

that is, \(42 < n \le 56\).

7. Small Numerical Examples

For \(n=10\) with initial guess \(x_0=7\), the sequence is

$7\to 4\to 3\to 3,$

so the iteration count is \(3\).

For \(n=24\), the sequence is

$7\to 5\to 5,$

so the count is \(2\).

For \(n=99\), the sequence is

$7\to 11\to 10\to 10,$

so the count is again \(3\).

These examples show that the iterate is not necessarily monotone in \(x\); what matters is that the interval recursion reproduces the exact same branches as the original per-number process.

8. Validation Checkpoint

The C++ file contains a checkpoint for 5-digit numbers:

$A_5=3.2102888889\ldots$

More precisely, the code asserts that the average lies strictly between

$3.2102888888 \quad\text{and}\quad 3.2102888890.$

This is a useful sanity check before computing the 14-digit case.

How the Code Works

The helper ceil_div computes \(\lceil a/b\rceil\). The recursive routine total_iterations(lo, hi, x) returns \(T([lo,hi],x)\): it first adds the interval size \((hi-lo+1)\) for the current step, then computes the possible next values \(y\), derives the corresponding \(n\)-subinterval for each \(y\), and recurses only when \(y\ne x\). The wrapper average_iterations_for_digits(d) builds the \(d\)-digit domain, chooses the correct \(x_0\), and divides the total by the number of \(d\)-digit inputs, namely

$9\cdot 10^{d-1}.$

Complexity Analysis

The naive method would inspect every one of the \(9\cdot 10^{d-1}\) integers separately. The interval method instead visits only the subintervals induced by the map \(n\mapsto y(n,x)\). Each recursive call partitions its interval into a small number of child intervals, and those children are disjoint. So the runtime is governed by the number of distinct interval/estimate pairs that actually occur, which is dramatically smaller than the size of the full input range.

Memory usage is just the recursion stack, because no large tables are stored.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=255
2. Integer square roots and Newton-like methods: Wikipedia — Integer square root
3. Floor and ceiling functions: Wikipedia — Floor and ceiling functions
4. Newton's method: Wikipedia — Newton's method
5. Interval arithmetic intuition: Wikipedia — Interval