Qfyilyi

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

Now 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.

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, seconds 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.

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.

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.

Initialize the variables

c = 0;

n = 1000000;

Below is a more or less literal translation of the Visual Basic code.

check[i_, j_] := If[(i*i + j*j)/(n*n) < 1, c++]



Timing[

 Scan[

  Function[i,

   Scan[

    Function[j,

     check[i, j]

     ],

    Range[n]

    ]

   ],

  Range[n]

  ]

 ]

A slightly faster verion gets rid of the overhead caused by the check function and the computation of Range[n] as well as n*n (suggested by dimachaerus comment).

nList = Range[n];

n2 = n*n;



Timing[

 Scan[

  Function[i,

   Scan[

    Function[j,

     If[(i*i + j*j)/n2 < 1, c++]

     ],

    nList

    ]

   ],

  nList

  ]

 ]

The answers from the two experiments are: (to be inserted when complete)