Problem 901: Well Drilling

View on Project Euler

Project Euler Problem 901 Solution

EulerSolve provides an optimized solution for Project Euler Problem 901, Well Drilling, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The depth sequence starts from $$d_0=0,\qquad d_1=x,$$ and then follows the nonlinear recurrence $$d_{k+1}=e^{d_k-d_{k-1}}\qquad (k\ge 1).$$ An admissible drilling schedule must be strictly increasing, so the required condition is $$d_0<d_1<d_2<\cdots.$$ The quantity to evaluate is $$C(x)=\sum_{k\ge 1} d_k e^{-d_{k-1}}.$$ The C++, Python, and Java implementations all treat the answer as the boundary point of the upper admissible branch: they locate the smallest \(x\) in that branch and then evaluate \(C(x)\) there. Numerically, this boundary is near \(x_\star\approx 0.746542014027\), and the corresponding value is \(C(x_\star)\approx 2.364497769\). Mathematical Approach Two identities make the search practical: one turns the monotonicity condition into a statement about successive gaps, and the other rewrites the objective as a rapidly convergent tail sum....

Detailed mathematical approach

Problem Summary

The depth sequence starts from

$$d_0=0,\qquad d_1=x,$$

and then follows the nonlinear recurrence

$$d_{k+1}=e^{d_k-d_{k-1}}\qquad (k\ge 1).$$

An admissible drilling schedule must be strictly increasing, so the required condition is

$$d_0<d_1<d_2<\cdots.$$

The quantity to evaluate is

$$C(x)=\sum_{k\ge 1} d_k e^{-d_{k-1}}.$$

The C++, Python, and Java implementations all treat the answer as the boundary point of the upper admissible branch: they locate the smallest \(x\) in that branch and then evaluate \(C(x)\) there. Numerically, this boundary is near \(x_\star\approx 0.746542014027\), and the corresponding value is \(C(x_\star)\approx 2.364497769\).

Mathematical Approach

Two identities make the search practical: one turns the monotonicity condition into a statement about successive gaps, and the other rewrites the objective as a rapidly convergent tail sum.

Step 1: Express Admissibility Through Gap Variables

Define the forward gaps

$$\Delta_k=d_k-d_{k-1}\qquad (k\ge 1).$$

The schedule is strictly increasing exactly when

$$\Delta_k>0\qquad \text{for all }k\ge 1.$$

From the recurrence we immediately get

$$d_{k+1}=e^{\Delta_k}.$$

For \(k\ge 2\), since \(d_k=e^{\Delta_{k-1}}\), the next gap satisfies

$$\Delta_{k+1}=d_{k+1}-d_k=e^{\Delta_k}-e^{\Delta_{k-1}}.$$

Because the exponential function is strictly increasing, this means

$$\Delta_{k+1}>0 \iff \Delta_k>\Delta_{k-1}\qquad (k\ge 2).$$

So admissibility is stronger than mere positivity: after the first two terms, the gap sequence must keep increasing as well.

Step 2: Rewrite the Objective in a Simpler Form

For \(k\ge 2\), substitute the recurrence into one summand of \(C(x)\):

$$d_k e^{-d_{k-1}}=e^{d_{k-1}-d_{k-2}}e^{-d_{k-1}}=e^{-d_{k-2}}.$$

Therefore the cost becomes

$$C(x)=x+\sum_{k\ge 2} e^{-d_{k-2}}=x+\sum_{j\ge 0} e^{-d_j}.$$

Expanding the first few terms gives

$$C(x)=x+1+e^{-x}+e^{-e^x}+e^{-d_3}+e^{-d_4}+\cdots.$$

This identity is very useful: once the depths become moderately large, the tail shrinks extremely fast.

Step 3: Why a Threshold Search Works

Every finite prefix \(d_0,d_1,\dots,d_n\) depends continuously on the starting value \(x\). If one sampled value of \(x\) eventually produces a non-positive gap while a nearby value keeps all tested gaps positive, then a boundary point must lie between them.

The implementations exploit this by searching for the left edge of the upper admissible branch. Below that edge, some future gap becomes non-positive and strict increase fails. Above that edge, the simulated trajectory stays increasing and the depths quickly escape upward.

This is why bisection is appropriate: once a failing value and a succeeding value are bracketed, repeated midpoint testing isolates the smallest admissible \(x\) on that branch.

Step 4: Why the Boundary Value Is the Relevant Candidate

On the upper admissible branch, any smaller starting value is forbidden because it eventually breaks monotonicity. That makes the left boundary \(x_\star\) the natural candidate explored by the implementations.

The C++ implementation also checks a separated low admissible branch numerically and confirms that its sampled costs remain above the value found at the upper boundary. In other words, the final reported value comes from the smallest admissible \(x\) on the upper branch, not from a larger interior point.

Step 5: Truncating the Infinite Sum Safely

Once some depth exceeds \(80\), every later tail term in

$$x+\sum_{j\ge 0} e^{-d_j}$$

is at most \(e^{-80}\approx 1.8\times 10^{-35}\). That is astronomically smaller than the \(10^{-9}\) precision of the printed answer. So the implementations can stop the simulation after the depth crosses \(80\), with a generous hard cap of \(200\) recurrence steps as a safety limit.

Worked Example

Take \(x=0.75\), which lies very close to the upper-branch boundary used by the implementations. The first depths are

$$\begin{aligned} d_0&=0,\\ d_1&=0.75,\\ d_2&=e^{0.75}\approx 2.117000017,\\ d_3&=e^{2.117000017-0.75}\approx 3.923562400,\\ d_4&=e^{3.923562400-2.117000017}\approx 6.089478119,\\ d_5&=e^{6.089478119-3.923562400}\approx 8.722585700. \end{aligned}$$

The cost identity now gives

$$\begin{aligned} C(0.75)&=0.75+1+e^{-0.75}+e^{-2.117000017}+e^{-3.923562400}+e^{-6.089478119}+\cdots\\ &\approx 0.75+1+0.472366553+0.120392262+0.019770539+0.002266591+\cdots\\ &\approx 2.3648. \end{aligned}$$

This already sits very close to the boundary value reported by the implementations, and it shows how quickly the tail becomes negligible.

How the Code Works

The C++, Python, and Java implementations all simulate the recurrence directly from \(d_0=0\) and a trial value of \(d_1=x\). During this simulation they check the gap \(d_k-d_{k-1}\) at each step, and they reject the trial immediately if any gap becomes non-positive.

To locate the upper admissible branch, the implementations first perform a coarse scan on the interval \(0.6\le x\le 1.5\) with step \(0.01\). This finds the first transition from failure to success on that branch.

After bracketing the transition, they run \(90\) bisection iterations. The midpoint is tested by the same recurrence simulation, so the boundary estimate tightens until the desired floating-point precision is reached.

Once the boundary value has been isolated, the implementations replay the recurrence and accumulate the cost. They stop when the depth exceeds \(80\), because the remaining tail is far below the displayed precision.

The C++ implementation adds extra sanity checks around the boundary and compares against sampled points on the low branch, while the Python and Java implementations keep only the core search-and-evaluate procedure.

Complexity Analysis

Let \(T\) be the maximum number of recurrence steps used in a single simulation. Here \(T\le 200\). One admissibility test or one cost evaluation therefore costs \(O(T)\) time and \(O(1)\) memory.

If the coarse scan uses \(S\) samples and the refinement uses \(B\) bisection steps, the total running time is

$$O((S+B)T).$$

Since \(B=O(\log(1/\varepsilon))\) for target precision \(\varepsilon\), this is also \(O(T\log(1/\varepsilon))\) once the initial scan interval is fixed. Because all constants are small and fixed in the implementations, the practical runtime is effectively constant and the memory usage remains \(O(1)\).

Footnotes and References

  1. Problem page: Project Euler 901
  2. Recurrence relation: Wikipedia — Recurrence relation
  3. Exponential function: Wikipedia — Exponential function
  4. Bisection method: Wikipedia — Bisection method
  5. Fixed-point iteration: Wikipedia — Fixed-point iteration

Problem 901 source code

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>

using namespace std;

static constexpr int MAX_STEPS = 200;
static constexpr long double STOP_DEPTH = 80.0L;

static bool is_increasing(long double d1, long double *min_diff_out = nullptr) {
    long double prev = 0.0L;
    long double curr = d1;
    long double min_diff = numeric_limits<long double>::infinity();
    for (int i = 0; i < MAX_STEPS; ++i) {
        long double diff = curr - prev;
        if (diff <= 0.0L) {
            if (min_diff_out) {
                *min_diff_out = diff;
            }
            return false;
        }
        if (diff < min_diff) {
            min_diff = diff;
        }
        if (curr > STOP_DEPTH) {
            break;
        }
        long double next = expl(curr - prev);
        prev = curr;
        curr = next;
    }
    if (min_diff_out) {
        *min_diff_out = min_diff;
    }
    return true;
}

static long double expected_cost(long double d1) {
    long double prev = 0.0L;
    long double curr = d1;
    long double sum = 0.0L;
    for (int i = 0; i < MAX_STEPS; ++i) {
        sum += curr * expl(-prev);
        if (curr > STOP_DEPTH) {
            break;
        }
        long double next = expl(curr - prev);
        prev = curr;
        curr = next;
    }
    return sum;
}

static long double find_high_threshold() {
    const long double scan_start = 0.6L;
    const long double scan_end = 1.5L;
    const long double scan_step = 0.01L;

    bool prev_inc = is_increasing(scan_start);
    long double low = scan_start;
    long double high = scan_start;
    for (long double d = scan_start + scan_step; d <= scan_end; d += scan_step) {
        bool inc = is_increasing(d);
        if (!prev_inc && inc) {
            low = d - scan_step;
            high = d;
            break;
        }
        prev_inc = inc;
    }
    if (high == scan_start) {
        return -1.0L;
    }

    for (int i = 0; i < 90; ++i) {
        long double mid = (low + high) * 0.5L;
        if (is_increasing(mid)) {
            high = mid;
        } else {
            low = mid;
        }
    }
    return high;
}

static bool run_validation(long double d1, long double answer) {
    bool ok = true;
    if (!(d1 > 0.7L && d1 < 0.8L)) {
        cerr << "Validation failed: threshold out of expected bracket.\n";
        ok = false;
    }
    if (!is_increasing(d1)) {
        cerr << "Validation failed: threshold not increasing.\n";
        ok = false;
    }
    long double eps = 1e-6L;
    if (is_increasing(d1 - eps)) {
        cerr << "Validation failed: threshold not minimal.\n";
        ok = false;
    }
    long double answer_up = expected_cost(d1 + eps);
    if (!(answer_up >= answer - 1e-12L)) {
        cerr << "Validation failed: expected cost not monotone near threshold.\n";
        ok = false;
    }

    long double best_low = numeric_limits<long double>::infinity();
    for (long double x = 0.02L; x <= 0.3L; x += 0.001L) {
        if (!is_increasing(x)) {
            continue;
        }
        long double v = expected_cost(x);
        if (v < best_low) {
            best_low = v;
        }
    }
    if (!(best_low > answer + 1e-4L)) {
        cerr << "Validation failed: low-branch minimum not higher than optimum.\n";
        ok = false;
    }
    if (!(answer > 2.35L && answer < 2.38L)) {
        cerr << "Validation failed: answer outside sanity range.\n";
        ok = false;
    }
    if (ok) {
        cerr << "Validation checkpoints passed.\n";
    }
    return ok;
}

int main() {
    long double d1 = find_high_threshold();
    if (d1 <= 0.0L) {
        cerr << "Failed to locate threshold.\n";
        return 1;
    }
    long double answer = expected_cost(d1);

    if (!run_validation(d1, answer)) {
        return 1;
    }

    cout.setf(ios::fixed);
    cout << setprecision(9) << answer << "\n";
    return 0;
}

Python

import math

MAX_STEPS = 200
STOP_DEPTH = 80.0

def is_increasing(d1):
    prev = 0.0
    curr = d1
    min_diff = float('inf')
    for _ in range(MAX_STEPS):
        diff = curr - prev
        if diff <= 0.0:
            return False, diff
        if diff < min_diff:
            min_diff = diff
        if curr > STOP_DEPTH:
            break
        next_val = math.exp(curr - prev)
        prev = curr
        curr = next_val
    return True, min_diff

def is_increasing_bool(d1):
    return is_increasing(d1)[0]

def expected_cost(d1):
    prev = 0.0
    curr = d1
    sum_cost = 0.0
    for _ in range(MAX_STEPS):
        sum_cost += curr * math.exp(-prev)
        if curr > STOP_DEPTH:
            break
        next_val = math.exp(curr - prev)
        prev = curr
        curr = next_val
    return sum_cost

def find_high_threshold():
    scan_start = 0.6
    scan_end = 1.5
    scan_step = 0.01

    prev_inc = is_increasing_bool(scan_start)
    low = scan_start
    high = scan_start
    
    d = scan_start + scan_step
    while d <= scan_end + 1e-9:
        inc = is_increasing_bool(d)
        if not prev_inc and inc:
            low = d - scan_step
            high = d
            break
        prev_inc = inc
        d += scan_step

    if high == scan_start:
        return -1.0

    for _ in range(90):
        mid = (low + high) * 0.5
        if is_increasing_bool(mid):
            high = mid
        else:
            low = mid
    return high

def solve():
    d1 = find_high_threshold()
    if d1 <= 0.0:
        return "Failed"
    answer = expected_cost(d1)
    return f"{answer:.9f}"

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

Java

public class Euler901 {
    static final int MAX_STEPS = 200;
    static final double STOP_DEPTH = 80.0;

    static boolean isIncreasing(double d1) {
        double prev = 0.0;
        double curr = d1;
        double minDiff = Double.POSITIVE_INFINITY;
        for (int i = 0; i < MAX_STEPS; ++i) {
            double diff = curr - prev;
            if (diff <= 0.0) {
                return false;
            }
            if (diff < minDiff) {
                minDiff = diff;
            }
            if (curr > STOP_DEPTH) {
                break;
            }
            double nextVal = Math.exp(curr - prev);
            prev = curr;
            curr = nextVal;
        }
        return true;
    }

    static double expectedCost(double d1) {
        double prev = 0.0;
        double curr = d1;
        double sum = 0.0;
        for (int i = 0; i < MAX_STEPS; ++i) {
            sum += curr * Math.exp(-prev);
            if (curr > STOP_DEPTH) {
                break;
            }
            double nextVal = Math.exp(curr - prev);
            prev = curr;
            curr = nextVal;
        }
        return sum;
    }

    static double findHighThreshold() {
        double scanStart = 0.6;
        double scanEnd = 1.5;
        double scanStep = 0.01;

        boolean prevInc = isIncreasing(scanStart);
        double low = scanStart;
        double high = scanStart;

        for (double d = scanStart + scanStep; d <= scanEnd + 1e-9; d += scanStep) {
            boolean inc = isIncreasing(d);
            if (!prevInc && inc) {
                low = d - scanStep;
                high = d;
                break;
            }
            prevInc = inc;
        }

        if (high == scanStart) {
            return -1.0;
        }

        for (int i = 0; i < 90; ++i) {
            double mid = (low + high) * 0.5;
            if (isIncreasing(mid)) {
                high = mid;
            } else {
                low = mid;
            }
        }
        return high;
    }

    public static String solve() {
        double d1 = findHighThreshold();
        if (d1 <= 0.0) {
            return "Failed";
        }
        double answer = expectedCost(d1);
        return String.format(java.util.Locale.US, "%.9f", answer);
    }

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