Skip to main content
SZCZ

Dark Alley Mathematics

Trouble in the alley

You're walking alone at night through a dark alley in the middle of the night[1] when a hooded figure stops you. They tell you: "Solve the following mathematical problem right now OR I WILL SHOOT YOU DEAD.". They give you a piece of paper with the problem statement.

Three points are chosen independently and uniformly at random from the interior of a unit circle. What is the probability that their circumcircle is entirely contained within the unit circle?

They have a gun so you have no choice but to solve it. Now, because you're super stressed you have no idea how to do this in some funky clever way by appealing to symmetry or whatever. You'll just have to do it the old fashioned way and painstakingly calculate the result. Luckily the hoodlum doesn't seem to care how long it takes you to arrive at the answer. He told you as much. He provides a notebook and pen.

The Cartesian start

You are dealing with 3 points, so 6 coordinates in total: 3 xs, 3 ys.

Now you know you're dealing with a circular problem (heh). So even though you're in grave danger (the ruffian keeps twirling the gun on his fingers as if to emphasize just how much danger) you've got the mental wherewithal to realise there's a potential variable transformation you can do to make it easier for yourself.

Let's express everything in terms of the circumcircle center O=(X,Y), its radius r as well as the three angles the 3 points make with the center O: ϕ1, ϕ2, ϕ3.

This gives you a mapping from (x1,y1,x2,y2,x3,y3) to (X,Y,r,ϕ1,ϕ2,ϕ3).

You have xi=X+rcosϕiyi=Y+rsinϕifor i=1,2,3

Now, you have to calculate the Jacobian since this is a non-trivial transformation, and you'll have to adjust the integral appropriately. The hoodlum reminds you the columns represent the partials with respect to the given output variable, and the rows are with respect to the same input variable. You proceed.

J=[x1Xx1Yx1rx1ϕ1x1ϕ2x1ϕ3y1Xy1Yy1ry1ϕ1y1ϕ2y1ϕ3x2Xx2Yx2rx2ϕ1x2ϕ2x2ϕ3y2Xy2Yy2ry2ϕ1y2ϕ2y2ϕ3x3Xx3Yx3rx3ϕ1x3ϕ2x3ϕ3y3Xy3Yy3ry3ϕ1y3ϕ2y3ϕ3]

The good news is that this is quite a regular Jacobian. The first two columns have a nice pattern, the third is straightforward (as long as you remember the cos and sin derivatives). The only slightly tricky columns are 4th, 5th and 6th - but they too have a cool pattern to them. You end up with:

J=[10cosϕ1rsinϕ10001sinϕ1rcosϕ10010cosϕ20rsinϕ2001sinϕ20rcosϕ2010cosϕ300rsinϕ301sinϕ300rcosϕ3]

Nice, you think to yourself. Now this is a 6x6 matrix. You really don't want to be calculating the determinant... But you need it. You begin trying to calculate the determinant but the masked figure stops you. They remind you of determinant properties that will make your job much easier:

  1. You can factor out a number from a whole row / column and multiply the determinant by the same factor instead.
  2. The determinant doesn't change if you add or subtract one row from another.
  3. For Block Triangular Matrices M=(AB0D) Where:

First, you notice the common factor of r in 3 columns. Factored out, they will give you a factor of r3 altogether. Now you have

J=[10cosϕ1sinϕ10001sinϕ1cosϕ10010cosϕ20sinϕ2001sinϕ20cosϕ2010cosϕ300sinϕ301sinϕ300cosϕ3] Now you think about the block triangular matrix thing the guy mentioned. Surely this is relevant here. You notice that you could clear the first two columns and turn the top left 2x2 section into an identity matrix. Then detA=1 which would make everything else much easier.

You continue. You want to clear the 1 in the first column, rows 3 and 5. Also the 1 in column 2, rows 4 and 6.

J=[10cosϕ1sinϕ10001sinϕ1cosϕ10000cosϕ2cosϕ1sinϕ1sinϕ2000sinϕ2sinϕ1cosϕ1cosϕ2000cosϕ3cosϕ1sinϕ10sinϕ300sinϕ3sinϕ1cosϕ10cosϕ3]

You're very happy. Since the top left A is now an identity matrix with detA=1 and you can skip the whole B submatrix in the top right. You've now boiled this all down to calculating:

detJ=detJ=detA×detD=detD

D is a 4×4 matrix, so not the easiest to calculate by hand. But it's tractable. Let's just remember to put that r3 factor out front at the end.

D=[cosϕ2cosϕ1sinϕ1sinϕ20sinϕ2sinϕ1cosϕ1cosϕ20cosϕ3cosϕ1sinϕ10sinϕ3sinϕ3sinϕ1cosϕ10cosϕ3]

You can now subtract row 1 from row 3, and subtract row 2 from row 4. You get D=[cosϕ2cosϕ1sinϕ1sinϕ20sinϕ2sinϕ1cosϕ1cosϕ20cosϕ3cosϕ20sinϕ2sinϕ3sinϕ3sinϕ20cosϕ2cosϕ3]

Now you can expand down colum 2 instead of the dense column 1.

det(D)=sinϕ1det(M12)cosϕ1det(M22) M12=[sinϕ2sinϕ1cosϕ20cosϕ3cosϕ2sinϕ2sinϕ3sinϕ3sinϕ2cosϕ2cosϕ3] M22=[cosϕ2cosϕ1sinϕ20cosϕ3cosϕ2sinϕ2sinϕ3sinϕ3sinϕ2cosϕ2cosϕ3] You go on. det(M12)=(sinϕ3)|sinϕ2sinϕ1cosϕ2sinϕ3sinϕ2cosϕ2|+(cosϕ3)|sinϕ2sinϕ1cosϕ2cosϕ3cosϕ2  sinϕ2|

det(M22)=(sinϕ3)|cosϕ2cosϕ1sinϕ2sinϕ3sinϕ2cosϕ2|+(cosϕ3)|cosϕ2cosϕ1sinϕ2cosϕ3cosϕ2  sinϕ2|

You start wondering if maybe you should have stuck with the original Cartesian product. Those sines look scary. You crack on nonetheless.

detJ=sinϕ1det(M12)cosϕ1det(M22)

And then, remembering about the r3 J=r3[sinϕ1(sinϕ1sin(ϕ3ϕ2)+cosϕ3cosϕ2)+cosϕ1(sinϕ2sinϕ3+cosϕ1sin(ϕ3ϕ2))]

Now, you know that cos2+sin2=1 and you also remember the formula for sin(AB). Also, the hoodlum helpfully points out that he's happy to give you a refresher on the sine of a difference in case you need it.

J=r3[(sin2ϕ1+cos2ϕ1)sin(ϕ3ϕ2)+(sinϕ1cosϕ3cosϕ1sinϕ3)+(cosϕ1sinϕ2sinϕ1cosϕ2)].

You're ecstatic. You love sines and cosines now. You loudly point out that sin2ϕ1+cos2ϕ1=1sinϕ1cosϕ3cosϕ1sinϕ3=sin(ϕ1ϕ3)cosϕ1sinϕ2sinϕ1cosϕ2=sin(ϕ2ϕ1)

and hence detJ=r3(sin(ϕ3ϕ2)+sin(ϕ1ϕ3)+sin(ϕ2ϕ1)).

Your thug gives you a supportive nod. You don't care about the sign of detJ so you flip the difference inside of the sines and add absolute value for good measure to remove the stupid minus sign. detJ=r3|sin(ϕ1ϕ2)+sin(ϕ2ϕ3)+sin(ϕ3ϕ1)|

At this point you notice that this kind of makes sense because three collinear points won't make any kind of circumcirle. You carry on but you make a mental note to remember this problem if you ever need to interview someone for a software role. This is trivial now after all.

Okay, the Jacobian is all sorted. You move on to the probability calculation. In general, since we're picking points from the unit circle (area π) and you do it 3 times ("independently" as the villain points out excitedly!) you have a "volume" of 1π3 to scale everything by.

Very generally P=1π31condition|J|,dX,dY,dr,dϕ1,dϕ2,dϕ3

and less generally

P=1π3dX,dY0drdϕ1,dϕ2,dϕ3|J|1condition

where the condition you're after is X2+Y2+r1 since you know you only care about points whose circumcirle fits nicely inside the unit circle. That is: R=X2+Y2 is the distance from the origin to the circumcircle, and you want to make sure that R+r1.

You hit this head on: you can split this into the angular part (the sine difference bit from the determinant of the Jacobian) and the geometric part.

For the geometric part, you realise that r can go from 0 to 1, and for any specific r the center (X,Y) must be within a smaller circle of radius (1r). "Area of the circle is pi times radius squared", the gangster points out. He's right.

So, the dX,dY0dr simplifies to 01r3π(1r)2dr (the r3 comes from that pesky |J| factor you promised not to forget. You did forget but the brigand helpfully reminded you to include it.).

Now, the angular part is just I=02πdϕ102πdϕ202πdϕ3|sin(ϕ1ϕ2)+sin(ϕ2ϕ3)+sin(ϕ3ϕ1)|

Now, you grab a spare bit of paper and try to tackle this triple angular integral that has cost you so much effort already.

You start scribbling on the second page to quickly work out this triple integral...

Let the integrand be S=sin(ϕ1ϕ2)+sin(ϕ2ϕ3)+sin(ϕ3ϕ1). Using the identity sinA+sinB+sin((A+B))=4sin(A/2)sin(B/2)sin((A+B)/2), you can rewrite the absolute value:

|S|=|4sin(ϕ1ϕ22)sin(ϕ2ϕ32)sin(ϕ1ϕ32)|

By shift invariance[2], you integrate over the differences u and v and multiply by the length of the interval 2π:

For more detail

Formally, you can perform a change of variables: {u=ϕ1ϕ2v=ϕ2ϕ3w=ϕ3 The Jacobian determinant of this transformation is 1, so dϕ1dϕ2dϕ3=du,dv,dw.

Substituting this into the integral: I=|S(u,v)|,du,dv,dw Notice that the integrand |S| only depends on u and v, not on w. This allows you to separate the integrals: I=(02πdw)×(|S(u,v)|,du,dv) The first integral simply evaluates to 2π. The remaining double integral accounts for the relative differences u and v over their periods.

I=2π02πdu02πdv|4sin(u2)sin(v2)sin(u+v2)|

Now, substitute x=u/2 and y=v/2. The limits become [0,π] and the differentials du,dv become 4,dx,dy:

I=2π440πdx0πdy|sin(x)sin(y)sin(x+y)|

The term sin(x+y) is positive when x+yπ and negative when x+yπ. Due to symmetry, the integral over these two triangular regions is identical. You calculate the positive region and multiply by 2.

I=32π20πdx0πxdy,sin(x)sin(y)sin(x+y)

I=64π0πsin(x)[0πxsin(y)sin(x+y),dy]dx

You use the product-to-sum formula sinAsinB=12[cos(AB)cos(A+B)] to integrate with respect to y.

Using the product-to-sum formula on the inner term:

0πx12(cos(x)cos(x+2y)),dy

=12[ycos(x)12sin(x+2y)]0πx

=12((πx)cos(x)+sin(x))

Substituting back into I:

I=64π0πsin(x)12((πx)cos(x)+sin(x)),dx

I=32π(0π(πx)sin(x)cos(x),dx+0πsin2(x),dx)

Evaluating the two integrals separately:

  1. 0π(πx)12sin(2x),dx=π4 (via integration by parts)
  2. 0πsin2(x),dx=π2

I=32π(π4+π2)=32π(3π4)

I=24π2

and you finally get the answer: 24π2.

You feel you're very close now. You go back to the integral expression for P and you plug in the angular integral value: P=1π3(24π2)01πr3(1r)2,dr

You simplify and expand the polynomial.

P=2401r3(12r+r2),dr

Distribute everything. The criminal seems very happy you're at this stage of the calculation.

P=2401(r32r4+r5),dr=24[r442r55+r66]01=24(1425+16)

Just the fractions now... P=24(3012048120+20120)=24(2120)=24(160)=2460=25

You've done it! The hooligan looks happy. He says you can go. He also says that everyone always calculates it this way, and he's yet to see a clever way that doesn't involve a sextuple integration. He invites you to share a clever solution later on by emailing it to him.

You wonder if you actually got the right answer. It was a nice fraction but did you not overcount? Undercount? Did you mess up the integral? The whole angular thing had like 6 permutations to it after all... Well, the villain seemed satisfied either way. You go home. You're very tired. You plan on not doing the "dark alley at night" thing any more.


  1. You don't remember why you were in this situation. ↩︎

  2. Since the integrand S depends only on the differences between angles, not their absolute values. ↩︎