On e-day, I came across this cool tweet from Fermat’s library
Happy e-day of the century! 2 7 18
Curious fact about Euler's Number: Pick a uniformly random number in [0,1] and repeat until the sum of the numbers picked is >1. You'll on average pick e≈2.718… numbers! pic.twitter.com/P0wx9hQu79
— Fermat's Library (@fermatslibrary) February 7, 2018
So I spend a few minutes coding this into Julia
function euler(n) m=0 for i=1:n the_sum=0.0 while true m+=1 the_sum+=rand() (the_sum>1.0) && break; end end m/n end
Timing this on my machine, I got
julia> @time euler(1000000000) 15.959913 seconds (5 allocations: 176 bytes) 2.718219862
Gave a little under 16 seconds.
Tried a c
implementation
#include <stdio.h> /* printf, NULL */ #include <stdlib.h> /* srand, rand */ #include <time.h> /* time */ double r2() { return (double)rand() / (double)((unsigned)RAND_MAX + 1); } double euler(long int n) { long int m=0; long int i; for(i=0; i<n; i++){ double the_sum=0; while(1) { m++; the_sum+=r2(); if(the_sum>1.0) break; } } return (double)m/(double)n; } int main () { printf ("Euler : %2.5f\n", euler(1000000000)); return 0; }
and compiling with either gcc
gcc -Ofast euler.c
or clang
clang -Ofast euler.c
gave a timing twice as long
$ time ./a.out Euler : 2.71829 real 0m36.213s user 0m36.238s sys 0m0.004s
For the curios, I am using this version of Julia
julia> versioninfo() Julia Version 0.6.3-pre.0 Commit 93168a6 (2017-12-18 07:11 UTC) Platform Info: OS: Linux (x86_64-linux-gnu) CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz WORD_SIZE: 64 BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell) LAPACK: libopenblas64_ LIBM: libopenlibm LLVM: libLLVM-3.9.1 (ORCJIT, haswell)
Now one should not put too much emphasis on such micro benchmarks. However, I found this a very curious examples when a high level language like Julia
could be twice as fast a c
. The Julia
language authors must be doing some amazing mojo.
You are benchmarking rand(), not C vs. Julia. Oh, and don’t use rand().
Indeed. Please check my follow on post.