<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Mon, Apr 3, 2017 at 3:20 PM, Pierre Haessig <span dir="ltr"><<a href="mailto:pierre.haessig@crans.org" target="_blank">pierre.haessig@crans.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
  
    
  
  <div bgcolor="#FFFFFF">
    <p>Hello, <br>
    </p>
    Le 30/03/2017 à 13:31, Pierre Haessig a écrit :<br>
    <blockquote type="cite">
      
      [....]
      <span class="gmail-"><p>But how come Julia is 4-5x faster since Numpy uses C
        implementation for the entire process ? (Mersenne Twister ->
        uniform double -> Box-Muller transform to get a Gaussian
        <a class="gmail-m_5503631979399465996moz-txt-link-freetext" href="https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/randomkit.c" target="_blank">https://github.com/numpy/<wbr>numpy/blob/master/numpy/<wbr>random/mtrand/randomkit.c</a>).
        Also I noticed that Julia uses a different algorithm (<span class="gmail-m_5503631979399465996pl-c">Ziggurat Method from Marsaglia and Tsang</span> ,
        <a class="gmail-m_5503631979399465996moz-txt-link-freetext" href="https://github.com/JuliaLang/julia/blob/master/base/random.jl#L700" target="_blank">https://github.com/JuliaLang/<wbr>julia/blob/master/base/random.<wbr>jl#L700</a>)
        but this doesn't explain the difference for uniform rng.<br>
      </p>
    </span></blockquote>
    Any ideas?<br></div></blockquote><div><br></div><div><a href="https://github.com/JuliaLang/julia/blob/7fb758a275a0b4cf0e3f4cbf482c065cb32f0011/doc/src/stdlib/numbers.md#L116">This</a> says that Julia uses <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT">this library</a>, which is different from the home brewed version of the Mersenne twister in NumPy. The second link I posted claims their speed comes from generating double precision numbers directly, rather than generating random bytes that have to be converted to doubles, as is the case of NumPy through this <a href="https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/randomkit.c#L514">magical incantation</a>. They also throw the SIMD acronym around, which likely means their random number generation is parallelized.</div><div><br></div><div>My guess is that most of the speed-up comes from the SIMD parallelization: the Mersenne algorithm <a href="https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/randomkit.c#L221">does a lot of work</a> to produce 32 random bits, so that likely dominates over a couple of arithmetic operations, even if divisions are involved.</div><div><br></div><div>Jaime</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div bgcolor="#FFFFFF"><p>Do you think Stackoverflow would be a better place for my
      question?</p>
    <p>best,</p>
    <p>Pierre</p>
  </div>

<br>______________________________<wbr>_________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@python.org">NumPy-Discussion@python.org</a><br>
<a href="https://mail.python.org/mailman/listinfo/numpy-discussion" rel="noreferrer" target="_blank">https://mail.python.org/<wbr>mailman/listinfo/numpy-<wbr>discussion</a><br>
<br></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">(\__/)<br>( O.o)<br>( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.</div>
</div></div>