diff --git a/spot/kripke/kripke.hh b/spot/kripke/kripke.hh
index 0784a6740..393a035e9 100644
--- a/spot/kripke/kripke.hh
+++ b/spot/kripke/kripke.hh
@@ -38,18 +38,19 @@ namespace spot
     public std::enable_shared_from_this<kripkecube<State, SuccIterator>>
   {
   public:
-    /// \brief Returns the initial State of the System
-    State initial();
+    /// \brief Returns the initial State of the System. The \a tid parameter
+    /// is used internally for sharing this structure among threads.
+    State initial(unsigned tid);
 
     /// \brief Provides a string representation of the parameter state
-    std::string to_string(const State) const;
+    std::string to_string(const State, unsigned tid) const;
 
     /// \brief Returns an iterator over the successors of the parameter state.
-    SuccIterator* succ(const State);
+    SuccIterator* succ(const State, unsigned tid);
 
     /// \brief Allocation and deallocation of iterator is costly. This
     /// method allows to reuse old iterators.
-    void recycle(SuccIterator*);
+    void recycle(SuccIterator*, unsigned tid);
 
     /// \brief This method allow to deallocate a given state.
     const std::vector<std::string> get_ap();
@@ -62,11 +63,11 @@ namespace spot
   bool is_a_kripkecube(kripkecube<State, SuccIter>&)
   {
     // Hardly waiting C++17 concepts...
-    State (kripkecube<State, SuccIter>::*test_initial)() =
+    State (kripkecube<State, SuccIter>::*test_initial)(unsigned) =
       &kripkecube<State, SuccIter>::initial;
     std::string (kripkecube<State, SuccIter>::*test_to_string)
-      (const State) const = &kripkecube<State, SuccIter>::to_string;
-    auto (kripkecube<State, SuccIter>::*test_recycle)(SuccIter*) =
+      (const State, unsigned) const = &kripkecube<State, SuccIter>::to_string;
+    auto (kripkecube<State, SuccIter>::*test_recycle)(SuccIter*, unsigned) =
       &kripkecube<State, SuccIter>::recycle;
     const std::vector<std::string>
       (kripkecube<State, SuccIter>::*test_get_ap)() =
diff --git a/spot/ltsmin/ltsmin.cc b/spot/ltsmin/ltsmin.cc
index 894323873..f439422ca 100644
--- a/spot/ltsmin/ltsmin.cc
+++ b/spot/ltsmin/ltsmin.cc
@@ -1260,8 +1260,9 @@ namespace spot
                     cube cond,
                     bool compress,
                     bool selfloopize,
-                    cubeset& cubeset, int dead_idx)
-      :  current_(0), cond_(cond)
+                    cubeset& cubeset,
+                    int dead_idx, unsigned tid)
+      :  current_(0), cond_(cond), tid_(tid)
     {
       successors_.reserve(10);
       inner.manager = &manager;
@@ -1307,8 +1308,9 @@ namespace spot
                  cube cond,
                  bool compress,
                  bool selfloopize,
-                 cubeset& cubeset, int dead_idx)
+                 cubeset& cubeset, int dead_idx, unsigned tid)
     {
+      tid_ = tid;
       cond_ = cond;
       current_ = 0;
       // Constant time since int* is is_trivially_destructible
@@ -1365,7 +1367,10 @@ namespace spot
 
     cspins_state state() const
     {
-      return successors_[current_];
+      if (SPOT_UNLIKELY(!tid_))
+        return successors_[current_];
+      return  successors_[(((current_+1)*primes[tid_])
+                           % ((int)successors_.size()))];
     }
 
     cube condition() const
@@ -1377,6 +1382,7 @@ namespace spot
     std::vector<cspins_state> successors_;
     unsigned int current_;
     cube cond_;
+    unsigned tid_;
   };
 
   ////////////////////////////////////////////////////////////////////////
@@ -1409,45 +1415,68 @@ namespace spot
   public:
     kripkecube(spins_interface_ptr sip, bool compress,
                std::vector<std::string> visible_aps,
-               bool selfloopize, std::string dead_prop)
-      : sip_(sip), d_(sip.get()), manager_(d_->get_state_size(), compress),
-      compress_(compress), cubeset_(visible_aps.size()),
-      selfloopize_(selfloopize), aps_(visible_aps)
+               bool selfloopize, std::string dead_prop,
+               unsigned int nb_threads)
+      : sip_(sip), d_(sip.get()),
+        compress_(compress), cubeset_(visible_aps.size()),
+        selfloopize_(selfloopize), aps_(visible_aps), nb_threads_(nb_threads)
     {
       map_.setSize(2000000);
-      recycle_.reserve(2000000);
-      inner_.compressed_ = new int[d_->get_state_size() * 2];
-      inner_.uncompressed_ = new int[d_->get_state_size()+30];
+      manager_ = static_cast<cspins_state_manager*>
+        (::operator new(sizeof(cspins_state_manager) * nb_threads));
+      inner_ = new inner_callback_parameters[nb_threads_];
+      for (unsigned i = 0; i < nb_threads_; ++i)
+        {
+          recycle_.push_back(std::vector<cspins_iterator*>());
+          recycle_.back().reserve(2000000);
+          new (&manager_[i])
+            cspins_state_manager(d_->get_state_size(), compress);
+          inner_[i].compressed_ = new int[d_->get_state_size() * 2];
+          inner_[i].uncompressed_ = new int[d_->get_state_size()+30];
+        }
       dead_idx_ = -1;
       match_aps(visible_aps, dead_prop);
+
     }
+
     ~kripkecube()
       {
         for (auto i: recycle_)
           {
-            cubeset_.release(i->condition());
-            delete i;
+            for (auto j: i)
+              {
+                cubeset_.release(j->condition());
+                delete j;
+              }
           }
-        delete inner_.compressed_;
-        delete inner_.uncompressed_;
+
+        for (unsigned i = 0; i < nb_threads_; ++i)
+          {
+            manager_[i].~cspins_state_manager();
+            delete inner_[i].compressed_;
+            delete inner_[i].uncompressed_;
+          }
+        delete[] inner_;
       }
 
-    cspins_state initial()
+    cspins_state initial(unsigned tid)
     {
-      d_->get_initial_state(inner_.uncompressed_);
-      return manager_.alloc_setup(inner_.uncompressed_, inner_.compressed_,
-                                  manager_.size() * 2);
+      d_->get_initial_state(inner_[tid].uncompressed_);
+      return manager_[tid].alloc_setup(inner_[tid].uncompressed_,
+                                  inner_[tid].compressed_,
+                                  manager_[tid].size() * 2);
     }
 
-    std::string to_string(const cspins_state s) const
+    std::string to_string(const cspins_state s, unsigned tid = 0) const
     {
       std::string res = "";
-      cspins_state out = manager_.unbox_state(s);
+      cspins_state out = manager_[tid].unbox_state(s);
       cspins_state ref = out;
       if (compress_)
         {
-          manager_.decompress(s, inner_.uncompressed_, manager_.size());
-          ref = inner_.uncompressed_;
+          manager_[tid].decompress(s, inner_[tid].uncompressed_,
+                                   manager_[tid].size());
+          ref = inner_[tid].uncompressed_;
         }
       for (int i = 0; i < d_->get_state_size(); ++i)
         res += (d_->get_state_variable_name(i)) +
@@ -1456,28 +1485,28 @@ namespace spot
       return res;
     }
 
-    cspins_iterator* succ(const cspins_state s)
+    cspins_iterator* succ(const cspins_state s, unsigned tid)
     {
-      if (SPOT_LIKELY(!recycle_.empty()))
+      if (SPOT_LIKELY(!recycle_[tid].empty()))
         {
-          auto tmp  = recycle_.back();
-          recycle_.pop_back();
-          compute_condition(tmp->condition(), s);
-          tmp->recycle(s, d_, manager_, map_, inner_, -1,
+          auto tmp  = recycle_[tid].back();
+          recycle_[tid].pop_back();
+          compute_condition(tmp->condition(), s, tid);
+          tmp->recycle(s, d_, manager_[tid], map_, inner_[tid], -1,
                        tmp->condition(), compress_, selfloopize_,
-                       cubeset_, dead_idx_);
+                       cubeset_, dead_idx_, tid);
           return tmp;
         }
       cube cond = cubeset_.alloc();
-      compute_condition(cond, s);
-      return new cspins_iterator(s, d_, manager_, map_, inner_,
+      compute_condition(cond, s, tid);
+      return new cspins_iterator(s, d_, manager_[tid], map_, inner_[tid],
                                  -1, cond, compress_, selfloopize_,
-                                 cubeset_, dead_idx_);
+                                 cubeset_, dead_idx_, tid);
     }
 
-    void recycle(cspins_iterator* it)
+    void recycle(cspins_iterator* it, unsigned tid)
     {
-      recycle_.push_back(it);
+      recycle_[tid].push_back(it);
     }
 
     const std::vector<std::string> get_ap()
@@ -1485,6 +1514,11 @@ namespace spot
       return aps_;
     }
 
+    unsigned get_threads()
+    {
+      return nb_threads_;
+    }
+
   private:
     // The two followings functions are too big  to be inlined in
     // this class. See below for more details
@@ -1495,34 +1529,36 @@ namespace spot
 
     /// Compute the cube associated to each state. The cube
     /// will then be given to all iterators.
-    void compute_condition(cube c, cspins_state s);
+    void compute_condition(cube c, cspins_state s, unsigned tid = 0);
 
 
     spins_interface_ptr sip_;
     const spot::spins_interface* d_; // To avoid numerous sip_.get()
-    cspins_state_manager manager_; // FIXME One per thread!
+    cspins_state_manager* manager_;
     bool compress_;
-    std::vector<cspins_iterator*> recycle_;
-    inner_callback_parameters inner_; // FIXME Should be an array for threads.
+    std::vector<std::vector<cspins_iterator*>> recycle_;
+    inner_callback_parameters* inner_;
     cubeset cubeset_;
     bool selfloopize_;
     int dead_idx_;
     std::vector<std::string> aps_;
     cspins_state_map map_;
+    unsigned int nb_threads_;
   };
 
   void
-  kripkecube<cspins_state,
-             cspins_iterator>::compute_condition(cube c, cspins_state s)
+  kripkecube<cspins_state, cspins_iterator>::compute_condition
+  (cube c, cspins_state s, unsigned tid)
   {
     int i = -1;
-    int *vars = manager_.unbox_state(s);
+    int *vars = manager_[tid].unbox_state(s);
 
     // State is compressed, uncompress it
     if (compress_)
       {
-        manager_.decompress(s, inner_.uncompressed_+2, manager_.size());
-        vars = inner_.uncompressed_ + 2;
+        manager_[tid].decompress(s, inner_[tid].uncompressed_+2,
+                                 manager_[tid].size());
+        vars = inner_[tid].uncompressed_ + 2;
       }
 
     for (auto& ap: pset_)
@@ -1943,7 +1979,8 @@ namespace spot
 
   ltsmin_kripkecube_ptr
   ltsmin_model::kripkecube(std::vector<std::string> to_observe,
-                           const formula dead, int compress) const
+                           const formula dead, int compress,
+                           unsigned int nb_threads) const
   {
     // Register the "dead" proposition.  There are three cases to
     // consider:
@@ -1973,7 +2010,7 @@ namespace spot
     // Finally build the system.
     return std::make_shared<spot::kripkecube<spot::cspins_state,
                                              spot::cspins_iterator>>
-      (iface, compress, to_observe, selfloopize, dead_ap);
+      (iface, compress, to_observe, selfloopize, dead_ap, nb_threads);
   }
 
   kripke_ptr
@@ -2008,7 +2045,7 @@ namespace spot
   {
   }
 
-  std::tuple<bool, std::string, istats>
+  std::tuple<bool, std::string, std::vector<istats>>
   ltsmin_model::modelcheck(ltsmin_kripkecube_ptr sys,
                            spot::twacube_ptr twa, bool compute_ctrx)
   {
@@ -2018,13 +2055,33 @@ namespace spot
     for (unsigned int i = 0; i < sys->get_ap().size(); ++i)
       assert(sys->get_ap()[i].compare(twa->get_ap()[i]) == 0);
 
-    ec_renault13lpar<cspins_state, cspins_iterator,
-                     cspins_state_hash, cspins_state_equal> ec(*sys, twa);
-    bool has_ctrx = ec.run();
+    bool stop = false;
+    std::vector<ec_renault13lpar<cspins_state, cspins_iterator,
+                                 cspins_state_hash, cspins_state_equal>> ecs;
+    for (unsigned i = 0; i < sys->get_threads(); ++i)
+      ecs.push_back({*sys, twa, i, stop});
+
+    std::vector<std::thread> threads;
+    for (unsigned i = 0; i < sys->get_threads(); ++i)
+      threads.push_back
+        (std::thread(&ec_renault13lpar<cspins_state, cspins_iterator,
+                     cspins_state_hash, cspins_state_equal>::run, &ecs[i]));
+
+    for (unsigned i = 0; i < sys->get_threads(); ++i)
+      threads[i].join();
+
+    bool has_ctrx = false;
     std::string trace = "";
-    if (has_ctrx && compute_ctrx)
-      trace = ec.trace();
-    return std::make_tuple(has_ctrx, trace, ec.stats());
+    std::vector<istats> stats;
+    for (unsigned i = 0; i < sys->get_threads(); ++i)
+      {
+        has_ctrx |= ecs[i].counterexample_found();
+        if (compute_ctrx && ecs[i].counterexample_found()
+            && trace.compare("") == 0)
+          trace = ecs[i].trace(); // Pick randomly one !
+        stats.push_back(ecs[i].stats());
+      }
+    return std::make_tuple(has_ctrx, trace, stats);
   }
 
   int ltsmin_model::state_size() const
diff --git a/spot/ltsmin/ltsmin.hh b/spot/ltsmin/ltsmin.hh
index 38cb3a2d6..da787c631 100644
--- a/spot/ltsmin/ltsmin.hh
+++ b/spot/ltsmin/ltsmin.hh
@@ -85,13 +85,14 @@ namespace spot
     // atomic propositions such as "P.a == P.c"
     ltsmin_kripkecube_ptr kripkecube(std::vector<std::string> to_observe,
                                      formula dead = formula::tt(),
-                                     int compress = 0) const;
+                                     int compress = 0,
+                                     unsigned int nb_threads = 1) const;
 
     /// \brief Check for the emptiness between a system and a twa.
     /// Return a pair containing a boolean indicating wether a counterexample
     /// has been found and a string representing the counterexample if the
     /// computation have been required
-    static std::tuple<bool, std::string, istats>
+    static std::tuple<bool, std::string, std::vector<istats>>
     modelcheck(ltsmin_kripkecube_ptr sys,
                spot::twacube_ptr twa, bool compute_ctrx = false);
 
diff --git a/spot/mc/ec.hh b/spot/mc/ec.hh
index 86f799ceb..e14ddc149 100644
--- a/spot/mc/ec.hh
+++ b/spot/mc/ec.hh
@@ -45,10 +45,10 @@ namespace spot
 
   public:
     ec_renault13lpar(kripkecube<State, SuccIterator>& sys,
-                     twacube_ptr twa)
+                     twacube_ptr twa, unsigned tid, bool stop)
       : intersect<State, SuccIterator, StateHash, StateEqual,
                   ec_renault13lpar<State, SuccIterator,
-                                   StateHash, StateEqual>>(sys, twa),
+                                   StateHash, StateEqual>>(sys, twa, tid, stop),
       acc_(twa->acc()), sccs_(0U)
       {
       }
@@ -114,6 +114,8 @@ namespace spot
         }
       roots_.back().acc |= cond;
       found_ = acc_.accepting(roots_.back().acc);
+      if (SPOT_UNLIKELY(found_))
+        this->stop_ = true;
       return found_;
     }
 
@@ -147,7 +149,7 @@ namespace spot
       acc_cond::mark_t acc = 0U;
 
       bfs.push(new ctrx_element({&this->todo.back().st, nullptr,
-              this->sys_.succ(this->todo.back().st.st_kripke),
+              this->sys_.succ(this->todo.back().st.st_kripke, this->tid_),
               this->twa_->succ(this->todo.back().st.st_prop)}));
       while (true)
         {
@@ -160,7 +162,7 @@ namespace spot
               while (!front->it_prop->done())
                 {
                   if (this->twa_->get_cubeset().intersect
-                      (this->twa_->trans_data(front->it_prop).cube_,
+                      (this->twa_->trans_data(front->it_prop, this->tid_).cube_,
                        front->it_kripke->condition()))
                     {
                       const product_state dst = {
@@ -181,7 +183,8 @@ namespace spot
                       // This is a valid transition. If this transition
                       // is the one we are looking for, update the counter-
                       // -example and flush the bfs queue.
-                      auto mark = this->twa_->trans_data(front->it_prop).acc_;
+                      auto mark = this->twa_->trans_data(front->it_prop,
+                                                         this->tid_).acc_;
                       if (!acc.has(mark))
                         {
                           ctrx_element* current = front;
@@ -213,7 +216,7 @@ namespace spot
                           const product_state* q = &(it->first);
                           ctrx_element* root = new ctrx_element({
                               q , nullptr,
-                              this->sys_.succ(q->st_kripke),
+                              this->sys_.succ(q->st_kripke, this->tid_),
                               this->twa_->succ(q->st_prop)
                           });
                           bfs.push(root);
@@ -224,7 +227,7 @@ namespace spot
                       const product_state* q = &(it->first);
                       ctrx_element* root = new ctrx_element({
                           q , nullptr,
-                          this->sys_.succ(q->st_kripke),
+                          this->sys_.succ(q->st_kripke, this->tid_),
                           this->twa_->succ(q->st_prop)
                       });
                       bfs.push(root);
@@ -243,7 +246,7 @@ namespace spot
     virtual istats stats() override
     {
       return {this->states(), this->trans(), sccs_,
-              (unsigned)roots_.size(), dfs_};
+          (unsigned) roots_.size(), dfs_, found_};
     }
 
   private:
diff --git a/spot/mc/intersect.hh b/spot/mc/intersect.hh
index b13187a97..b47853812 100644
--- a/spot/mc/intersect.hh
+++ b/spot/mc/intersect.hh
@@ -34,6 +34,7 @@ namespace spot
     unsigned sccs;
     unsigned instack_sccs;
     unsigned instack_item;
+    bool is_empty;
   };
 
   /// \brief This class explores (with a DFS) a product between a
@@ -59,8 +60,8 @@ namespace spot
   {
   public:
     intersect(kripkecube<State, SuccIterator>& sys,
-              twacube_ptr twa):
-      sys_(sys), twa_(twa)
+              twacube_ptr twa, unsigned tid, bool& stop):
+      sys_(sys), twa_(twa), tid_(tid), stop_(stop)
     {
       SPOT_ASSERT(is_a_kripkecube(sys));
       map.reserve(2000000);
@@ -87,10 +88,10 @@ namespace spot
     bool run()
     {
       self().setup();
-      product_state initial = {sys_.initial(), twa_->get_initial()};
+      product_state initial = {sys_.initial(tid_), twa_->get_initial()};
       if (SPOT_LIKELY(self().push_state(initial, dfs_number+1, 0U)))
         {
-          todo.push_back({initial, sys_.succ(initial.st_kripke),
+          todo.push_back({initial, sys_.succ(initial.st_kripke, tid_),
                 twa_->succ(initial.st_prop)});
 
           // Not going further! It's a product with a single state.
@@ -100,7 +101,7 @@ namespace spot
           forward_iterators(true);
           map[initial] = ++dfs_number;
         }
-      while (!todo.empty())
+      while (!todo.empty() && !stop_)
         {
           // Check the kripke is enough since it's the outer loop. More
           // details in forward_iterators.
@@ -113,8 +114,8 @@ namespace spot
                                              is_init,
                                              newtop,
                                              map[newtop])))
-              {
-                  sys_.recycle(todo.back().it_kripke);
+                {
+                  sys_.recycle(todo.back().it_kripke, tid_);
                   // FIXME a local storage for twacube iterator?
                   todo.pop_back();
                   if (SPOT_UNLIKELY(self().counterexample_found()))
@@ -126,9 +127,9 @@ namespace spot
               ++transitions;
               product_state dst = {
                 todo.back().it_kripke->state(),
-                twa_->trans_storage(todo.back().it_prop).dst
+                twa_->trans_storage(todo.back().it_prop, tid_).dst
               };
-              auto acc = twa_->trans_data(todo.back().it_prop).acc_;
+              auto acc = twa_->trans_data(todo.back().it_prop, tid_).acc_;
               forward_iterators(false);
               auto it  = map.find(dst);
               if (it == map.end())
@@ -136,7 +137,7 @@ namespace spot
                   if (SPOT_LIKELY(self().push_state(dst, dfs_number+1, acc)))
                     {
                       map[dst] = ++dfs_number;
-                      todo.push_back({dst, sys_.succ(dst.st_kripke),
+                      todo.push_back({dst, sys_.succ(dst.st_kripke, tid_),
                             twa_->succ(dst.st_prop)});
                       forward_iterators(true);
                     }
@@ -167,7 +168,7 @@ namespace spot
 
     virtual istats stats()
     {
-      return {dfs_number, transitions, 0U, 0U, 0U};
+      return {dfs_number, transitions, 0U, 0U, 0U, false};
     }
 
   protected:
@@ -190,7 +191,7 @@ namespace spot
       // There is no need to move iterators forward.
       SPOT_ASSERT(!(todo.back().it_prop->done()));
       if (just_pushed && twa_->get_cubeset()
-          .intersect(twa_->trans_data(todo.back().it_prop).cube_,
+          .intersect(twa_->trans_data(todo.back().it_prop, tid_).cube_,
                      todo.back().it_kripke->condition()))
           return;
 
@@ -207,7 +208,7 @@ namespace spot
           while (!todo.back().it_prop->done())
             {
               if (SPOT_UNLIKELY(twa_->get_cubeset()
-                .intersect(twa_->trans_data(todo.back().it_prop).cube_,
+                .intersect(twa_->trans_data(todo.back().it_prop, tid_).cube_,
                            todo.back().it_kripke->condition())))
                 return;
               todo.back().it_prop->next();
@@ -263,5 +264,7 @@ namespace spot
     visited_map map;
     unsigned int dfs_number = 0;
     unsigned int transitions = 0;
+    unsigned tid_;
+    bool& stop_; // Do not need to be atomic.
   };
 }
diff --git a/spot/mc/reachability.hh b/spot/mc/reachability.hh
index 3c52a75f2..abd43517c 100644
--- a/spot/mc/reachability.hh
+++ b/spot/mc/reachability.hh
@@ -32,8 +32,8 @@ namespace spot
   class SPOT_API seq_reach_kripke
   {
   public:
-    seq_reach_kripke(kripkecube<State, SuccIterator>& sys):
-      sys_(sys)
+    seq_reach_kripke(kripkecube<State, SuccIterator>& sys, unsigned tid):
+      sys_(sys), tid_(tid)
     {
       SPOT_ASSERT(is_a_kripkecube(sys));
       visited.reserve(2000000);
@@ -54,15 +54,15 @@ namespace spot
     void run()
     {
       self().setup();
-      State initial = sys_.initial();
-      todo.push_back({initial, sys_.succ(initial)});
+      State initial = sys_.initial(tid_);
+      todo.push_back({initial, sys_.succ(initial, tid_)});
       visited[initial] = ++dfs_number;
       self().push(initial, dfs_number);
       while (!todo.empty())
         {
           if (todo.back().it->done())
             {
-              sys_.recycle(todo.back().it);
+              sys_.recycle(todo.back().it, tid_);
               todo.pop_back();
             }
           else
@@ -76,7 +76,7 @@ namespace spot
                   self().push(dst, dfs_number);
                   self().edge(visited[todo.back().s], dfs_number);
                   todo.back().it->next();
-                  todo.push_back({dst, sys_.succ(dst)});
+                  todo.push_back({dst, sys_.succ(dst, tid_)});
                 }
               else
                 {
@@ -113,5 +113,6 @@ namespace spot
     visited_map visited;
     unsigned int dfs_number = 0;
     unsigned int transitions = 0;
+    unsigned int tid_;
   };
 }
diff --git a/spot/misc/hash.hh b/spot/misc/hash.hh
index f0f55eea1..cc5e93ef5 100644
--- a/spot/misc/hash.hh
+++ b/spot/misc/hash.hh
@@ -87,4 +87,34 @@ namespace spot
                          static_cast<size_t>(uh(p.second)));
     }
   };
+
+  // From primes.utm.edu shuffled. This primes are used when we simulate
+  // transition shuffling using "primitive root modulo n" technique.
+  static const unsigned primes[144] =
+    {
+      295075531, 334214857, 314607103, 314607191, 334214891, 334214779,
+      295075421, 472882073, 256203329, 275604599, 314606953, 256203397,
+      275604547, 256203631, 275604617, 472882141, 472882297, 472882219,
+      256203229, 256203469, 275604643, 472882169, 275604803, 472882283,
+      295075463, 334214593, 295075213, 256203373, 314607019, 314607193,
+      295075399, 256203523, 314607001, 295075289, 256203293, 256203641,
+      256203307, 314607047, 295075373, 314607053, 314606977, 334214681,
+      275604691, 275604577, 472882447, 314607031, 275605019, 472882477,
+      472882499, 314607173, 295075241, 295075471, 295075367, 256203559,
+      314607233, 275604881, 334214941, 472882103, 275604821, 472882511,
+      295075357, 472882517, 314607023, 256203317, 295075337, 275605007,
+      472882391, 256203223, 334214723, 295075381, 295075423, 275604733,
+      314607113, 256203257, 472882453, 256203359, 295075283, 314607043,
+      256203403, 472882259, 314606891, 472882253, 314606917, 256203349,
+      256203457, 295075457, 472882171, 314607229, 295075513, 472882429,
+      334214953, 275604841, 295075309, 472882099, 334214467, 334214939,
+      275604869, 314607077, 314607089, 275604947, 275605027, 295075379,
+      334214861, 314606909, 334214911, 314607199, 275604983, 314606969,
+      256203221, 334214899, 256203611, 256203679, 472882337, 275604767,
+      472882199, 295075523, 472882049, 275604817, 334214561, 334214581,
+      334214663, 295075489, 295075163, 334214869, 334214521, 295075499,
+      275604913, 334214753, 334214687, 256203491, 295075153, 334214893,
+      472882411, 472882117, 275604793, 334214833, 334214591, 314607091,
+      256203419, 275604727, 256203659, 275604961, 334214557, 275604677
+    };
 }
diff --git a/spot/twacube/twacube.hh b/spot/twacube/twacube.hh
index 50a7d22b6..b8f91b5e7 100644
--- a/spot/twacube/twacube.hh
+++ b/spot/twacube/twacube.hh
@@ -22,6 +22,7 @@
 #include <vector>
 #include <iosfwd>
 #include <spot/graph/graph.hh>
+#include <spot/misc/hash.hh>
 #include <spot/twa/acc.hh>
 #include <spot/twacube/cube.hh>
 #include <spot/twacube/fwd.hh>
@@ -93,8 +94,14 @@ namespace spot
       // swarming, we expect swarming is activated.
       if (SPOT_UNLIKELY(!seed))
         return idx_;
-      // swarming.
-      return (((idx_-st_.succ+1)*seed) % (st_.succ_tail-st_.succ+1)) + st_.succ;
+      // Here swarming performs a technique called "primitive
+      // root modulo n", i. e.  for i in [1..n]: i*seed (mod n). We
+      // also must have seed prime with n: to solve this, we use
+      // precomputed primes and seed access one of this primes. Note
+      // that the chosen prime must be greater than n.
+      SPOT_ASSERT(primes[seed] > (st_.succ_tail-st_.succ+1));
+      return (((idx_-st_.succ+1)*primes[seed]) % (st_.succ_tail-st_.succ+1))
+        + st_.succ;
     }
 
   private:
diff --git a/tests/core/twacube.test b/tests/core/twacube.test
index 69f82f478..be3f4e213 100755
--- a/tests/core/twacube.test
+++ b/tests/core/twacube.test
@@ -1,6 +1,6 @@
 #!/bin/sh
 # -*- coding: utf-8 -*-
-# Copyright (C) 2014, 2015 Laboratoire de Recherche et Développement de
+# Copyright (C) 2014, 2015, 2016 Laboratoire de Recherche et Développement de
 # l'Epita (LRDE).
 #
 # This file is part of Spot, a model checking library.
@@ -61,10 +61,10 @@ init : 0
 2->1 : p1&p2 {}
 2->2 : 1 {0,1}
 -----------
-2 1  !p1 {}
 2 2  1 {0,1}
-2 0  p1 {0,1}
 2 1  p1&p2 {}
+2 1  !p1 {}
+2 0  p1 {0,1}
 2 0  !p1&p2 {0,1}
 digraph G {
   rankdir=LR
diff --git a/tests/ltsmin/modelcheck.cc b/tests/ltsmin/modelcheck.cc
index 6d0a65e3a..c164dc2c9 100644
--- a/tests/ltsmin/modelcheck.cc
+++ b/tests/ltsmin/modelcheck.cc
@@ -410,7 +410,8 @@ static int checked_main()
       try
         {
            modelcube = spot::ltsmin_model::load(mc_options.model)
-            .kripkecube(propcube->get_ap(), deadf, mc_options.compress);
+             .kripkecube(propcube->get_ap(), deadf, mc_options.compress,
+                         mc_options.nb_threads);
         }
       catch (std::runtime_error& e)
         {
@@ -432,32 +433,57 @@ static int checked_main()
         }
 
       // Display statistics
-      std::cout << std::get<2>(res).states << " unique states visited\n";
-      std::cout << std::get<2>(res).instack_sccs
-                << " strongly connected components in search stack\n";
-      std::cout << std::get<2>(res).transitions
-                << " transitions explored\n";
-      std::cout << std::get<2>(res).instack_item
-                << " items max in DFS search stack\n";
-      std::cout << memused << " pages allocated for emptiness check\n";
-
-      if (!std::get<0>(res))
-        std::cout << "no accepting run found";
-      else if (!mc_options.compute_counterexample)
+      unsigned smallest = 0;
+      for (unsigned i = 0; i < std::get<2>(res).size(); ++i)
         {
-          std::cout << "an accepting run exists "
-                    << "(use -c to print it)" << std::endl;
-          exit_code = 1;
+          if (std::get<2>(res)[i].states < std::get<2>(res)[smallest].states)
+            smallest = i;
+
+          std::cout << "\n---- Thread number : " << i << '\n';
+          std::cout << std::get<2>(res)[i].states << " unique states visited\n";
+          std::cout << std::get<2>(res)[i].instack_sccs
+                    << " strongly connected components in search stack\n";
+          std::cout << std::get<2>(res)[i].transitions
+                    << " transitions explored\n";
+          std::cout << std::get<2>(res)[i].instack_item
+                    << " items max in DFS search stack\n";
+
+          // FIXME: produce walltime for each threads.
+          if (mc_options.csv)
+            {
+              std::cout << "Find following the csv: "
+                        << "thread_id,walltimems,type,"
+                        << "states,transitions\n";
+              std::cout << "@th_" << i << ','
+                        << tm.timer("emptiness check").walltime() << ','
+                        << (!std::get<2>(res)[i].is_empty ?
+                            "EMPTY," : "NONEMPTY,")
+                        << std::get<2>(res)[i].states << ','
+                        << std::get<2>(res)[i].transitions
+                        << std::endl;
+            }
         }
-      else
-        std::cout << "an accepting run exists!\n" << std::get<1>(res)
-                  << std::endl;
+
 
       if (mc_options.csv)
         {
-          std::cout << "\nFind following the csv: "
-                    << "model,formula,walltimems,memused,type,"
+          std::cout << "\nSummary :\n";
+          if (!std::get<0>(res))
+            std::cout << "no accepting run found\n";
+          else if (!mc_options.compute_counterexample)
+            {
+              std::cout << "an accepting run exists "
+                        << "(use -c to print it)" << std::endl;
+              exit_code = 1;
+            }
+          else
+            std::cout << "an accepting run exists!\n" << std::get<1>(res)
+                      << '\n';
+
+          std::cout << "Find following the csv: "
+                    << "model,formula,walltimems,memused,type"
                     << "states,transitions\n";
+
           std::cout << '#'
                     << split_filename(mc_options.model)
                     << ','
@@ -465,9 +491,9 @@ static int checked_main()
                     << tm.timer("emptiness check").walltime() << ','
                     << memused << ','
                     << (!std::get<0>(res) ? "EMPTY," : "NONEMPTY,")
-                    << std::get<2>(res).states << ','
-                    << std::get<2>(res).transitions
-                    << std::endl;
+                    << std::get<2>(res)[smallest].states << ','
+                    << std::get<2>(res)[smallest].transitions
+                    << '\n';
         }
     }