8 ms·

An example of fast numerical computation using Haskell

• xanth 13y ago
  I think the functional languages are always going to be slower than procedural languages for most task, because they are an abstraction away from the bare metal. However the ease of maintaining functional code may out way the 19% (and diminishing) faster procedural code.
  • friendly_chap 13y ago
    Unless bare metal changes to something which matches the functional approach better. Or just simply changes and higher level concepts protect you because you are not tied to implementation details.
  • toolslive 13y ago
    C is indeed close to hardware. 1970s hardware that is. That being said, the rule of thumb for Haskell,OCaml,... is that in terms of performance they are somewhere between 2 times slower and 2 times faster, depending on details, strategy, etc. In this particular case, it seems that the runs might be too short to see a significant difference. btw, try -O3 iso -O2. it makes a difference...
    • bitkrieg 13y ago
      Well, today's hardware is still all Von Neumann, just like in the 70s.
  • profquail 13y ago
    There's nothing requiring functional languages to abstract away from the bare metal. For example, read up on Typed Assembly Language (TAL) -- it's a typed, functional implementation of x86 assembly. I think the main reason functional languages are currently slower, in general, than procedural languages is that they haven't been popular for long enough yet to get all of the necessary optimizations implemented in the compilers. Think about it -- millions of person-hours have been invested into implementing optimizations in C and Fortran compilers over the past few decades. I expect the performance gap will get shrink steadily as time goes on, and as you said, at some point a cost/benefit threshold is reached where it makes more sense to switch to a functional language even if it's a little slower.
  • freyrs3 13y ago
    > because they are an abstraction away from the bare metal. This really is a myth, there's been plenty of work showing you can compile the lambda calculus to abstract machines which map nearly one-to-one on to assembly. High level functional languages make the same performance compromises that high level imperative languages make in performance ( garbage collectors, runtime dispatch, etc ), but there's nothing inherently about slower about compiling functional languages.
• mjn 13y ago
  The main takeaway here is that Haskell's stream fusion library is very nice. You can write code in the "naive" modular style that generates intermediate lists, and they get automagically fused away. There's a little work on that in Lispland as well. Some discussion, a prototype, and a link to the older SERIES package: http://pvk.ca/Blog/Lisp/Pipes/ http://pvk.ca/Blog/Lisp/Pipes/
  • DanWaterworth 13y ago
    Yes, this is an example where purity is really helpful.
• antimora 13y ago
  I have reimplemented the function in Julia and it took 16.51 seconds to run, compared to 0.19 seconds when running C version. Here is the code, perhaps someone else would know how to optimize it. CODE: function main(max_a0) longest = max_len = len = a = 0 for a0 in 1:max_a0 a = a0 len = 0 while a != 1 len += 1 a = ((a%2==0) ? a : 3*a+1)/2 end if len > max_len max_len = len longest = a0 end end println("($max_len, $longest)") end @time main(1000000) CMD: $ julia collatz.jl (329, 837799) elapsed time: 16.436623922 seconds (4206873264 bytes allocated) $ julia --version julia version 0.3.0-prerelease+143
  • antimora 13y ago
    WOW, I was able to optimize down to 0.88 seconds by specifying type for a variable like this a::Int64. It's by a factor of 18.6! Here is the new code: function main(max_a0) longest = max_len = len = a = 0 for a0 in 1:max_a0 a = a0 len = 0 while a != 1 len += 1 a::Int64 = ((a%2==0) ? a : 3*a+1)/2 end if len > max_len max_len = len longest = a0 end end println("($max_len, $longest)") end @time main(1000000) CMD: $ julia collatz.jl (329, 837799) elapsed time: 0.882312881 seconds (1377276 bytes allocated)
    • mhath 13y ago
      You can avoid allocating so much memory by changing a::Int64 = ((a%2==0) ? a : 3*a+1)/2 to a = ((a%2 == 0) ? a : 3a+1) >> 1 or a = div((a%2 == 0) ? a : 3a+1, 2) for me this gives (329, 837799) elapsed time: 0.449473338 seconds (544 bytes allocated) For reference the C version gives ./a.out 1000000 0.40s user 0.00s system 99% cpu 0.407 total
      • toolslive 13y ago
        what about: a = ((a & 0x01)? (a*3)+1 : a)>>1; anyway, when you think about it, you can skip everything between 0 and max_a0/2 because for any value there there is one with a longer chain (its double). that's a factor 2 in all languages.
• killercup 13y ago
  I was bored and curious to know how the Collatz conjecture works. Then, I implemented it in Lua and JS and compared their performance to C and Haskell. It's not meant to be a (good) benchmark. Here you go: https://gist.github.com/killercup/7718804 https://gist.github.com/killercup/7718804 tl;dr Implementation | Time -------------- | ----- clang 3.3 | 0.40s luajit 2.0.2 | 0.86s GHC 7.6.3 | 0.49s Node 0.10.22 | 2.43s
  • xdenser 13y ago
    my 2 cents to defend NodeJs 329 837799 naive 3.06 sec 329 837799 decompozed 3.03 sec 329 837799 memoized decompozed 0.52 sec The program https://gist.github.com/xdenser/7722083 https://gist.github.com/xdenser/7722083
    • dbaupp 13y ago
      Presumably the Lua, c and Haskell programs would also be 6× faster if you use the memoised algorithm, so that's a pretty meaningless comparison.
      • xdenser 13y ago
        Ok this intArithm variant (w/o memoization) for NodeJS is 6 times faster (0.54 sec) comparing to naive function intArithm(){ var max_a0 = +process.argv[2], longest = 0, max_len = 0, a0,len, a,z; for(a0=1;a0<=max_a0;++a0){ a = a0; len = 0; while( a != 1){ len++; z = a >>> 1; a = (!(a&1))?z:(a+z+1); } if(len>max_len){ max_len = len; longest = a0; } } console.log(max_len,longest); } the problem was with division by 2 operation which is not integer in js. So replacing it with shift >>> makes the trick. Formula optimization gives additional 200 ms. And I doubt memoized c or lua variant will be 6 times faster than non-memoized. As it is not that linear.
• nnq 13y ago
  Cython 0.47s less ...ignoring the cdefs, it's just as readable or more, at least for someone used to reading imperative code as opposed to functional code or mathematical equations. Plus that it takes you 5 minutes to quickly get from a quick and dirty Python hackoff to Cython. And you can stay blissful and ignorant about what "stream fusion" is, focusing on the problem you try to solve and not on the computer science of the tools. Haskell really is far far away from the "don't make me think" philosophy (or "don't make me think of anything else besides the real problem that I'm trying to solve"). No wonder Python and Cython are so loved in the scientific computing field :)
  • AndreasFrom 13y ago
    As you acknowledge yourself, whether Haskell makes you think or not relies on what you're trained to read.
    • nnq 13y ago
      Yes, you're right that Haskell is probably easier to read for someone with the mathematical background. But I'm referring to the cognitive overload incurred when you try to optimize it, thinking about list creation costs, then stream fusion to avoid it: most physicists or mathematicians will cringe at the thought of having to think about these "computer sciency optimization stuff". As opposed to this, a crude tool such as Cython or C or even Fortran allows one to just "skip thinking" about these and still have fast enough code. Haskell is great if you want to expand your CS knowledge, but not if your math/physics/engineering problem is already so complicated that it uses up all your available brain cycles and working memory and you just can't squeeze thinking about language-specific optimizations too. I guess this is why it's only used by companies like Galois and such that are already neck deep in advanced CS problems and used to thinking about them all the time.
      • PaulAJ 13y ago
        There is a saying: first make it work, then make it work fast. The problem with doing the second in imperative languages is that you need to pull all your working code apart and rewrite it. Haskell separates the "make it work" from the "make it work fast".
• lindig 13y ago
  Here is an implementation in OCaml. Time taken on a 1.3 GHz MacBook Air: :tmp $ time ./collatz.native 1000000 start: 837799 length 329 real 0m0.782s user 0m0.776s sys 0m0.004s exception Error of string let (@@) f x = f x let error fmt = Printf.kprintf (fun msg -> raise (Error msg)) fmt let max (a, x) (b, y) = if x > y then (a, x) else (b, y) let length start = let rec loop len a = let next = if a mod 2 = 0 then a/2 else (3*a+1)/2 in if a != 1 then loop (len+1) next else len in loop 0 start let collatz limit = let rec loop start longest = if start >= limit then longest else loop (start+1) (max longest (start, length start)) in loop 1 (1,0) let report (start, len) = Printf.printf "start: %d length %d\n" start len let main () = let argv = List.tl @@ Array.to_list Sys.argv in match argv with | [] -> report @@ collatz 10000 | [n] -> report @@ collatz @@ int_of_string n | _ -> error "usage: collatz [n]" let () = main ()
  • jds375 13y ago
    As a big OCaml fan, this is pretty cool. The speeds aren't that bad either. I hope OCaml gets more competitive in terms of running time in the future.
    • lindig 13y ago
      The code probably could be made faster by avoiding passing tuples (start, length) around as these are allocated on the heap. However, I found it a little bit more abstract.
• bcrack 13y ago
  It is interesting to point out that using -O3 optimization drops the C (and Cython) code exec time by roughly 45%.
• lelf 13y ago
  Here we go, Data.Vector and the same somewhat flawed comparison as in the post: next a = Just . join (,) $ (if even a then a else 3*a+1) `div` 2 len = (+1) . V.length . V.takeWhile (/=1) . V.unfoldr next result = V.maximum . V.map (\x -> (len x, x)) . V.enumFromTo 1 main = do [a0] <- getArgs let max_a0 = read a0 :: Word32 print . result $ max_a0 Results: haskell 5.0s, C 4.8s. Good enough I guess? PS: and this is how you do a proper benchmarking: import Criterion.Main main = defaultMain [ bench "data.vector/1M" $ whnf result 10000000 ] with result like http://lelf.lu/tmp/criterion.html http://lelf.lu/tmp/criterion.html
  • lelf 13y ago
    resultQ :: W -> (W,W) resultQ n = foldAllS max (0,0) $ runIdentity $ R.computeUnboxedP vec where fun (Z :. i) = len &&& id $ fromIntegral i vec = R.fromFunction (Z :. fromIntegral n) fun :: Array D DIM1 (W,W) Not as straight forward, but http://lelf.lu/files/fast-num-haskell-bench-N2.html http://lelf.lu/files/fast-num-haskell-bench-N2.html (now you now number of cores on my macbook I guess). Not bad for 3 lines of code?
• tsahyt 13y ago
  I've already posted this on /r/haskell on reddit but since the article turned up here as well, I might just as well post this again. Since we're talking optimization I had a short look at the C code and found it to be a very naive implementation and have managed to cut the execution time quite a bit by employing a sort of memoization. Once you get to a number of the sequence which you have already seen before, the length from there to 1 remains the same as before. So here is the code: http://lpaste.net/96397 http://lpaste.net/96397 Here are my measurements: time ./original 1000000 (329, 837799) real 0m0.183s user 0m0.180s sys 0m0.001s time ./new 1000000 (329, 837799) real 0m0.017s user 0m0.015s sys 0m0.001s
  • deletes 13y ago
    The idea is to optimize while keeping the algorithm the same throughout the languages, at least if you want any meaningful comparison. You are comparing languages not different algorithms. +1 for effort though.
    • anon1385 13y ago
      It doesn't seem that meaningful to compare a totally naive C implementation to an optimised Haskell one either. (Well, other than to show how fast C is even without any effort to optimise the code). Saying you have to keep the algorithm the same seems like a rather arbitrary limitation on the comparisons, especially since some algorithms are much less of a natural fit in some languages (see: quicksort in Haskell vs C [1]). I would suggest it's best to show (at least) 4 versions of the code to compare two language: the naive implementation in both languages, and the really optimised version in both languages (plus intermediate levels of optimisation if you feel dedicated). That gives a picture of how fast the typical code one would write will be, as well as the potential for optimisation if you are experiencing bottlenecks, and how much effort it takes to optimise the code. [1] http://augustss.blogspot.co.uk/2007/08/quicksort-in-haskell-quicksort-is.html http://augustss.blogspot.co.uk/2007/08/quicksort-in-haskell-...
      • lelf 13y ago
• ghshephard 13y ago
  The python equivalent: real 0m21.891s On machine in which the C Variant: real 0m0.570s But - I learned about Cython today, so great article. Interesting that the Python Version takes so much longer. import sys check=int(sys.argv[1]) lnum=1 lcount=1 for x in range(1,check): len=1 a=x while (a>1): len+=1 if a%2==0: a=a/2 else: a=(a*3+1)/2 if len > lcount: lcount=len lnum=x print lcount,lnum
• IsTom 13y ago
  I don't see how C and Cython versions are different. It's almost word-for-word C.
• wereHamster 13y ago
  Why are the stream fusion functions in a separate module? Why were the functions in prelude not replaced by these more efficient ones?
  • PhineasRex 13y ago
    Partly for historical reasons and partly because Prelude is really intended for beginners.
    • tel 13y ago
      Also for definitional reasons!
      • DanWaterworth 13y ago
        Could you elaborate?
        • cantankerous 13y ago
          He means that there is a generally understood semantic meaning for the functions as they exist in the Haskell Prelude. Replacing them under the hood with something potentially different makes the "meaning" of the basic functions not quite the same. Haskell has a history rooted in this kind of thing so it continues in the community now.
          • tel 13y ago
            Yup! Thanks
          • wereHamster 13y ago
            But they are not different. The functions in stream-fusion are a drop-in replacement for the ones in Data.List. To quote the package description: "To use, simply import Data.List.Stream in place of Data.List, and hide list functions from the Prelude." So to put simply: you have two functions. They are semantically equivalent, but one is faster than the other. Why not use the faster one?
• lelf 13y ago
  http://lelf.lu/files/fast-num-haskell-bench.html http://lelf.lu/files/fast-num-haskell-bench.html Without start-up overhead. C versions via FFI. https://gist.github.com/llelf/7721211 https://gist.github.com/llelf/7721211 (quick'n'shitty, sorry)
• glogla 13y ago
  I made quick Clojure version [1] it's really slow, and probably terrible, but I wonder if someone more knowledgeable in Clojure could look at it and optimize it the way it happened with the Haskell version [1]: http://pastebin.com/PH6DYbmr http://pastebin.com/PH6DYbmr
  • cbp 13y ago
    Since you know the division in collatz-next is always exact you can change / for quot or a bit-shift-right. That should improve it somewhat. You should also typehint with ^long the arguments to your functions. Further than that I would probably compile the code and then use a decompiler to check for boxing. But I can't say I've done much math intensive programming in clojure.
    • wging 13y ago
      Nice. 4.5 -> 2.7 seconds on my machine when using 'quot' instead of '/'. 2.4 seconds when hinting ^long in argument of collatz-length. 1.6 when also hinting ^long in collatz-next. Interestingly, when I hint only in collatz-next it takes 1.5.
    • glogla 13y ago
      Thanks! I was completely unaware that such type-hinting existed. I was also unsure if loop/recur was a good loop paradigm for this kind of problem (maybe for with variable would be better?) but I'd say in 1/3 speed of C it's pretty good for Clojure :)
• brihat 13y ago
  A naive implementation in Nimrod: proc main(max_a0: int): int = var a, longest, len, max_len : int for a0 in countup(1, max_a0): a = a0 len = 0 while a != 1: len += 1 if (a mod 2 != 0): a = (3*a + 1) a = a div 2 if len > max_len: max_len = len longest = a0 return longest # Main program starts here echo(main(1000000)) Takes 0.76s (the C program in TFA takes 0.58s)
• Nican 13y ago
  I just started learning haskell, can someone explain what: (\a0 -> (collatzLen a0, a0)) is doing?
  • wisty 13y ago
    > print $ maximum $ map (\a0 -> (collatzLen a0, a0)) [1..max_a0] I think the last bit of the line is: for a0 in range(1,max_a0): collatzLen(a0,a0) a0 is basically the loop variable.
    • evincarofautumn 13y ago
      Function application has high precedence in Haskell, so it would be (collatzLen(a0), a0).
  • chas 13y ago
    That is a function that takes a starting parameter and returns a pair with the length of the sequence as the first element and the starting value as the second element. This is mapped over a list of integers from 1 to max_a0 to find the length of the sequence starting with each integer. Maximum finds the largest component of a list and compares on the first element of a pair first, so the pair with the longest sequence length as its first element and the starting value that produced that sequence length as it's second element will be returned by `maximum $ map (\a0 -> (collatzLen a0, a0)) [1..max_a0]`. This function could also be expressed as `collatzLen >>= (,)`, if you were into that sort of thing.
  • tsahyt 13y ago
    chas explained it very well already, I'd just like to add that this is called a "lambda function" (or lambda expression) or an anonymous function, just in case you want to look it up. There's more on it here: http://www.haskell.org/haskellwiki/Anonymous_function http://www.haskell.org/haskellwiki/Anonymous_function
• carterschonwald 13y ago
  Ok. I'm going to cook up a variant that should match the c version in dirty dirty style Haskell. Will post it after lunch. The list fusion stuff is great for getting performance out of composing systems, but leaves those last dangling constant factors vs the ideal tight inner loop.
• craigyk 13y ago
  what is the preferred python speed-up route these days? I've tried: numpy C API cython weave numba parakeet
• iskander 13y ago
  Timings on my machine (2.67ghz Xeon W3520): $time ./c_compiled_with_O2 1000000 real 0m0.362s user 0m0.360s $time ./c_compiled_with_O3 1000000 real 0m0.189s user 0m0.184s I also tried dropping all the cdefs from the Cython code and running the Python version in Parakeet (http://www.parakeetpython.com http://www.parakeetpython.com) and it took 0.57s (whereas CPython without the @jit decorator took 19.57s).
• bjorg 13y ago
  Is it just me or does nobody care about the performance times being so low to begin with? Wouldn't it be better if the algorithm ran for 10+ seconds at least? Maybe closer to a minute. That would reduce variance due to operating system interference, as well as hide the process startup cost, which shouldn't be relevant here.
• rakoo 13y ago
  Interestingly, this naive go version is faster than the naive C version: package main import ( "log" ) func main() { max := 1000000 var maxa int var maxlength int for a0 := 0; a0 < max; a0++ { var length int a := a0 for a > 1 { if a%2 == 0 { a = a/2 } else { a = (3 * a + 1) / 2 } length++ } if length > maxlength { maxa = a0 maxlength = length } } log.Println(maxlength, maxa) } Compiled with go 1.1: > go version go version go1.1.2 linux/amd64 Result: > time ./collatz 2013/12/01 13:51:31 329 837799 ./collatz 0,34s user 0,00s system 99% cpu 0,344 total
  • rakoo 13y ago
    Oh, and I just tested this on play.golang.org, and I had wrong results because a was too small, so we need to manually set it to an uint64: package main import ( "log" ) func main() { max := 1000000 var maxa int var maxlength int for a0 := 0; a0 < max; a0++ { var length int a := uint64(a0) for a > 1 { if a%2 == 0 { a = a/2 } else { a = (3 * a + 1) / 2 } length++ } if length > maxlength { maxa = a0 maxlength = length } } log.Println(maxlength, maxa) } Even more interestingly, this is even faster: > time ./collatz 2013/12/01 13:59:06 329 837799 ./collatz 0,27s user 0,00s system 99% cpu 0,276 total