Problem 737: Coin Loops

Project Euler Problem 737 Solution

EulerSolve provides an optimized solution for Project Euler Problem 737, Coin Loops, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a target of $L$ complete loops, we want the smallest number of coins $N(L)$ whose cumulative turning angle exceeds $2\pi L$. The checkpoint values used by the implementations are $N(1)=31$, $N(2)=154$, and $N(10)=6947$. The key observation is that the geometry can be encoded by the harmonic numbers $$H_n=\sum_{k=1}^{n}\frac{1}{k},\qquad r_n=\sqrt{\frac{H_n}{n}},$$ so the problem becomes a summation problem for a slowly decreasing sequence of angles rather than a brute-force geometric simulation. Mathematical Approach The fast solution separates the angle contributed by the current outer boundary from the much smaller interior angle increments. That makes it possible to keep an exact formulation for small cases and a block-accelerated asymptotic formulation for the large target. Step 1: Express the Coin Chain as an Angle Sum For stage $n$, define two geometric angles by $$\beta_n=\arctan\left(\frac{\sqrt{4-r_n^2}}{r_n}\right),\qquad \gamma_n=\arctan\left(\frac{n\,r_n\sqrt{4-r_n^2}}{n r_n^2+2}\right).$$ The turning contributed by the first transition is simply $\theta_1=\beta_1$....

Detailed mathematical approach

Problem Summary

For a target of $L$ complete loops, we want the smallest number of coins $N(L)$ whose cumulative turning angle exceeds $2\pi L$. The checkpoint values used by the implementations are $N(1)=31$, $N(2)=154$, and $N(10)=6947$.

The key observation is that the geometry can be encoded by the harmonic numbers

$$H_n=\sum_{k=1}^{n}\frac{1}{k},\qquad r_n=\sqrt{\frac{H_n}{n}},$$

so the problem becomes a summation problem for a slowly decreasing sequence of angles rather than a brute-force geometric simulation.

Mathematical Approach

The fast solution separates the angle contributed by the current outer boundary from the much smaller interior angle increments. That makes it possible to keep an exact formulation for small cases and a block-accelerated asymptotic formulation for the large target.

Step 1: Express the Coin Chain as an Angle Sum

For stage $n$, define two geometric angles by

$$\beta_n=\arctan\left(\frac{\sqrt{4-r_n^2}}{r_n}\right),\qquad \gamma_n=\arctan\left(\frac{n\,r_n\sqrt{4-r_n^2}}{n r_n^2+2}\right).$$

The turning contributed by the first transition is simply $\theta_1=\beta_1$. For every later transition, the previously accumulated contact geometry subtracts $\gamma_{n-1}$, so

$$\theta_n=\beta_n-\gamma_{n-1}\qquad (n\ge 2).$$

$$T_m=\sum_{n=1}^{m}\theta_n,$$

then $m$ is the number of completed transitions, and the required number of coins is

$$N(L)=m+1.$$

Here $m$ is the first index for which $T_m>2\pi L$.

Step 2: Simplify the Local Increment

Introduce the smaller angle

$$\phi_n=\beta_n-\gamma_n.$$

Then the total turn can be rewritten as

$$T_m=\beta_m+\sum_{n=1}^{m-1}\phi_n.$$

This identity is exactly what the accelerated routine exploits: $\beta_m$ is a single boundary term, while the long sum consists only of the small increments $\phi_n$.

Using the tangent subtraction formula and $r_n^2=H_n/n$, we obtain

$$\tan(\phi_n)=\frac{\sqrt{4-r_n^2}}{r_n(2n+1)}=\frac{\sqrt{n/H_n-1/4}}{n+1/2}.$$

Therefore

$$\phi_n=\arctan\left(\frac{\sqrt{n/H_n-1/4}}{n+1/2}\right),\qquad \beta_n=\arctan\left(\sqrt{\frac{4n}{H_n}-1}\right).$$

This is the main closed form behind the implementation.

Step 3: Estimate the Harmonic Numbers for Large $n$

Directly updating $H_n$ term by term is exact, but too slow when the target loop count is large. The code therefore replaces $H_n$ by its Euler-Maclaurin expansion once $n$ is big:

$$H_n=\log n+\gamma+\frac{1}{2n}-\frac{1}{12n^2}+\frac{1}{120n^4}-\frac{1}{252n^6}+O(n^{-8}),$$

where $\gamma$ is the Euler-Mascheroni constant.

From this we see that

$$\phi_n\sim \frac{1}{\sqrt{nH_n}}\sim \frac{1}{\sqrt{n\log n}},$$

so the total turn keeps growing, but extremely slowly. That slow divergence is the reason a naive term-by-term scan becomes expensive for the full problem.

Step 4: Replace $\arctan$ by a Short Odd Polynomial

Let

$$x_n=\frac{\sqrt{n/H_n-1/4}}{n+1/2}.$$

For large $n$, this quantity is small, so

$$\arctan(x_n)=x_n-\frac{x_n^3}{3}+\frac{x_n^5}{5}+O(x_n^7).$$

The bulk summation therefore uses the polynomial

$$P(x)=x-\frac{x^3}{3}+\frac{x^5}{5}$$

instead of calling the exact inverse tangent at every large index. The implementations also begin the accelerated accumulation with a tiny fixed angular correction, compensating for the small bias introduced by truncating the power series.

Step 5: Sum Large Ranges in Blocks

Define

$$\widetilde{H}(t)=\log t+\gamma+\frac{1}{2t}-\frac{1}{12t^2}+\frac{1}{120t^4}-\frac{1}{252t^6},$$

and

$$f(t)=P\left(\frac{\sqrt{t/\widetilde{H}(t)-1/4}}{t+1/2}\right).$$

Over a symmetric block of width $2k+1$ centered at $c$, the long sum

$$\sum_{j=-k}^{k} f(c+j)$$

is approximated by a quadratic model for the interior plus Euler-Maclaurin-style endpoint corrections:

$$\sum_{j=-k}^{k} f(c+j)\approx 2k\,f(c)+\frac{k^3}{3}f''(c)+\frac{f(c-k)+f(c+k)}{2}+\frac{f'(c+k)-f'(c-k)}{12}.$$

The second derivative and the endpoint derivatives are estimated from the five sampled values $f(c-2k)$, $f(c-k)$, $f(c)$, $f(c+k)$, and $f(c+2k)$ via central differences. This lets the program jump over tens of thousands of indices at a time.

Because $\beta_m<\pi/2$ for every $m$, the block stage only needs to push the partial sum beyond

$$2\pi L-\frac{\pi}{2}.$$

After that, a short forward scan with the exact boundary term $\beta_m$ is enough to finish.

Worked Example: One Full Loop

The first few exact turns are

$$\theta_1=\frac{\pi}{3}\approx 1.04719755,\qquad \theta_2\approx 0.59936515,\qquad \theta_3\approx 0.44076494.$$

Continuing the exact recurrence gives

$$T_{29}\approx 6.24142680<2\pi,\qquad T_{30}\approx 6.33370271>2\pi.$$

So 30 transitions are not enough, but 31 coins are enough. Hence

$$N(1)=31.$$

How the Code Works

The C++, Python, and Java implementations all use the same mathematics. The exact recurrence is retained for small loop counts: it updates $H_n$, $r_n$, the carried boundary angle, and the accumulated turn until the threshold $2\pi L$ is crossed.

For the large target, the implementation first sums a fixed prefix term by term. It then switches to the polynomial model $P(x)$ together with the harmonic approximation $\widetilde{H}(n)$, advancing in wide symmetric blocks and estimating the block contribution from five sampled points. Once the reduced target $2\pi L-\pi/2$ is exceeded, it rolls back one block, rebuilds the local harmonic value, and performs a short final scan while checking the exact boundary term $\beta_n$.

The C++ version also includes explicit checkpoint assertions for the small cases $L=1$, $L=2$, and $L=10$, confirming that both branches match the same mathematical model.

Complexity Analysis

If the exact branch needs $m$ transitions, its running time is $O(m)$ and its memory usage is $O(1)$.

In the accelerated branch, let $n_0$ be the fixed prefix length and let $k$ be the half-width of a block. Then the total work is

$$O\!\left(n_0+\frac{m-n_0}{2k+1}+k\right),$$

because each wide block is processed in constant time and the tail scan is short. The memory usage remains $O(1)$. In practice this is the crucial improvement: the program scales with the number of sampled blocks, not with the full number of coins.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=737
Harmonic number: Wikipedia — Harmonic number
Euler-Maclaurin formula: Wikipedia — Euler-Maclaurin formula
Inverse trigonometric functions: Wikipedia — Inverse trigonometric functions
Finite difference: Wikipedia — Finite difference

Problem 737 source code

C++

#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>

namespace {

using i64 = std::int64_t;

constexpr long double kPi = 3.141592653589793238462643383279502884L;

long double harmonic_approx(const long double n) {
    const long double n2 = n * n;
    return std::logl(n) + 0.5772156649L + 1.0L / (2.0L * n) - 1.0L / (12.0L * n2) +
           1.0L / (120.0L * n2 * n2) - 1.0L / (252.0L * n2 * n2 * n2);
}

long double tan_phi(const long double n, const long double hn) {
    return std::sqrtl(n / hn - 0.25L) / (n + 0.5L);
}

long double atan_poly(const long double x) {
    const long double x2 = x * x;
    return x - x * x2 / 3.0L + x * x2 * x2 / 5.0L;
}

long double integral_block(const long double k, const long double f_n, const long double f_n_der2) {
    return 2.0L * k * f_n + f_n_der2 * k * k * k / 3.0L;
}

i64 coins_needed_exact(const i64 loops) {
    const long double target = 2.0L * kPi * static_cast<long double>(loops);
    long double harmonic = 1.0L;
    long double radius = 1.0L;
    long double tan_gamma = 0.0L;
    long double accumulated = 0.0L;
    i64 n = 1;

    while (accumulated <= target) {
        const long double r2 = radius * radius;
        const long double root = std::sqrtl(4.0L - r2);
        const long double tan_beta = root / radius;
        const long double tan_theta = (tan_beta - tan_gamma) / (1.0L + tan_beta * tan_gamma);
        accumulated += std::atanl(tan_theta);

        const long double nl = static_cast<long double>(n);
        tan_gamma = (nl * radius * root) / (nl * r2 + 2.0L);

        const i64 m = n + 1;
        harmonic += 1.0L / static_cast<long double>(m);
        radius = std::sqrtl(harmonic / static_cast<long double>(m));
        ++n;
    }

    return n;
}

i64 coins_needed_approx(const i64 loops) {
    long double total_angle = -0.16103076705762L / 180.0L * kPi;
    i64 n = 0;
    long double harmonic = 0.0L;
    constexpr i64 initial = 1'000'000;
    constexpr i64 k = 20'000;

    for (i64 i = 0; i < initial; ++i) {
        ++n;
        harmonic += 1.0L / static_cast<long double>(n);
        total_angle += atan_poly(tan_phi(static_cast<long double>(n), harmonic));
    }

    i64 center = n + k + 1;
    long double block_sum = 0.0L;
    const long double stage1_target = static_cast<long double>(loops) * 2.0L * kPi - kPi / 2.0L;

    while (total_angle < stage1_target) {
        const long double f_nm2k =
            atan_poly(tan_phi(static_cast<long double>(center - 2 * k), harmonic_approx(center - 2.0L * k)));
        const long double f_nmk =
            atan_poly(tan_phi(static_cast<long double>(center - k), harmonic_approx(center - 1.0L * k)));
        const long double f_n = atan_poly(tan_phi(static_cast<long double>(center), harmonic_approx(center)));
        const long double f_npk =
            atan_poly(tan_phi(static_cast<long double>(center + k), harmonic_approx(center + 1.0L * k)));
        const long double f_np2k =
            atan_poly(tan_phi(static_cast<long double>(center + 2 * k), harmonic_approx(center + 2.0L * k)));

        const long double f_n_der2 = (f_npk + f_nmk - 2.0L * f_n) / (static_cast<long double>(k) * k);
        const long double f_npk_der = (-f_n + f_np2k) / (2.0L * static_cast<long double>(k));
        const long double f_nmk_der = (f_n - f_nm2k) / (2.0L * static_cast<long double>(k));
        block_sum = integral_block(static_cast<long double>(k), f_n, f_n_der2) + (f_nmk + f_npk) / 2.0L +
                    (f_npk_der - f_nmk_der) / 12.0L;

        total_angle += block_sum;
        center += 2 * k + 1;
    }

    center -= 2 * k + 1;
    total_angle -= block_sum;
    n = center - k;
    harmonic = harmonic_approx(static_cast<long double>(n));

    long double d_theta = std::atanl(std::sqrtl(4.0L * static_cast<long double>(n) / harmonic - 1.0L));
    long double d_phi = atan_poly(tan_phi(static_cast<long double>(n), harmonic));
    const long double stage2_target = static_cast<long double>(loops) * 2.0L * kPi;

    while (total_angle + d_theta < stage2_target) {
        total_angle += d_phi;
        ++n;
        harmonic += 1.0L / static_cast<long double>(n);
        d_phi = atan_poly(tan_phi(static_cast<long double>(n), harmonic));
        d_theta = std::atanl(std::sqrtl(4.0L * static_cast<long double>(n) / harmonic - 1.0L));
    }

    return n + 1;
}

i64 coins_needed(const i64 loops) {
    if (loops <= 50) {
        return coins_needed_exact(loops);
    }
    return coins_needed_approx(loops);
}

}  // namespace

int main() {
    assert(coins_needed(1) == 31);
    assert(coins_needed(2) == 154);
    assert(coins_needed(10) == 6947);
    std::cout << coins_needed(2020) << '\n';
    return 0;
}

Python

import math

kPi = 3.141592653589793238462643383279502884

def harmonic_approx(n):
    n2 = n * n
    return math.log(n) + 0.5772156649 + 1.0 / (2.0 * n) - 1.0 / (12.0 * n2) + \
           1.0 / (120.0 * n2 * n2) - 1.0 / (252.0 * n2 * n2 * n2)

def tan_phi(n, hn):
    return math.sqrt(n / hn - 0.25) / (n + 0.5)

def atan_poly(x):
    x2 = x * x
    return x - x * x2 / 3.0 + x * x2 * x2 / 5.0

def integral_block(k, f_n, f_n_der2):
    return 2.0 * k * f_n + f_n_der2 * k * k * k / 3.0

def coins_needed_exact(loops):
    target = 2.0 * kPi * loops
    harmonic = 1.0
    radius = 1.0
    tan_gamma = 0.0
    accumulated = 0.0
    n = 1
    
    while accumulated <= target:
        r2 = radius * radius
        root = math.sqrt(4.0 - r2)
        tan_beta = root / radius
        tan_theta = (tan_beta - tan_gamma) / (1.0 + tan_beta * tan_gamma)
        accumulated += math.atan(tan_theta)
        
        nl = float(n)
        tan_gamma = (nl * radius * root) / (nl * r2 + 2.0)
        
        m = n + 1
        harmonic += 1.0 / m
        radius = math.sqrt(harmonic / m)
        n += 1
        
    return n

def coins_needed_approx(loops):
    total_angle = -0.16103076705762 / 180.0 * kPi
    n = 0
    harmonic = 0.0
    initial = 1000000
    k = 20000
    
    for i in range(initial):
        n += 1
        harmonic += 1.0 / n
        total_angle += atan_poly(tan_phi(n, harmonic))
        
    center = n + k + 1
    block_sum = 0.0
    stage1_target = loops * 2.0 * kPi - kPi / 2.0
    
    while total_angle < stage1_target:
        f_nm2k = atan_poly(tan_phi(center - 2 * k, harmonic_approx(center - 2.0 * k)))
        f_nmk = atan_poly(tan_phi(center - k, harmonic_approx(center - 1.0 * k)))
        f_n = atan_poly(tan_phi(center, harmonic_approx(center)))
        f_npk = atan_poly(tan_phi(center + k, harmonic_approx(center + 1.0 * k)))
        f_np2k = atan_poly(tan_phi(center + 2 * k, harmonic_approx(center + 2.0 * k)))
        
        f_n_der2 = (f_npk + f_nmk - 2.0 * f_n) / (k * k)
        f_npk_der = (-f_n + f_np2k) / (2.0 * k)
        f_nmk_der = (f_n - f_nm2k) / (2.0 * k)
        
        block_sum = integral_block(k, f_n, f_n_der2) + (f_nmk + f_npk) / 2.0 + (f_npk_der - f_nmk_der) / 12.0
        
        total_angle += block_sum
        center += 2 * k + 1
        
    center -= 2 * k + 1
    total_angle -= block_sum
    n = center - k
    harmonic = harmonic_approx(n)
    
    d_theta = math.atan(math.sqrt(4.0 * n / harmonic - 1.0))
    d_phi = atan_poly(tan_phi(n, harmonic))
    stage2_target = loops * 2.0 * kPi
    
    while total_angle + d_theta < stage2_target:
        total_angle += d_phi
        n += 1
        harmonic += 1.0 / n
        d_phi = atan_poly(tan_phi(n, harmonic))
        d_theta = math.atan(math.sqrt(4.0 * n / harmonic - 1.0))
        
    return n + 1

def solve():
    return str(coins_needed_approx(2020))

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

Java

public class Euler737 {
    static final double kPi = 3.141592653589793238462643383279502884;

    static double harmonicApprox(double n) {
        double n2 = n * n;
        return Math.log(n) + 0.5772156649 + 1.0 / (2.0 * n) - 1.0 / (12.0 * n2) +
                1.0 / (120.0 * n2 * n2) - 1.0 / (252.0 * n2 * n2 * n2);
    }

    static double tanPhi(double n, double hn) {
        return Math.sqrt(n / hn - 0.25) / (n + 0.5);
    }

    static double atanPoly(double x) {
        double x2 = x * x;
        return x - x * x2 / 3.0 + x * x2 * x2 / 5.0;
    }

    static double integralBlock(double k, double f_n, double f_n_der2) {
        return 2.0 * k * f_n + f_n_der2 * k * k * k / 3.0;
    }

    static long coinsNeededExact(long loops) {
        double target = 2.0 * kPi * loops;
        double harmonic = 1.0;
        double radius = 1.0;
        double tanGamma = 0.0;
        double accumulated = 0.0;
        long n = 1;

        while (accumulated <= target) {
            double r2 = radius * radius;
            double root = Math.sqrt(4.0 - r2);
            double tanBeta = root / radius;
            double tanTheta = (tanBeta - tanGamma) / (1.0 + tanBeta * tanGamma);
            accumulated += Math.atan(tanTheta);

            double nl = (double) n;
            tanGamma = (nl * radius * root) / (nl * r2 + 2.0);

            long m = n + 1;
            harmonic += 1.0 / m;
            radius = Math.sqrt(harmonic / m);
            n++;
        }

        return n;
    }

    static long coinsNeededApprox(long loops) {
        double totalAngle = -0.16103076705762 / 180.0 * kPi;
        long n = 0;
        double harmonic = 0.0;
        long initial = 1000000;
        long k = 20000;

        for (long i = 0; i < initial; ++i) {
            n++;
            harmonic += 1.0 / n;
            totalAngle += atanPoly(tanPhi(n, harmonic));
        }

        long center = n + k + 1;
        double blockSum = 0.0;
        double stage1Target = loops * 2.0 * kPi - kPi / 2.0;

        while (totalAngle < stage1Target) {
            double fnm2k = atanPoly(tanPhi(center - 2 * k, harmonicApprox(center - 2.0 * k)));
            double fnmk = atanPoly(tanPhi(center - k, harmonicApprox(center - 1.0 * k)));
            double fn = atanPoly(tanPhi(center, harmonicApprox(center)));
            double fnpk = atanPoly(tanPhi(center + k, harmonicApprox(center + 1.0 * k)));
            double fnp2k = atanPoly(tanPhi(center + 2 * k, harmonicApprox(center + 2.0 * k)));

            double fnDer2 = (fnpk + fnmk - 2.0 * fn) / (k * k);
            double fnpkDer = (-fn + fnp2k) / (2.0 * k);
            double fnmkDer = (fn - fnm2k) / (2.0 * k);

            blockSum = integralBlock(k, fn, fnDer2) + (fnmk + fnpk) / 2.0 + (fnpkDer - fnmkDer) / 12.0;

            totalAngle += blockSum;
            center += 2 * k + 1;
        }

        center -= 2 * k + 1;
        totalAngle -= blockSum;
        n = center - k;
        harmonic = harmonicApprox(n);

        double dTheta = Math.atan(Math.sqrt(4.0 * n / harmonic - 1.0));
        double dPhi = atanPoly(tanPhi(n, harmonic));
        double stage2Target = loops * 2.0 * kPi;

        while (totalAngle + dTheta < stage2Target) {
            totalAngle += dPhi;
            n++;
            harmonic += 1.0 / n;
            dPhi = atanPoly(tanPhi(n, harmonic));
            dTheta = Math.atan(Math.sqrt(4.0 * n / harmonic - 1.0));
        }

        return n + 1;
    }

    public static String solve() {
        return Long.toString(coinsNeededApprox(2020));
    }

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