(类似的讨论以往有过多次，不过当前老帖几乎全被百度隐藏，发出来权当补充吧。)
因为新电脑里预装了正版的Office，就稍微折腾了下VBA，四处翻找相关文章的时候，知乎的某主题引起了我的兴趣：
《为什么 C++ 只比 VBA 快 4倍？》（链接就不贴了，否则八成又要被吞。直接搜索这个标题就可以找到。）
那么Mathematica（在不改变算法实质的情况下）的速度怎么样？——基本可以肯定比C++慢（即便编译，Mathematica在大规模低级运算上也往往不及精心编写过的C/C++/Fortran），但和VBA比呢？（当年学VB时我曾不止一次听到过“VB速度慢”的说法，但我从没正经研究过。）
总之先试试C++和VBA的速度。C++快得意料之中，51毫秒（编译器是TDM-GCC 5.1.0-2 64位，编译选项-Ofast）：
VBA约279毫秒，和文中所述速度基本一致：
然后是Mathematica（测试版本9.0.1）：
（自认为）写的比较好的非编译版本：
(*
可以借助
Simplify`PWToUnitStep@PiecewiseExpand@If[x < 1, 1, 0]
计算出所需的UnitStep组合
*)
With[{a = RandomReal[1, {25 10^5}], b = RandomReal[1, {25 10^5}]},
1/(2.5 10^6) ((Total[1. - UnitStep[-1 + (a^2 + b^2)]] +
Total[1. - UnitStep[-1 + ((1 - a)^2 + (1 - b)^2)]])) 2] // AbsoluteTiming
(* {0.295614, 3.14008} *)
和VBA纯循环版本的速度差不多。
然后是编译版：
Compile[{},
With[{a = RandomReal[1, {25 10^5}], b = RandomReal[1, {25 10^5}]},
1/(2.5 10^6) ((Total[1. - UnitStep[-1 + (a^2 + b^2)]] +
Total[1. - UnitStep[-1 + ((1 - a)^2 + (1 - b)^2)]])) 2]][] // AbsoluteTiming
(* {0.126434, 3.14154} *)
比VBA快了，不过还是只有C++的一半。
让我有些惊讶的是，高维列表会导致代码略微变慢：
Compile[{},
With[{n = 25 10^5},
With[{lst = RandomReal[1, {2, n}]},
1./n ((Total[1. - UnitStep[-1. + Total[lst^2]]] +
Total[1. - UnitStep[-1. + Total[(1. - lst)^2]]])) 2]]][] // AbsoluteTiming
(* {0.152793, 3.14172} *)
编译为C的话又能榨出一点速度来，当然了，调用C编译器挺慢的……：
cf = Compile[{},
Module[{c = 0., n = 25 10^5},
Do[With[{a = RandomReal[], b = RandomReal[]}, If[a^2 + b^2 < 1., c += 1.];
If[(1. - a)^2 + (1. - b)^2 < 1., c += 1.];], {n}]; c 2/(n)],
CompilationTarget -> C, RuntimeOptions -> "Speed"];//AbsoluteTiming
(* {0.296880, Null} *)
cf[] // AbsoluteTiming
(* {0.110118, 3.14154} *)
上面的代码还有优化余地吗？

不编译的版本，使用Clip或Unitize比UnitSteps略快
With[{a=RandomReal[1,{25 10^5}],b=RandomReal[1,{25 10^5}],f=Clip[#,{1,0},{1,0}]&},
1/(2.5 10^6) (Total[f[a^2+b^2]]+Total[f[((1-a)^2+(1-b)^2)]]) 2]//AbsoluteTiming
(* f = Unitize[1, #] & *)

生成随机数的方法可以加快一点
SeedRandom[Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 3}}]

用在Compile里Random比RandomReal略快，只是生成的随机数可能质量不高，不太建议使用

之前回复“数值计算还是 MATLAB 比较快”被删了，理由是没给代码，那就给一下好了
上边每种方法都试过了（包括下边回复提到的改进方法），且都跑三次，取最快的一次, 所有方法中最快的是：
SeedRandom[Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 3}}]
Compile[{},
With[{a = RandomReal[1, {25 10^5}], b = RandomReal[1, {25 10^5}]},
1/(2.5 10^6) ((Total[1. - UnitStep[-1 + (a^2 + b^2)]] +
Total[1. -
UnitStep[-1 + ((1 - a)^2 + (1 -
b)^2)]])) 2]][] // AbsoluteTiming
(* {0.0415543, 3.14067} *)
（编译到 C 的方法在我这里反而更慢何况编译也挺慢就不贴了，编译器为 vs 2019）
环境是 win10 1909 x64, Mathematica 12.0.0, Intel 9900k @ 4.7GHz
相同环境的 MATLAB R2020a prerelease：
% foo.m
function foo
timeit(@璐村惂鐢ㄦ埛_05NCUJN馃惥 computePi)
function p = computePi
a = rand(25e5, 1);
b = rand(25e5, 1);
p = (sum(a .^ 2 + b .^ 2 < 1) + sum((1 - a) .^ 2 + (1 - b) .^ 2 < 1)) / 25e5 * 2;
>> rng(0, 'threefry')
>> foo
ans =
0.0108385555
用时大约是 Mathematica 的四分之一
考虑到这里随机数生成器不同，单独对比下随机数生成性能：
SeedRandom[
Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 3}}]
a = RandomReal[1, {25 10^5}]; // AbsoluteTiming
(*{0.004575, Null}*)
>> rng(0, 'threefry'); timeit(@() rand(25e5, 1))
ans =
0.0038504555
可以看到两者随机数生成性能是差不多的，所以应该主要差别在计算上：
0.0415543-0.004575*2 = 0.0324043
0.0108385555 - 0.0038504555*2 = 0.0031376445
这样估算，排除随机数生成器的影响后，计算性能可能相差 10 倍以上。
当然也可以写个对于提前计算得到的 a,b 进行计算的实现进行比较：
function foo
a = rand(25e5, 1);
b = rand(25e5, 1);
timeit(@() computePi(a, b))
function p = computePi(a, b)
p = (sum(a .^ 2 + b .^ 2 < 1) + sum((1 - a) .^ 2 + (1 - b) .^ 2 < 1)) / numel(a) * 2;
>> foo
ans =
0.00272264175
比上边计算出来的 0.0031376445 还快一些

计算卷积
欧拉计划249，matlab代码很难在10秒内完成，zhuanlan.zhihu.com/p/53416615