Problem 144 of Project Euler is once again a geometry problem, just like the previous. However, it is completely different. The problem reads

In laser physics, a “white cell” is a mirror system that acts as a delay line for the laser beam. The beam enters the cell, bounces around on the mirrors, and eventually works its way back out.

The specific white cell we will be considering is an ellipse with the equation 4x^{2}+y^{2}= 100

The section corresponding to -0.01≤x≤+0.01 at the top is missing, allowing the light to enter and exit through the hole.

The light beam in this problem starts at the point (0.0,10.1) just outside the white cell, and the beam first impacts the mirror at (1.4,-9.6).

Each time the laser beam hits the surface of the ellipse, it follows the usual law of reflection “angle of incidence equals angle of reflection.” That is, both the incident and reflected beams make the same angle with the normal line at the point of incidence.

In the figure on the left, the red line shows the first two points of contact between the laser beam and the wall of the white cell; the blue line shows the line tangent to the ellipse at the point of incidence of the first bounce.

The slopemof the tangent line at any point (x,y) of the given ellipse is:m= -4x/y

The normal line is perpendicular to this tangent line at the point of incidence.

The animation on the right shows the first 10 reflections of the beam.

How many times does the beam hit the internal surface of the white cell before exiting?

We will simply brute force our way through this problem, by calculating the laser beams path through the cell, and check if it hits the exit. In order to do that, we need to calculate how the laserbeam reflects. Once we know the angle of the reflecting beam, we can calculate the corresponding line, since we have the point of reflection. Once we have a line parameterization of the reflecting line, it is simply a matter of finding out where the line and the ellipse intersect. This will be the next point out laser beam hits. Confused yet? Don’t be. I will elaborate on it. Lets start by finding the slope of the reflecting beam.

## Slope of reflecting beam

Let us draw the situation and put some names on the different angles.

As you can see we have a horizontal line which we will base everything off. Then we have the ellipse tangent where the slope (which we will denote ) is already know from the problem description. Then we have the line AO, which we also know the slope of (denoted ). And then we have the line OB which is the one we are interested in finding the slope of (denoted ). So now we have already established a relation between the slopes and the angles.

One thing we can deduce from the image is that

and

We are interested in rewriting these relations into something based on the slopes of the lines.. Before continue let me just mention that we have

A relation which we will use in just a few seconds. For the first relation we have

For the second relation we get that

In this equation we can isolate such that we have

where tan(a) is given by the first equation.

So know we have the slop of the line OB, and since we know a point on the line (O) we can describe it as

where is the slope and d is the intercept. Which having a point becomes

## Intersecting a line and an ellipse

Now that we have a full description of the line OB all we need is to figure out where the line intersects with the ellipse. The ellipse have the equation

Plugging the line equation into this gives us that

So all we have to do is solve this equation and take the one of the two solutions which is not , then we have the x value of B and the y value follows from the line equation.

## Coding it all

For once I have actually managed to comment my C# code, so you should be able to read it based on the problem description. At least I wont comment further on it.

int result = 0; double xA = 0.0; double yA = 10.1; double xO = 1.4; double yO = -9.6; while(xO > 0.01 || xO < -0.01 || yO < 0){ //Calculate the slope of A double slopeA = (yO - yA) / (xO - xA); //Calculate the slope of the ellipse tangent double slopeO = -4*xO/yO; //Calculate the slope of B double tanA = (slopeA - slopeO)/(1+slopeA*slopeO); double slopeB = (slopeO- tanA)/ (1+ tanA*slopeO); //calculate intercept of line B double interceptB = yO - slopeB * xO; //solve the quadratic equation for finding // the intersection of B and the ellipse // a*x^2 + b*x + c = 0 double a = 4 + slopeB*slopeB; double b = 2 * slopeB * interceptB; double c = interceptB * interceptB - 100; double ans1 = (-b + Math.Sqrt(b * b - 4 * a * c)) / (2 * a); double ans2 = (-b - Math.Sqrt(b * b - 4 * a * c)) / (2 * a); xA = xO; yA = yO; //Take the solution which is furtherst from x0 xO = (Math.Abs(ans1 - xO) > Math.Abs(ans2 - xO)) ? ans1 : ans2; yO = slopeB * xO + interceptB; result++; }

This gives me the answer in 0.02ms or so. So this is by far fast enough. As always you can find the full source code here.

Hi.

One suggestion and one question regarding this post, which I now mostly understand, after many failed attempts over the past few years.

First the suggestion.

When you are explaining tab(a-b), I thought that ‘a’ was referring to the ‘a’ in the diagram and earlier calculations, not just a generic variable in the tab(a-b) formula, and it wasn’t for a while that it twigged. It might be worth explaining that somewhere.

Secondly, the question.

I’m probably being dumb, but there again as I type this I’m still half asleep, but I can see how you extract a=(a1+pi)-a2 from the diagram, but not a=a0-a1. How do you do this?

Hi Scott

I take the first part as a huge compliment, so thank you. For the second part, I have changed it to p and q to avoid confusing now. Thanks for that suggestion.

The question you have. Take a look at the illustration again. It is apparent in the right side of the illustration

a is the angle between the ellipse tangent and the laser beam

a0 is the line between the laser beam and the horizontal line

a1 is the angle between the ellipse tangent and the horizontal line

so if we take the a0 – a1 then we get

(laser beam – horizontal) – (ellipse – horizontal) -> laserbeam – horizontal

which is exactly the angle a.

Hi Kristian

🙂

My attempts in the past have all ended in failure, either from taking so long I killed it, giving me the wrong answer, or working only if I multiple something by -1 sometimes and not others, something I could never understand. I still have the code floating around somewhere (One of my old attempts is http://pastebin.com/BiHj1ACu). Yet this solution using tan(a-b) just works and is easy to understand.

Using your explanation of a=a0-a1, I can now see where I was going wrong when trying to understand the diagram you provided.

Firstly, I was only paying attention to the a on the left, not noticing a on the right which is the same angle due to the nature of the problem.

Secondly, I was only seeing a0 as the portion between the laser beam and the ellipse tangent (which happens to actually be the right a), not realizing the mark on the diagram extended through the ellipse tangent to the horizontal line.

Table 1

Discrete variations in the limits of the parameters j and k and its effects on the number of reflections of the laser beam.

White cell equations:

4*(x^2) + y^2 = 100 [1]

j*(x^2) + y^2 = k [2]

http://img4.hostingpics.net/pics/556817pe144ihtable1fig12.jpg

As j is kept constant at 4, we observe that

during the variation of k from 99 to 100.2, refl.

will vary between 19 and 896. The value of 354

reflections is reached when k is closer to 100.

Similar type of observations when k=100 and j varies

between 3.9 and 4.2. The value of 354 reflections is

obtained when j is closer to 4.

___

Sources:

1) Project Euler 144

Kristian’s algorithm

http://www.mathblog.dk/files/euler/Problem144.cs

2) Microsoft Visual C# 2010 Express

3) Microsoft Office Excel (2007)

//Take the solution which is furthest from x0

xO = (Math.Abs(ans1 – xO) > Math.Abs(ans2 – xO)) ? ans1 : ans2;

yO = slopeB * xO + interceptB;

Why you took furthest ?

There are two intersections between the oval and the line. One of these intersections is the one where the beam is hitting after the reflection. The other is the point of reflection.

Due to numerical issues I cannot be sure that x0 does not vary slightly from the answer. So instead I just take the point which is farthest away.

Hi Kristian,

You don’t need to solve a quadratic equation. Since you know (xO,yO) is a solution, you can factor out x-xO and remain with a linear equation for the other solution. I programmed it and it works.

https://sites.google.com/site/orenlivne/

Oren

That sounds right. However, I cannot wrap my head around how to do it immidatly. Can you show how you did it?

Hi Kristian,

Since OB passes through (xO,yO), write it as y=c(x-xO)+yO. Hence

4x^2 + (c(x-xO)+yO)^2 = 100

4x^2 + c^2(x-xO)^2 + 2c(x-xO)yO + yO^2 – 100 = 0

But yO^2 – 100 = -4xO^2, since (xO,yO) also lies on the ellipse.

4x^2 + c^2(x-xO)^2 + 2c(x-xO)yO – 4xO^2 = 0

4(x^2-xO^2) + c^2(x-xO)^2 + 2c(x-xO)yO = 0

4(x-xO)(x+xO) + c^2(x-xO)^2 + 2c(x-xO)yO = 0

Now we can finally factor out x-xO.

4(x+xO) + c^2(x-xO) + 2c yO = 0

(c^2+4)x + 4 xO-c^2 xO + 2c yO = 0

x = (xO(c^2-4) – 2c yO)/(c^2+4).

Please double-check this computation as I have a tendency to make mistakes. Thank you very much for putting this website together. It saved me many sleepless nights over several of the Euler problems.

Ah yes, now I see it. And that means we do indeed reduce x to be a linear equation based on values we already know. Nice one.

I tried this problem instead in parametric form of an ellipse (x=acos(t),y=bsin(t)) and it’s MUCH simpler. No need of a quadratic equation. Finding the third parametric angle if first two parametric angles are known doesn’t even need inverse trigonometric functions. It’s almost elementary after you derive the expression.

It seems to me that you went of on a tangent here. (excuse the pun)

This problem however is very well suited to be solved in vector form.

Let u_k be the direction of the incomming ray and x_k be the point where it strikes the ellipse.

u_0 = (1.4, -19.7), x_0 = (1.4, -9.6)

Since we are told the formula for the slope of the tangent is -4*x/y it is easy to show that the direction of the normal vector as point x_k can be written as

n_k = (-4*x_k[0], -x_k[1])

By the law or reflection the vector differnce between the incoming and outgoing rays must be a multiple of the normal vector:

u_(k+1) = u_k + p*n_k

where p is an unknown scalar.

The fact that both rays must make the same angle with the tangent line is equivalent to requirering that thier inner products with the normal vector be of same size but opposite sign. This can be used to solve for p:

p = -2*dot(u_k, n_k)/dot(n_k, n_k)

Next we need to find the new point where the outgoing ray strikes the ellipse. This point must lie on the line:

x_(k+1) = x_k + q*u_(k+1)

where q is an unknown scalar.

Define a matrix:

A = [[2,0],[0,1]]

Note that the function f(x_k) = 4*x_k[0]^2 + x_k[1]^2

can be written as: f(x_k) = dot(A*x_k, A*x_k)

Therfore:

f(x_(k+1)) = dot(A*(x_k + q*u_(k+1)), A*(x_k + q*u_(k+1)))

= dot(A*x_k + q*A*u_(k+1), A*x_k + q*A*u_(k+1))

= dot(A*x_k, A*x_k) + 2*q*dot(A*x_k, A*u_(k+1))

+ q^2*dot(A*u_k(k+1), A*u_(k+1))

The first term is the function value of the current point, which should be unchanged. Therefore the sum last two terms must equal zero.

The trivial solution, q = 0, is the current point.

The non-trivial solution gives the new point:

q = -2*dot(A*u_k(k+1), A*u_(k+1))/dot(A*x_k, A*u_(k+1))

This became a lot of text, but for the code it boils down to needing just five lines inside the loop:

calculate n_k at x_k

calculate p

update u_k

calculate q

update x_k