<p>It is sort of expected that on numerical applications, Haskell code tends to be say 2X to 5X slower than equivalent code written in classic imperative languages such as C/C++/Fortran. However, Haskell code being 100 times slower is very unexpected !</p>
<h2>Sanity check: reproducing the result</h2>
<p>Using a somewhat oldish notebook with an Intel Celeron N2840 64 bits cpu, running Linux kernel 5.3, GLIBC 2.28, GCC 8.3, GHC 8.2.2, GHC random-1.1, we get the timings:</p>
<pre>
C code,        gcc -O3:    950 msec
Haskell code,  ghc -O3:    104 sec
</pre>
<p><strong>So indeed, with these configurations, the Haskell code runs about 100 times slower than the C code.</strong></p>
<h2>Why is that ?</h2>
<p>A first remark is that the π computing arithmetics atop random number generation looks quite simple, hence the runtimes are probably dominated by the generation of 2*10 millions pseudo-random numbers. Furthermore, we have every reason to expect that <strong>Haskell and C are using completely unrelated random number generation algorithms.</strong></p>
<p>So, rather than comparing Haskell with C, we could be comparing the random number generation algorithms these 2 languages happen to be using, as well as their respective implementations.</p>
<p>In C, the language standard does not specify which algorithm function <code>srand()</code> is supposed to use, but that tends to be old, simple and fast algorithms.</p>
<p>On the other hand, in Haskell, we have traditionally seen a lot of complaints about the poor efficiency of <code>StdGen</code>, like here back in 2006:</p>
<p><a href="https://gitlab.haskell.org/ghc/ghc/-/issues/427" rel="nofollow noreferrer">https://gitlab.haskell.org/ghc/ghc/-/issues/427</a></p>
<p>where one of the leading Haskell luminaries mentions that <code>StdGen</code> could possibly be made 30 times faster.</p>
<p>Fortunately, help is already on the way. This <a href="https://alexey.kuleshevi.ch/blog/2021/01/29/random-interface/" rel="nofollow noreferrer">recent blog post</a> explains how the <a href="https://hackage.haskell.org/package/random-1.2.0" rel="nofollow noreferrer"><strong>new Haskell random-1.2 package</strong></a> is solving the problem of <code>StdGen</code> lack of speed, using a completely different algorithm known as <em>splitmix</em>.</p>
<p>The field of pseudo-random number generation (PRNG) is a rather active one. Algorithms routinely get obsoleted by newer and better ones. For the sake of perspective, here is a relatively recent (2018) <a href="https://arxiv.org/pdf/1811.04035" rel="nofollow noreferrer">review paper on the subject</a>.</p>
<h2>Moving to more recent, better Haskell components:</h2>
<p>Using another, slightly more powerful machine with an Intel Core i5-4440 3GHz 64 bits cpu, running Linux kernel 5.10, GLIBC 2.32, GCC 10.2, and critically <strong>Haskell package random-1.2</strong>:</p>
<pre>
C code,        gcc -O3:    164 msec
Haskell code,  ghc -O3:    985 msec
</pre>
<p>So Haskell is now “just” 6 times slower instead of 100 times.</p>
<p>And we still have to address the injustice of Haskell having to use <code>x**2+y**2</code> versus C getting <code>x*x+y*y</code>, a detail which did not really matter before with <strong>random-1.1</strong>. This gives us 379 msec ! <strong>So we are back into the usual 2X-5X ballpark for Haskell to C speed comparisons.</strong></p>
<p>Note that if we ask the Haskell executable to run with statistics on, we get the following output:</p>
<pre>
$ time q66441802.x +RTS -s -RTS 10000000
3.1415616
     925,771,512 bytes allocated in the heap
     ...
  Alloc rate    2,488,684,937 bytes per MUT second

  Productivity  99.3% of total user, 99.3% of total elapsed

real   0m0,379s
user   0m0,373s
sys    0m0,004s
</pre>
<p>so Haskell is found to allocate close to one gigabyte of memory along the way, which helps to understand the speed difference.</p>
<h2>A code fix:</h2>
<p>We note that the C code uses a single random serie while the Haskell code is using two, with two calls to <code>mkStdGen</code> with numbers 1 and 2 as seeds. This is not only unfair but also incorrect. Applied mathematicians who design PRNG algorithms take great care to ensure any single serie has the right statistical properties, but give essentially no guarantees about possible correlations between different series. It is not unheard of to even use the seed as an offset into a single global sequence.</p>
<p>This is fixed in the following code, which does not alter performance significantly:</p>
<pre><code>computeNorms :: [Double] -&gt; [Double]
computeNorms []  = []
computeNorms [x] = []
computeNorms (x:y:xs2) = (x*x + y*y) : (computeNorms xs2)

monteCarloCircle :: Int -&gt; Double
monteCarloCircle nn =
    let
         randomSeed   =  1
         gen0         =  mkStdGen randomSeed
                      -- use a single random serie:
         uxys         =  (randomRs  (-1.0, 1.0)  gen0) :: [Double]
         norms        =  take nn (computeNorms uxys)
         insiderCount =  length  $  filter (&lt;= 1.0) norms
    in
         (4.0::Double) * ((fromIntegral insiderCount) / (fromIntegral nn))

</code></pre>
<h2>Side note:</h2>
<p>The new Haskell random-1.2 package has been mentioned in this <a href="https://stackoverflow.com/questions/66404604/using-system-random-over-multiple-data-types-the-haskell-way">recent SO question</a>, though in the context of the new monadic interface.</p>
<h2>A word of conclusion:</h2>
<p>Assuming the goal of the exercise is to gauge the relative runtime speeds of the C and Haskell languages, there are essentially two possibilities:</p>
<p>One is to avoid using PRNG altogether, because of the different algorithms in use.</p>
<p>The other one is to control random number generation by manually providing one and the same algorithm into both languages. For example, one could use the publicly available <a href="https://www.iro.umontreal.ca/%7Elecuyer/myftp/papers/streams00.pdf" rel="nofollow noreferrer">MRG32k3a</a> algorithm designed by Pierre L'Ecuyer. <a href="https://hackage.haskell.org/package/mrg-random-0.1.0.0/candidate" rel="nofollow noreferrer">Candidate Haskell MRG32k3a implementation here</a> (as of March 2021). Left as an exercise for the reader.</p>
Why is this Monte Carlo Haskell program so slow?
<p><em>update:</em></p>
<p>My compilation command is <code>ghc -O2 Montecarlo.hs</code>. My random version is <code>random-1.1</code>, ghc version is <code>8.6.4</code>, and my system is macOS Big Sur 11.1 (Intel chip). The command I used to test the speed is <code>time ./Montecarlo 10000000</code>, and the result it returns is <code>real 0m17.463s, user 0m17.176s, sys 0m0.162s</code>.</p>
<hr />
<p>The following is a Haskell program that uses Monte Carlo to calculate <code>pi</code>. However, when the input is 10 million, the program ran for 20 seconds. The C program written in the same logic took only 0.206 seconds. Why is this, and how can I speed it up? Thank you all.</p>
<p>This is the Haskell version:</p>
<pre><code>import System.Random
import Data.List
import System.Environment

montecarloCircle :: Int -&gt; Double
montecarloCircle x
   = 4*fromIntegral
        (foldl' (\x y -&gt; if y &lt;= 1 then x+1 else x) 0
           $ zipWith (\x y -&gt; (x**2 + y**2))
                     (take x $ randomRs (-1.0,1) (mkStdGen 1) :: [Double])
                     (take x $ randomRs (-1.0,1) (mkStdGen 2) :: [Double]) )
        / fromIntegral x

main = do
    num &lt;- getArgs
    let n = read (num !! 0)
    print $ montecarloCircle n
</code></pre>
<p>This is the C version:</p>
<pre><code>#include &lt;stdio.h&gt;
#include &lt;math.h&gt;
#include &lt;time.h&gt;
#include &lt;stdlib.h&gt;

#define N       10000000

#define int_t   long // the type of N and M

// Rand from 0.0 to 1.0
double rand01()
{
    return rand()*1.0/RAND_MAX;
}

int main()
{
        srand((unsigned)time(NULL));

        double x,y;
        int_t M = 0;

        for (int_t i = 0;i &lt; N;i++)
        {
                x = rand01();
                y = rand01();
                if (x*x+y*y&lt;1) M++;
        }
        double pi = (double)4*M/N;
        printf(&quot;%lf\n&quot;, pi);
}
</code></pre>