Problem 144: Investigating Multiple Reflections of a Laser Beam

Project Euler Problem 144 Solution

EulerSolve provides an optimized solution for Project Euler Problem 144, Investigating Multiple Reflections of a Laser Beam, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The mirror is the ellipse $4x^2+y^2=100$, except for a very small opening at the top where $|x|\le 0.01$. The beam starts at $(0,10.1)$, enters through that opening, and its first actual impact on the mirror is $(1.4,-9.6)$. The task is to count how many reflections occur before the beam reaches the opening again and escapes. This is a deterministic billiard problem inside an ellipse. Once two consecutive points on the path are known, the next segment is forced by the law of specular reflection, so the mathematics is about deriving a stable update rule for one bounce. Mathematical Approach Let $F(x,y)=4x^2+y^2-100$. Every reflection point lies on the curve $F(x,y)=0$. The key is to express one bounce entirely in terms of the previous point and the current point. The Right State Variables A single point on the ellipse is not enough, because the next bounce depends on the incoming direction. If the last two points are $P_{k-1}$ and $P_k$, then the incoming direction at $P_k$ is $$v_{\mathrm{in}}=P_k-P_{k-1}.$$ So the true state of the system is the ordered pair $(P_{k-1},P_k)$. The implementation stores exactly that pair and updates it after every reflection. The Normal Vector and the Reflection Law For an implicit curve $F(x,y)=0$, the gradient is normal to the curve....

Detailed mathematical approach

Problem Summary

The mirror is the ellipse $4x^2+y^2=100$, except for a very small opening at the top where $|x|\le 0.01$. The beam starts at $(0,10.1)$, enters through that opening, and its first actual impact on the mirror is $(1.4,-9.6)$. The task is to count how many reflections occur before the beam reaches the opening again and escapes.

This is a deterministic billiard problem inside an ellipse. Once two consecutive points on the path are known, the next segment is forced by the law of specular reflection, so the mathematics is about deriving a stable update rule for one bounce.

Mathematical Approach

Let $F(x,y)=4x^2+y^2-100$. Every reflection point lies on the curve $F(x,y)=0$. The key is to express one bounce entirely in terms of the previous point and the current point.

The Right State Variables

A single point on the ellipse is not enough, because the next bounce depends on the incoming direction. If the last two points are $P_{k-1}$ and $P_k$, then the incoming direction at $P_k$ is

$$v_{\mathrm{in}}=P_k-P_{k-1}.$$

So the true state of the system is the ordered pair $(P_{k-1},P_k)$. The implementation stores exactly that pair and updates it after every reflection.

The Normal Vector and the Reflection Law

For an implicit curve $F(x,y)=0$, the gradient is normal to the curve. Here

$$\nabla F(x,y)=(8x,2y).$$

Therefore the normal at $P_k=(x_k,y_k)$ is $n=(8x_k,2y_k)$. If $v_{\mathrm{in}}$ is the incoming direction, the reflected direction is

$$v_{\mathrm{out}}=v_{\mathrm{in}}-2\frac{v_{\mathrm{in}}\cdot n}{n\cdot n}\,n.$$

Equivalently, one may normalize $n$ first and write the usual unit-normal formula. Both forms describe the same line. Geometrically, the tangential component is preserved and the normal component is reversed, which is exactly the law of reflection.

The Next Intersection Without a Full Quadratic Solve

Starting from the current impact point $P_k=(x_k,y_k)$, the outgoing ray is

$$P_k+t\,v_{\mathrm{out}}=(x_k+t d_x,\ y_k+t d_y).$$

Substituting this into the ellipse equation gives

$$4(x_k+t d_x)^2+(y_k+t d_y)^2=100.$$

Because $P_k$ already lies on the ellipse, the constant term vanishes, so the quadratic factorizes as

$$t\Big((4d_x^2+d_y^2)t+(8x_k d_x+2y_k d_y)\Big)=0.$$

The root $t=0$ is the current point. The other root is the next impact point:

$$t=-\frac{8x_k d_x+2y_k d_y}{4d_x^2+d_y^2}.$$

This simplification is the reason the implementations do not need a general quadratic formula. They only compute the nonzero root and step directly to the next reflection.

Worked Example: The First Reflected Segment

The beam arrives from $P_{-1}=(0,10.1)$ to $P_0=(1.4,-9.6)$, so

$$v_{\mathrm{in}}=(1.4,-19.7).$$

At $P_0$, the normal is

$$n=(8\cdot 1.4,\ 2\cdot (-9.6))=(11.2,-19.2).$$

Applying the reflection formula gives the outgoing direction

$$v_{\mathrm{out}}\approx(-16.4590673575,\ 10.9155440415).$$

Then

$$t\approx 0.3275153751,$$

and the next point is

$$P_1=P_0+t\,v_{\mathrm{out}}\approx(-3.9905976194,\ -6.0249914989).$$

This concrete first step is a useful checkpoint because all three implementations reproduce it.

Detecting the Exit

The mirror is missing a tiny arc near the top. On the ellipse, any point with $|x|\le 0.01$ and $y>0$ lies on that upper opening rather than on reflective material. So the simulation stops as soon as the next would-be impact satisfies

$$|x|\le 0.01,\qquad y>0.$$

The last counted reflection is the bounce that sends the beam toward that gap. For the given initial ray, the total number of reflections is $354$.

How the Code Works

Initialization

The C++, Python, and Java implementations store the previous point $(0,10.1)$, the current point $(1.4,-9.6)$, and a counter initialized to zero. That setup means the first loop iteration computes the reflection at the first interior impact point.

One Bounce Per Loop Iteration

Each iteration performs the same four geometric steps: build the incoming vector from the two stored points, compute the normal at the current impact point, reflect the direction, and use the nonzero root above to jump to the next point on the ellipse. After that, the state is shifted forward by replacing the previous point with the current one and the current point with the newly found point.

The C++ implementation uses the unnormalized form $v-2(v\cdot n)n/(n\cdot n)$, while the Python and Java implementations normalize the normal first. These are algebraically equivalent. After each update the counter is incremented, and the new point is tested against the opening condition. When the test succeeds, the loop stops and prints $354$.

Complexity Analysis

If the beam reflects $B$ times before leaving, the running time is $O(B)$, because every iteration uses only a fixed number of arithmetic operations. Memory usage is $O(1)$: the algorithm keeps only a few floating-point coordinates, one direction vector, and the reflection counter.

For this specific problem $B=354$, so the computation is tiny. The important point is not asymptotic difficulty but deriving a numerically stable geometric update that follows the exact reflection law.

Footnotes and References

Problem page: Project Euler 144
Ellipse: Wikipedia - Ellipse
Gradient and implicit curves: Wikipedia - Gradient
Specular reflection: Wikipedia - Specular reflection
Analytic geometry: Wikipedia - Analytic geometry

Problem 144 source code

C++

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

namespace {

struct Options {
    bool run_checkpoints = true;
};

bool nearly_equal(const double a, const double b, const double eps = 1e-9) {
    return std::abs(a - b) <= eps;
}

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;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

struct Point {
    double x = 0.0;
    double y = 0.0;
};

Point reflect_direction(const Point& incoming, const Point& point) {
    const double nx = 8.0 * point.x;
    const double ny = 2.0 * point.y;

    const double dot = incoming.x * nx + incoming.y * ny;
    const double norm2 = nx * nx + ny * ny;

    return {
        incoming.x - 2.0 * dot * nx / norm2,
        incoming.y - 2.0 * dot * ny / norm2,
    };
}

Point next_intersection(const Point& current, const Point& direction) {
    const double num = -(8.0 * current.x * direction.x + 2.0 * current.y * direction.y);
    const double den = 4.0 * direction.x * direction.x + direction.y * direction.y;
    const double t = num / den;

    return {
        current.x + t * direction.x,
        current.y + t * direction.y,
    };
}

int solve() {
    Point previous{0.0, 10.1};
    Point current{1.4, -9.6};

    int reflections = 0;
    while (true) {
        const Point incoming{current.x - previous.x, current.y - previous.y};
        const Point outgoing = reflect_direction(incoming, current);
        const Point next = next_intersection(current, outgoing);

        ++reflections;
        previous = current;
        current = next;

        if (std::abs(current.x) <= 0.01 && current.y > 0.0) {
            break;
        }
    }

    return reflections;
}

bool run_checkpoints() {
    const Point previous{0.0, 10.1};
    const Point current{1.4, -9.6};
    const Point incoming{current.x - previous.x, current.y - previous.y};
    const Point outgoing = reflect_direction(incoming, current);
    const Point next = next_intersection(current, outgoing);

    if (!nearly_equal(next.x, -3.990597619361616, 1e-9) ||
        !nearly_equal(next.y, -6.024991498863843, 1e-9)) {
        std::cerr << "Checkpoint failed for first reflected intersection" << '\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() << '\n';
    return 0;
}

Python

# Problem 144: Investigating multiple reflections of a laser beam
# Laser bounces inside ellipse 4x²+y²=100 until it exits through opening |x|<=0.01, y>0.

import math

def solve():
    # Ellipse: 4x²+y²=100
    px, py = 0.0, 10.1  # previous point (start)
    cx, cy = 1.4, -9.6  # current point (first hit)
    bounces = 0
    
    while True:
        # incoming direction
        dx, dy = cx - px, cy - py
        # normal to ellipse at (cx, cy): gradient of 4x²+y²-100 = (8x, 2y)
        nx, ny = 8*cx, 2*cy
        nlen = math.sqrt(nx*nx + ny*ny)
        nx, ny = nx/nlen, ny/nlen
        # reflect: d' = d - 2*(d·n)*n
        dot = dx*nx + dy*ny
        rdx = dx - 2*dot*nx
        rdy = dy - 2*dot*ny
        # find next intersection with 4x²+y²=100
        # parametric: (cx+t*rdx, cy+t*rdy) on ellipse
        # 4(cx+t*rdx)²+(cy+t*rdy)²=100
        a_coef = 4*rdx*rdx + rdy*rdy
        b_coef = 2*(4*cx*rdx + cy*rdy)
        # c_coef = 4*cx*cx + cy*cy - 100 = 0 (already on ellipse)
        t = -b_coef / a_coef
        nx2 = cx + t*rdx
        ny2 = cy + t*rdy
        bounces += 1
        px, py = cx, cy
        cx, cy = nx2, ny2
        if abs(cx) <= 0.01 and cy > 0:
            break
    
    print(bounces)

solve()

Java

public class Euler144 {
    public static void main(String[] args) {
        double px = 0, py = 10.1, cx = 1.4, cy = -9.6;
        int bounces = 0;
        while (true) {
            double dx = cx - px, dy = cy - py;
            double nx = 8 * cx, ny = 2 * cy, nl = Math.sqrt(nx * nx + ny * ny);
            nx /= nl;
            ny /= nl;
            double dot = dx * nx + dy * ny;
            double rdx = dx - 2 * dot * nx, rdy = dy - 2 * dot * ny;
            double a = 4 * rdx * rdx + rdy * rdy, b = 2 * (4 * cx * rdx + cy * rdy);
            double t = -b / a;
            px = cx;
            py = cy;
            cx += t * rdx;
            cy += t * rdy;
            bounces++;
            if (Math.abs(cx) <= 0.01 && cy > 0)
                break;
        }
        System.out.println(bounces);
    }
}