# vectorized operations is not efficient?

Consider the following example:

4.*Mean@Table[ Boole[Random[]^2 + Random[]^2 < 1], {10^6}] // AbsoluteTiming
(*{0.408023, 3.14105}*)

With[{n = 10^6},
4. Mean@Boole@Thread[(RandomReal[1, n]^2 + RandomReal[1, n]^2) < 1]] // AbsoluteTiming
(*{1.460084, 3.14276}*)


I expected the second method to be more efficient since it used vectorized operations, but it works slowly, and I don't know why. Is it possible to make it faster? I used Mathematica 9 on Windows 7.

First, the old Random function uses an inferior PRN generator and should not be used other than for legacy compatibility, if I recall correctly.

On my machine (version 7, Windows 7) the second method is faster than the first, though not by a great amount:

4.*Mean@Table[Boole[RandomReal[]^2 + RandomReal[]^2 < 1], {10^6}] // AbsoluteTiming

With[{n = 10^6},
4. Mean@Boole@Thread[(RandomReal[1, n]^2 + RandomReal[1, n]^2) < 1]] // AbsoluteTiming


{1.2840735, 3.14154}

{0.4570262, 3.1401}

Numeric methods will be much faster:

With[{n = 10^6},
With[{r = RandomReal[1, n]^2 + RandomReal[1, n]^2},
4. * Mean[1 - UnitStep[r - 1]]
]
] // AbsoluteTiming


{0.0510029, 3.14232}

This is partly because the array is not unpacked, as it is in your code as can be seen by first running On["Packing"]:

On["Packing"]

With[{n = 10^6},
4. Mean@Boole@Thread[(RandomReal[1, n]^2 + RandomReal[1, n]^2) < 1]] // AbsoluteTiming


DeveloperFromPackedArray::unpack: Unpacking array in call to Less. >>

DeveloperFromPackedArray::punpack1: Unpacking array with dimensions {1000000} to level 1. >>

• "First, the old Random function uses an inferior PRN generator and should not be used other than for legacy compatibility, if I recall correctly." - right, the old Marsaglia-Zaman method (which is now also available with SeedRandom[seed, Method -> "Legacy"]) has been shown to be statistically inferior to the current default method. The OP might also consider switching to a different PRNG if need be, using also SeedRandom[]. – J. M. is away May 18 '13 at 13:47
• Thank you. Is this right? With[{n = 10^6}, With[{r = RandomReal[1, n]^2 + RandomReal[1, n]^2}, 4.*Mean[UnitStep[1 - r]]]] – expression May 19 '13 at 2:55
• @explorer I believe you are asking: is UnitStep[1 - r] equivalent to 1 - UnitStep[r - 1]? Not precisely, because of the case of 1: {UnitStep[1 - r], 1 - UnitStep[r - 1]} /. r -> 1. It is unlikely to make a difference in your application however, IMHO, but to be rigorous I used the longer form. – Mr.Wizard May 19 '13 at 5:49
• @Mr.Wizard Thanks again for being so helpful. – expression May 19 '13 at 8:14

Can't compete with Mr Wizard, but here's a vectorized version that's not too shabby:

n = 10^6;
-Mean[Sign@Clip[Total[RandomReal[1, {2, n}]^2] - 1, {-100, 0}]] // AbsoluteTiming


{0.133305, 3.14322}