Probability a random triangle on a circle contains the origin

Question

Another great video from 3Blue1Brown describes a neat proof (at about 5:45) that this has to be 25%. I wanted to do the experiment in Mathematica as a way of showing how exploring the problem computationally might give you a clue about the proof.

While I'm happy I got there in the end, the timing is really bad, about 10 seconds for 3000 triangles. Any hints?

RegionMember[#, 
   origin] & /@ (Table[{Cos[RandomReal[2 Pi]], 
       Sin[RandomReal[2 Pi]]}, 9000] // Partition[#, 3] & // 
    Map[Polygon, #] &) // Count[#, True] &

Answer 1

Some issues with your code are that the list of random data you generate is not a PackedArray; moreover creating random numbers is more efficient if you generate all at once like angles = RandomReal[2 Pi, {n, 3}]. Even if the data were packed, it would be unpacked when mapping Polygon to it. Moreover, testing against triangles with RegionMember is like "shooting sparrows with a cannon". In general: Map on big data is not so fast, even if it is considerably better than using For loops. Best efficiency is obtained with CompiledFunctions with the option RuntimeAttributes -> {Listable}.

The following CompiledFunction computes the signed areas spanned by the origin and each pair of corners of the input triangle and returns 1 if all the signs of these signed areas are equal (that's only possible if the origin lies within the triangle).

Quiet[Block[{P, PP},
  PP = Table[P[[i, j]], {i, 1, 3}, {j, 1, 2}] /. Part -> Compile`GetElement;
  With[{
    det3 = Det[{PP[[1]], PP[[2]]}],
    det1 = Det[{PP[[2]], PP[[3]]}],
    det2 = Det[{PP[[3]], PP[[1]]}]},
   testTriangle = Compile[{{P, _Real, 2}},
     Block[{s1, s2, s3},
      s1 = Sign[det1];
      s2 = Sign[det2];
      s3 = Sign[det3];
      Boole[s1 == s2 && s2 == s3]],
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
   ]
  ]]

On my machine, this can test about 25000000 triangles per second. That's about 10000 times faster than the original code. Note that generating the random triangles is more expensive than checking them.

n = 25000000;
triangles = With[{angles = RandomReal[2 Pi, {n, 3, 1}]},
   Join[Cos[angles], Sin[angles], 3]
   ];
N[Total[testTriangle[triangles]]/n] // AbsoluteTiming

{1.01982, 0.249998}

So, the experimental probability is pretty close to 25 %.

Answer 2

This is not efficient but I post it as another way and also for the visualization.

A query function:

qf = Compile[{{p, _Real}, {q, _Real}}, Module[{v = {p, q}, j, k, b},
    {j, k} = Sort@v;
    If[j > Pi, 0, Boole[Pi < k < j + Pi]]], CompilationTarget -> "C", 
   Parallelization -> True, RuntimeOptions -> "Speed"];

Simulation:

p[v_] := Module[{s = {#2 - #1, #3 - #1} & @@ (Sort@v)},
   qf @@ s];
nu[n_] := With[{rnd = RandomReal[2 Pi, {n, 3}]}, N[Total[p /@ rnd]/n]]

Visualization functions:

nug[v_] := 
 With[{res = 
    p@v}, {Graphics[{Circle[], PointSize[0.04], Point[{0, 0}], 
     RGBColor@RotateRight[{1, 0, 0}, res], 
     Polygon[{Cos@#, Sin@#} & /@ v]}], res}]
nugf[n_] := 
 With[{y = Transpose[nug /@ RandomReal[2 Pi, {n^2, 3}]]}, 
  Column[{N[Total@y[[2]]/n^2], 
    GraphicsGrid[Partition[y[[1]], n], Spacings -> {0, 0}, 
     ImageSize -> 800]}, Alignment -> Center]]

Some simulations:

Increasing sample size

After a long time (minutes) on my old machine:

Some visualizations:

• sample of 25 triangles

• sample of 100 triangles

• sample of 1600 triangles

Answer 3

This is a further extension of @Henrik's answer, a little bit too long for a comment.

When playing with Henric's great solution, I observed that the two settings Parallelization->True and Parallelization->False did not make the difference I hoped for. My parallel kernel configuration shows that automatically 4 local kernels will be launched, whereas the difference between two two settings was only a factor 1.5. So I tried to do the parallelization myself.

The next function is the function testTriangle of Henric, but with Parallelization set to False.

origininside=
  With[{Part=Compile`GetElement},
    Compile[{{tr, _Real,2}},
      Module[{s1, s2, s3},
        s1=Sign[tr[[2,1]]tr[[3,2]]-tr[[2,2]]tr[[3,1]]];
        s2=Sign[tr[[3,1]]tr[[1,2]]-tr[[3,2]]tr[[1,1]]];
        s3=Sign[tr[[1, 1]]tr[[2,2]]-tr[[1,2]]tr[[2,1]]];
        Boole[s1==s2==s3]],
      CompilationTarget->"C",
      Parallelization->False,
      RuntimeAttributes->{Listable}]
    ];

On my quadcore computer, I launch 8 kernels:

LaunchKernels[8];

I seem to remember that a C-compiled function is not immediately available in the subkernels, so I expected that now I had to do the compilation of the function origininside on each of the subkernels. Maybe that idea dates from a long time ago, maybe I am wrong, but I was pleasantly surprised to see that, having compiled the function in the main kernel, it is immediately avaible on the subkernels as well.

Therefore, with the following command I create and test 25,000,000 triangles:

nn=2500;
  N@Total[ParallelTable[
    Total[origininside @  With[{angles=RandomReal[2 Pi,{nn,3,1}]},
      Join[Cos[angles],Sin[angles],3]] ], {10000}]]/10000/nn // AbsoluteTiming

{0.770372, 0.250022}

Since this timing includes the time required for the construction of the triangles, this construction of the parallelization makes the computation about 2.5 times as fast.

Addendum

This addendum contains some mathematical refinements, including a proof that the probability indeed is 1/4.

The first refinement has to do with the construction of the triangles. There is no loss in generality when we assume the first pont of the triangle is at the top of the circle: A={0,1}. We then choose only two random points B={b1,b2} and C={c1,c2} on the unit circle. The condition that the origin is inside the triangle then amounts to b1 and c1 having opposite sign and the intersection of the line BC with the vertical axis is negative. The latter gives:

Assuming[b1 c1<0, Simplify[(z /. Solve[ (1-\[Lambda]){b1,b2}+\[Lambda]{c1,c2}=={0,z}, {\[Lambda], z}][[1]]) <0]]

(b1-c1) (-b2 c1+b1 c2)<0

which can be further simplified to

b1(-b2 c1+b1 c2)<0

The function origininside is now simpler:

origininside=With[{Part=Compile`GetElement},
  Compile[{{tr, _Real, 2}},
    If[tr[[1,1]]tr[[2,1]]>0, 0,
      Boole[tr[[1,1]](-tr[[1,2]] tr[[2,1]]+tr[[1,1]] tr[[2,2]])<0]],
    CompilationTarget->"C",
    Parallelization->False,
    RuntimeAttributes->{Listable}]]

The testing on 8 subkernels:

LaunchKernels[8]

(In the following command I use Henric's construction for random points on the circle; this turns out to be faster than using RandomPoint as mentioned by @Sjoerd Smit.)

nn=2500;
N@Total[ParallelTable[
  Total[origininside@With[{angles=RandomReal[2 Pi,{nn,2,1}]},
    Join[Cos[angles],Sin[angles],3]]],{10000}]]/10000/nn // AbsoluteTiming

{0.546887,0.25012}

This is about 25% faster.

The fact that the probability that the origin is inside a random triangle on the unit circle is 1/4 can be seen in the following way.

(Draw by hand a circle with a point A at the top and a point B at the right hand side. Let A' and B' be the opposite points of A and B. I would have placed this picture here in this answer, but I have no idea how to do that.)

The first point A is the top of the circle. Then we choose randomly the point B on the circle. Without loss of generality we may assume that B is at the right hand side of the circle, as in in the figure.

Now we choose the point C randomly on the unit circle. For the origin being inside the trangle, the point C must be on the arc A'B', This arc is as long as the arc AB.

Otherwise stated: we choose the point x randomly in the interval [0,π]. Then we choose y randomly in the interval[0,2π] and the criterion is that y in a subinterval of length x.

Therefore, the problem of finding the probability that the origin is inside a random triangle on the unit circle is equivalent with the following one:

What is the probability that when we randomly choose a number x in the interval [0,1], the random number yin the interval [0,1] will be less than x/2?

Obviously, the answer is

Integrate[x/2, {x, 0, 1}]

1/4

For the testing, we construct a list of 250000001 random numbers, so that we have 25000000 pairs. Then we can proceed in a lot ways, for example by partitioning the list into pairs with an increment of 1, and then count how many pairs have the property that the first element is larger than two times the second argument.

Here are two ways:

nn = 25000000; lst = RandomReal[1, {nn + 1}];

N[Total[UnitStep[ListCorrelate[{1, -2}, lst]]]/nn] // AbsoluteTiming

{0.45397, 0.250025}

lst1 = Most[lst];
Rest[lst];
N[Total[UnitStep[lst1 - 2 lst2]]/nn ] // AbsoluteTiming

{0.20445, 0.250025}

This is indeed very fast, without having used any C-compilation or parallel computing.

Answer 4

cf = With[{code = 
    Sign@Det@{{x1,y1},{x2,y2}}==Sign@Det@{{x2,y2},{x3,y3}}==Sign@Det@{{x3,y3},{x1,y1}}//Boole
    }, 
   Compile[{{A, _Real, 1}}, 
    Block[{x1 = A[[1]], x2 = A[[2]], x3 = A[[3]], y1 = A[[4]], y2 = A[[5]], y3 = A[[6]]},
     code
     ], RuntimeAttributes -> {Listable}
    ]
   ];

n = 10^7;

With[{angles = RandomReal[2 Pi, {n, 3}]},
  Join[Cos[angles], Sin[angles], 2]] // cf // Mean // N // AbsoluteTiming

{2.04003, 0.249968}

Module[{n, angles, x1, x2, x3, y1, y2, y3, s1, s2, s3},
    n = 10^7;
    angles = RandomReal[2 Pi, {3, n}];
    {x1, x2, x3} = Cos[angles];
    {y1, y2, y3} = Sin[angles];
    s1 = Sign[x1 y2 - x2 y1];
    s2 = Sign[x2 y3 - x3 y2];
    s3 = Sign[x3 y1 - x1 y3];
    Unitize[s1 - s2]~BitOr~Unitize[s2 - s3]
    ] // Mean // 1.0 - # & // AbsoluteTiming

{1.44645, 0.250057}

Answer 5

Since I wanted to try this out on a computer with very limited memory, generating the whole list of triangles all at once was out of the question. Instead, I use a Do[] loop to accumulate a running count of the triangles that passed the criterion, and then divide by the total number of samples to estimate the probability. Thus:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
            n = 1*^6; (* number of samples *) s = 0;
            Do[s += Boole[Graphics`PolygonUtils`InPolygonQ[
                          Transpose[Through[{Cos, Sin}[RandomReal[2 π, 3]]]], {0, 0}]],
               {n}];
            N[s/n, 20]]
   0.25003500000000000000

where I had used the undocumented point-in-polygon test featured in this answer.