Return to Answer

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

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

This is another way:

nn = 2500000;
(lst = RandomReal[1, nn + 1];
  N[Total[UnitStep[ListCorrelate[{1, -2}, lst]]]/nn]) // AbsoluteTiming

{0.0656782, 0.249927}

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}

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.

This is another way:

nn = 2500000;
(lst = RandomReal[1, nn + 1];
  N[Total[UnitStep[ListCorrelate[{1, -2}, lst]]]/nn]) // AbsoluteTiming

{0.0656782, 0.249927}

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