Fast bitarray in OCaml

Question

Yet another synthetic benchmark: Sieve of Eratosthenes

C++

#include <vector>
#include <cmath>

void find_primes(int n, std::vector<int>& out)
{
   std::vector<bool> is_prime(n + 1, true);
   int last = sqrt(n);
   for (int i = 2; i <= last; ++i)
   {
      if (is_prime[i])
      {
         for (int j = i * i; j <= n; j += i)
         {
            is_prime[j] = false;
         }
      }
   }

   for (unsigned i = 2; i < is_prime.size(); ++i)
   {
      if (is_prime[i])
      {
         out.push_back(i);
      }
   }
}

OCaml (using Jane Street's Core and Res libraries)

open Core.Std
module Bits = Res.Bits
module Vect = Res.Array

let find_primes n =
  let is_prime = Bits.make (n + 1) true in
  let last = float n |! sqrt |! Float.iround_exn ~dir:`Zero in
  for i = 2 to last do
    if not (Bits.get is_prime i) then () else begin
      let j = ref (i * i) in
      while !j <= n; do
        Bits.set is_prime !j false;
        j := !j + i;
      done;
    end;
  done;
  let ar = Vect.empty () in
  for i = 2 to n do
    if Bits.get is_prime i then Vect.add_one ar i else ()
  done;
  ar

I was surprised that OCaml version (native) is about 13 times slower than C++. I replaced Res.Bits with Core_extended.Bitarray, but it became ~18 times slower. Why it is so slow? Doesn't OCaml provide fast operations for bit manipulation? Is there any alternative fast implementation of bit arrays?

To be clear: I'm from C++ world and consider OCaml as a possible alternative for writing performance critical code. Actually, I'm a bit scary with such results.

EDIT:

Profiling results

Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total           
 time   seconds   seconds    calls  ms/call  ms/call  name    
 50.81      1.26     1.26                             camlRes__pos_1113
  9.72      1.50     0.24                             camlRes__unsafe_get_1117
  6.68      1.66     0.17                             camlRes__unsafe_set_1122
  6.28      1.82     0.16                             camlNopres_impl__set_1054
  6.07      1.97     0.15                             camlNopres_impl__get_1051
  5.47      2.10     0.14 47786824     0.00     0.00  caml_apply3
  3.64      2.19     0.09 22106943     0.00     0.00  caml_apply2
  2.43      2.25     0.06   817003     0.00     0.00  caml_oldify_one
  2.02      2.30     0.05        1    50.00   265.14  camlPrimes__find_primes_64139
  1.21      2.33     0.03                             camlRes__unsafe_get_1041
...

Answer 1

Did you try using simple datastructure first before jumping on the sophisticated ones?

On my machine, the following code is only 4x slower than you C++ version (note that I made the minimal changes to use an Array as the cache, and a list to accumulate results; you could use the array get/set syntactic sugar):

let find_primes n =
  let is_prime = Array.make (n + 1) true in
  let last = int_of_float (sqrt (float n)) in
  for i = 2 to last do
    if not (Array.get is_prime i) then () else begin
      let j = ref (i * i) in
      while !j <= n; do
        Array.set is_prime !j false;
        j := !j + i;
      done;
    end;
  done;
  let ar = ref [] in
  for i = 2 to n do
    if Array.get is_prime i then ar := i :: !ar else ()
  done;
  ar

(4x slower: it takes 4s to compute the 10_000_000 first primes, vs. 1s for g++ -O1 or -O2 on your code)

Realizing that the efficiency of your bitvector solution probably comes from the economic memory layout, I changed the code to use strings instead of arrays:

let find_primes n =
  let is_prime = String.make (n + 1) '0' in
  let last = int_of_float (sqrt (float n)) in
  for i = 2 to last do
    if not (String.get is_prime i = '0') then () else begin
      let j = ref (i * i) in
      while !j <= n; do
        String.set is_prime !j '1';
        j := !j + i;
      done;
    end;
  done;
  let ar = ref [] in
  for i = 2 to n do
    if String.get is_prime i = '0' then ar := i :: !ar else ()
  done;
  ar

This now takes only 2s, which makes it 2x slower than your C++ solution.

Answer 2

It seems Jeffrey Scofield is right. Such terrible performance degradation is due to div and mod operations.

I prototyped small Bitarray module

module Bitarray = struct
  type t = { len : int; buf : string }

  let create len x =
    let init = (if x = true then '\255' else '\000') in
    let buf = String.make (len / 8 + 1) init in
    { len = len; buf = buf }

  let get t i =
    let ch = int_of_char (t.buf.[i lsr 3]) in
    let mask = 1 lsl (i land 7) in
    (ch land mask) <> 0

  let set t i b =
    let index = i lsr 3 in
    let ch = int_of_char (t.buf.[index]) in
    let mask = 1 lsl (i land 7) in
    let new_ch = if b then (ch lor mask) else (ch land lnot mask) in
    t.buf.[index] <- char_of_int new_ch
end

It uses string as byte array (8 bits per char). Initially I used x / 8 and x mod 8 for bit extraction. It was 10x slower than C++ code. Then I replaced them with x lsr 3 and x land 7. Now, it is only 4x slower than C++.

Answer 3

It's not often useful to compare micro-benchmarks like this, but the basic conclusion is probably correct. This is a case where OCaml is at a distinct disadvantage. C++ can access a more or less ideal representation (vector of machine integers). OCaml can make a vector, but can't get at the machine integers directly. So OCaml has to use div and mod where C++ can use shift and mask.

I reproduced this test (using a different bit vector library) and found that appreciable time in OCaml was spent constructing the result, which isn't a bit array. So the test might not be measuring exactly what you think.

Update

I tried some quick tests packing 32 booleans into a 63-bit int. It does seem to make things go faster, but only a little bit. It's not a perfect test, but it suggests gasche is right that the non-power-of-2 effect is minor.

Answer 4

Please make sure that you install Core including the .cmx file (.cmxa is not enough!), otherwise cross-module inlining will not work. Your profile suggests that certain calls may not have been inlined, which would explain a dramatic loss of efficiency.

Sadly, the Oasis packaging tool, which a lot of OCaml projects use, currently has a bug that prevents it from installing the .cmx file. The Core package is also affected by this problem, probably irrespective of which package manager (Opam, Godi) you use.