Xrfgtjtk

I am testing the speed of Mathematica's latest version against a VB macro for simple loops to evaluate Pi.

So the VB macro is:

c = 0: n = 1000000: m = n * n



For i = 0 To n



    For j = 0 To n



        If (i * i + j * j) / m < 1 Then



            c = c + 1



        End If



    Next



Next



Selection.Text = 4 * c / m

I know how to write the corresponding code in Mathematica, but I do not know the intricacies of Mathematica regarding faster program execution. Can you help please by giving me an example of best optimized code for this example?

There are approximately 4 trillion operations in this example and the VB macro takes less than a day to execute...

This program partitions a quarter of a circular area of radius = 1 to smaller squares of area $1/n^2$ then counts them and adds them together. Essentially, a Monte Carlo Method but with randomness removed.

With your help I have created this faster code which uses parallel computing:

arg = Range[ 0, 1000000]; calcCompiled = 

Compile[{i}, 

Module[{c = 0}, Do[If[(i*i + j*j)/1000000^2 < 1, ++c], {j, 1000000}]; 

c], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 

Parallelization -> True, 

RuntimeOptions -> "Speed"];

AbsoluteTiming@N[4 Total[calcCompiled[arg]]/1000000^2, 20]

Any suggestions how to make this faster are welcome....

This is an algorithm of complexity $O(n)$ instead of $O(n^2)$. It runs through in about 8 milliseconds.

n = 1000000;

piapprox = 

    Total[Floor[

       Sqrt[Ramp[Subtract[N[n^2] - 1., Range[1., n]^2]]]]] (4./n^2); //

   RepeatedTiming // First



piapprox - Pi

0.00812

-4.00401*10^-6

Moreover, you get more bang for your bucks with the trapezoidal rule as follows:

n = 1000000;

weights = ConstantArray[1./n, n + 1];

weights[[{1, -1}]] = 0.5/n;

piapprox2 = 4. (weights.Sqrt[1. - Subdivide[0., 1., n]^2]);

piapprox2 - Pi

-1.17598*10^-9

And this computes Pi with 16-digit precision by approximating the intersection of the unit disk with the first quadrant by $2^k$ congruent triangles; the area of a single triangle is computed (one half of the result of Det) and multiplied by $4, 2^k$.

RepeatedTiming[

  iters = 25;

  {p, q} = N[IdentityMatrix[2]];

  piapprox = 2^(k + 1) Det[{p, Nest[x [Function] Normalize[p + x], q, iters]}];

  ][[1]]

piapprox - Pi

0.000048

4.44089*10^-16

This uncompiled versions needs about 0.04 milliseconds.

You're comparing a compiled language with an interpreted one. Or, at least, you are using the Wolfram language in interpreted mode. You'll have better results if you use a Compiled function.

Example without Compile — directly translated code:

calcNaive = Function[{n},

   Module[{c = 0., m = N[n^2]},

    Do[

     Do[

      If[(i*i + j*j)/m < 1.,

       ++c

       ]

      , {i, 0., n}]

     , {j, 0., n}];

    4. c/m]];

AbsoluteTiming@calcNaive[10000]

{264.352184, 3.14199016}

First value in the output is timing in seconds, second is the function return value.

Note here the N[n^2] instead of simple n^2. This way we avoid using arbitrary-length integers in further code, which would result in a 17% longer computation for no real benefit. In general, if you don't need advanced precision, you'd better enter all your numbers as machine-precision numbers, e.g. 1000.0 instead of 1000.

Now the code using Compile (which always uses machine-precision numbers, in addition to C-language-compiler optimizations):

calcCompiled = Compile[{n},

   Module[{c = 0, m = n^2},

    Do[

     Do[

      If[(i*i + j*j)/m < 1,

       ++c

       ]

      , {i, 0, n}]

     , {j, 0, n}];

    4 c/m]

    , CompilationTarget -> "C"

    , RuntimeOptions -> "Speed"];

AbsoluteTiming@calcCompiled[10000]

{0.918117, 3.14199016}

Notice more than two orders of magnitude faster calculation.

Here's another comparison. First let's look at the raw, inefficient version that someone who comes to Mathematica from another language might try:

badForMathematica[n_] :=

 Module[{i, j, c = 0., m = n^2},

  Do[

   If[(i^2 + j^2)/m < 1,

    c = c + 1

    ],

   {i, 0, n - 1},

   {j, 0, n - 1}

   ];

  (4*c)/m

  ]



badForMathematica[1000] // AbsoluteTiming



{3.07878, 3.14552}

Eeesh. But this is the dumb way. Here's a new approach. First we take the your algorithm and realize you're just subdividing and counting hits inside the circle. We can do this much faster using vectorized operations. Here's a way to get all the hits:

countHits[{i1_, i2_}, {j1_, j2_}, m_] :=



  With[{r1 = Range[i1, i2]^2, r2 = Range[j1, j2]^2},

   Total@UnitStep[m - Join @@ Outer[Plus, r1, r2]]

   ];

We generate the Outer product you're really generating, then totaling the counts. Simple as that. Now we cook this into a larger algorithm where we work in chunks because I don't have enough memory to work all at once:

goodForMathematica[n_, chunkSize_: 1000] :=

 Module[{m, c = 0, chunks},

  m = n^2;

  chunks =

   If[chunkSize >= n,

    {{0, n}},

    If[#[[-1]] =!= n, 

       Append[#, n],

       #

       ] &@Range[0, n, chunkSize]

    ];

  Do[c += countHits[c1, c2, m], {c1, chunks}, {c2, chunks}];

  (4.*c)/m

  ]

If you have lots of memory just do it straight out (but you probably don't have enough). Here's that test:

goodForMathematica[10000] // AbsoluteTiming



{3.37823, 3.14199}

We increased the count by an order of magnitude but had similar performance. Note that this algorithm has trash scaling, so we can expect that n = 20000 will take ~12s because we've got four blocks to work with:

1000000 will therefore take about ~35000s or ~10 hours. Not gonna actually do it.

Finally, here's probably the first thing to think of, which is to just use Compile:

compiledVersion =

  Compile[{{n, _Integer}},

   Module[{i, j, c = 0., m = n^2},

    Do[

     If[(i^2 + j^2)/m < 1,

      c = c + 1

      ],

     {i, 0, n - 1},

     {j, 0, n - 1}

     ];

    (4*c)/m

    ],

   RuntimeOptions -> "Speed",

   CompilationTarget -> "C"

   ];



compiledVersion[100000] // AbsoluteTiming



{14.8691, 3.14163}

And this blows everything else out of the water (notice I use 100000 per grid subdivision), but that makes sense, because we're really using C, not Mathematica.