diff --git a/benchmark/make.sh b/benchmark/make.sh
index 8368e2c9b107fbb3062eb4ad8a5c2705d4f6b884..c88d389373933b8cd3725ae8f41f18c243132832 100755
--- a/benchmark/make.sh
+++ b/benchmark/make.sh
@@ -641,6 +641,7 @@ run_reach () {
   run_postlogue
 }
 # more experiments: -seq Array,PSek_reach,PSek_get,ParrayNarytree,StackVectorNoFunctor,NumeriChunk,FmlibRbtrees \
+#     -seq ParrayNarytree,PdarrayNarytree \
 
 plot_reach () {
   result="$1"
@@ -705,7 +706,7 @@ run_set () {
   prun -prog $prog \
     -n 10$K,20$K,40$K,80$K,160$K,320$K,640$K,2$M,8$M \
     -k 100$K \
-    -seq Array,ESek,NumeriChunk,ParrayNarytree \
+    -seq Array,ESek,NumeriChunk,ParrayNarytree5,PdarrayNarytree8 \
     -runs 3 \
     -distance -1 \
     -output $result
@@ -718,15 +719,19 @@ run_set_debug () {
   run_prologue
   prun -prog $prog \
     -n 10$K,20$K,40$K,80$K,160$K,320$K \
-    -k 10$K \
-    -seq Array,ESek,NumeriChunk,ParrayNarytree \
+    -k 50$K \
+    -seq ParrayNarytree5 \
     -runs 1 \
     -distance -1 \
     -output $result
   run_postlogue
 }
+# Fmarray --> very slow
 
+# ,ParrayNarytree8,PdarrayNarytree5,PdarrayNarytree8
+# Array,NumeriChunk,ESek,PSek,ParrayNarytree,,PdarrayNarytree7,PdarrayNarytree10
 #    -seq Array,PSek,ESek,FmlibRbtrees \
+#    -seq PdarrayNarytree,PdarrayNarytree8,PdarrayNarytree7,PdarrayNarytree10 \
 
 plot_set () {
   result="$1"
diff --git a/benchmark/src/darray.ml b/benchmark/src/darray.ml
index bb1c9defc9fcad657fcdcfbbe3c7758fa2ae0b55..cae0e0af277bb5ce7eedc707d93148c167f1d0fe 100644
--- a/benchmark/src/darray.ml
+++ b/benchmark/src/darray.ml
@@ -38,8 +38,6 @@ and 'a desc =
 
 let cast = Obj.magic
 
-let repr_none : 'a array = [||]
-
 let index_repr = max_int
 
 let value_none : 'a = cast ()
@@ -114,19 +112,11 @@ let mk_repr (type a) (dist : int) (a : a array) : a t =
     value = cast value_none;
     name = fresh_name() }
 
-let mk_diff (dist : int) (index : int) (value : 'a) (parent : 'a t) : 'a t =
-  register {
-    dist = dist;
-    index = index; (* never flipped when created *)
-    base = cast parent;
-    value = value;
-    name = fresh_name() }
-
 (** Parameters *)
 
 let dist_max (t : 'a t) : int =
   (* assert (is_repr t); *)
-  Array.length (get_repr t) - 1 (* TODO: /2 *)
+  Array.length (get_repr t) / 4
   (* could be any constant *)
 
 
@@ -232,12 +222,20 @@ let iter (f : 'a -> unit) (t : 'a t) : unit =
   ensure_repr t;
   Array.iter f (get_repr t)
 
+let of_array_with_transfer (a : 'a array) : 'a t =
+  mk_repr 0 a
+
+let to_array_with_transfer (t : 'a t) : 'a array =
+  if is_repr t
+    then get_repr t
+    else get_fresh_array t
+
 
 (*------------------------------------------*)
 (*
 module Test = struct
 
-let _test () =
+let _test  =
   let pr = Printf.printf in
 
   let print_array a =
@@ -249,6 +247,7 @@ let _test () =
   let print_node t =
     pr "Node %d (dist %d) = " (t.name) (get_dist t);
     if is_repr t then begin
+      pr "Repr ";
       print_array (get_repr t)
     end else begin
       let p = get_parent t in
@@ -262,6 +261,10 @@ let _test () =
     pr "-------------------\n";
     List.iter print_node (List.rev !nodes) in
 
+  (*
+    for this test, use
+    let dist_max (t : 'a t) : int = Array.length (get_repr t) - 1
+
   let n = 4 in
   let t0 = init n (fun i -> i) in (* [0;1;2;3] *)
   assert (get t0 2 = 2);
@@ -282,6 +285,14 @@ let _test () =
   p();
   let t8 = set t0 2 0 in  (* [0;1;0;3] flip *)
   p();
+  *)
+
+  let n = 30 in
+  let r = ref (init n (fun i -> i)) in
+  for i = 0 to n-1 do
+    r := set !r i (-i);
+    p();
+  done
 
 
 end
@@ -291,4 +302,16 @@ end
 module type NullSig = sig val null : 'a end
 module NullImpl : NullSig = struct let null = cast () end
 let value_none : 'a = NullImpl.null
-*)
\ No newline at end of file
+*)
+
+
+(* DEPRECATED
+let mk_diff (dist : int) (index : int) (value : 'a) (parent : 'a t) : 'a t =
+  register {
+    dist = dist;
+    index = index; (* never flipped when created *)
+    base = cast parent;
+    value = value;
+    name = fresh_name() }
+
+*)
diff --git a/benchmark/src/darray.mli b/benchmark/src/darray.mli
index 82eb5a83cabd662eba0c27223014c13a97f7411b..e00bd07ada7288af0ecc5e4709d25d5bdf1227c0 100644
--- a/benchmark/src/darray.mli
+++ b/benchmark/src/darray.mli
@@ -20,4 +20,5 @@ val init: int -> (int -> 'a) -> 'a t
 val iter: ('a -> unit) -> 'a t -> unit
 
 
-
+val of_array_with_transfer : 'a array -> 'a t
+val to_array_with_transfer : 'a t -> 'a array
diff --git a/benchmark/src/numerichunk.ml b/benchmark/src/numerichunk.ml
index 30d0bb8464a090f153baffb2f902a279b8c4f71a..2629b08bb70b67a188c67d39b4d50ac18317b252 100644
--- a/benchmark/src/numerichunk.ml
+++ b/benchmark/src/numerichunk.ml
@@ -270,17 +270,17 @@ end
 
 (*--------------------------------------------------------------*)
 
-(*
+
 module BITS_PER_DIGIT_5 : BITS_PER_DIGIT = struct
   let bits_per_digit : int = 5
 end
 include Make(BITS_PER_DIGIT_5)
-*)
+(*
 module BITS_PER_DIGIT_8 : BITS_PER_DIGIT = struct
   let bits_per_digit : int = 8
 end
 include Make(BITS_PER_DIGIT_8)
-
+*)
 
 (*--------------------------------------------------------------
 
diff --git a/benchmark/src/parray_narytree.ml b/benchmark/src/parray_narytree.ml
index d97cd47960ac10a177d8d0715c75ff5645b73478..2ac2833b1a7a0328b2f019c34715f2031a793a4f 100644
--- a/benchmark/src/parray_narytree.ml
+++ b/benchmark/src/parray_narytree.ml
@@ -1,4 +1,4 @@
-let default_base = 32
+let default_base = 5 (* chunk of size 32 *)
 
 let cast = Obj.magic
 
@@ -133,7 +133,7 @@ let iter (type a) (f : a -> unit) (s : a t) : unit =
   (* note: if it wasn't for performance, the two for loops could be factorized with a if inside the body *)
 
 (** [print s] prints the sequence [s] for debugging purpose. *)
-let print (type a) (f:a -> unit) (s : a t) : unit =
+let _print (type a) (f:a -> unit) (s : a t) : unit =
   let pr = Printf.printf in
   let b, d = params s in
   pr "<base=%d, depth=%d>" b d;
@@ -156,11 +156,12 @@ let print (type a) (f:a -> unit) (s : a t) : unit =
   pr "\n"
 
 
+(*
 module Test = struct
 
 let _test () =
   let pr = Printf.printf in
-  let print = print (fun x -> pr "%d" x) in
+  let print = _print (fun x -> pr "%d" x) in
 
  (*
   for n = 0 to 20 do
@@ -203,3 +204,4 @@ let _test () =
 
 end
 
+*)
\ No newline at end of file
diff --git a/benchmark/src/pdarray_narytree.ml b/benchmark/src/pdarray_narytree.ml
new file mode 100644
index 0000000000000000000000000000000000000000..64cd7f64c6e3d54ae1c605df2ec228816aee486f
--- /dev/null
+++ b/benchmark/src/pdarray_narytree.ml
@@ -0,0 +1,211 @@
+let default_base = 8 (* chunk of size 256 *)
+
+let cast = Obj.magic
+
+(** This data structure represents a tree of arity [2^base].
+    Each node is represented as an array of elements or as an array of subtrees.
+    The array at the first level stores two extra cells at its front,
+    one for storing the [base], and one for storing the [depth].
+    Other cells store either elements or pointers to arrays, each of
+    which stores elements or arrays of ...etc.
+    All arrays in depth store between 1 and 2^base elements, inclusive.
+    The root array stores the base and the depth, then between 0 and 2^base elements. *)
+
+type 'a t = 'a Darray.t  (* or ['a array array] or ['a array array array], etc, using Obj.magic casts *)
+
+let params (*[@inline]*) (s : 'a t) : int * int = (* returns [base] and [depth] *)
+  let s = cast s in (* cast for reading parameters *)
+  Darray.get s 0, Darray.get s 1
+
+let length (s : 'a t) : int =
+  let b, d = params s in
+  (* d is the depth of each item in s, in the sense that items have type ['a array^d] *)
+  let rec aux : 'a. ?offset:int -> int -> 'a t -> int =
+   fun (type a) ?(offset:int=0) (d:int) (s:a t) : int ->
+    (* Printf.printf "aux offset=%d depth=%d\n" offset d; *)
+    let nb_subtrees = Darray.length s - offset in
+    if d = 0 then nb_subtrees else begin
+      let size_of_full_subtree = 1 lsl (d * b) in
+      let items_in_full_subtrees = (nb_subtrees - 1) * size_of_full_subtree in
+      let subtree = Darray.get (cast s) (offset + nb_subtrees - 1) in
+      let items_in_last_subtree = aux (d - 1) subtree in
+      items_in_full_subtrees + items_in_last_subtree
+    end in
+  aux ~offset:2 d s
+
+(** [get s i] returns the element at index [i].
+    Requires [0 <= i < length t]. *)
+let get (s : 'a t) (i : int) : 'a =
+  (* assert (0 <= i && i < length s) *)
+  let b, d = params s in
+  let rec aux : 'a. ?offset:int -> int -> 'a t -> int -> 'a =
+   fun (type a) ?(offset:int=0) (d:int) (s:a t) (i:int) : a ->
+    (* Printf.printf "aux offset=%d depth=%d i=%d\n" offset d i; *)
+    if d = 0 then Darray.get s (offset + i) else begin
+      let w = d * b in (* [w] is the number of bits required to represent an index stored in one element of [s] *)
+      let id_subtree = i lsr w in
+      let size_of_full_subtree = 1 lsl w in
+      let index_in_subtree = i - id_subtree * size_of_full_subtree in
+      let subtree = Darray.get (cast s) (offset + id_subtree) in
+      aux (d-1) subtree index_in_subtree
+    end in
+  aux ~offset:2 d s i
+
+(** [set s i v] updates the element at index [i] with [v].
+    Requires [0 <= i < length t]. *)
+let set (s : 'a t) (i : int) (v : 'a) : 'a t =
+  (* assert (0 <= i && i < length s) *)
+  let b, d = params s in
+  let rec aux : 'a. ?offset:int -> int -> 'a t -> int -> 'a -> 'a t =
+    fun (type a) ?(offset:int=0) (d:int) (s:a t) (i:int) (v:a) : a t ->
+    (* Printf.printf "aux offset=%d depth=%d i=%d\n" offset d i;*)
+    if d = 0 then begin
+      Darray.set s (offset + i) v
+    end else begin
+      let w = d * b in (* [w] is the number of bits required to represent an index stored in one element of [s] *)
+      let id_subtree = i lsr w in
+      let index_in_subtree = i land ((1 lsl w) - 1) in
+      let subtree = Darray.get (cast s) (offset + id_subtree) in
+      let updated_subtree = aux (d-1) subtree index_in_subtree (cast v) in
+      Darray.set s (offset + id_subtree) (cast updated_subtree)
+    end in
+  aux ~offset:2 d s i v
+
+let depth_of_base_and_length (base:int) (n:int) : int =
+  (* returns the least value of [d] such that [n <= 2^((d+1)*b)], that is, [d = log_(2*b)(n)-1].
+     returns 0 if n=0. *)
+  let rec aux (i:int) : int =
+    if i = 0 then 0 else 1 + (aux (i lsr base)) in
+  if n = 0 then 0 else aux ((n-1) lsr base)
+
+(** [init ?base n f] return a sequence of length [n], with [f i] as [i]-th element.
+    The optional argument [base] can be used to specify that arrays used to
+    represent every node in the tree should have length [2^base].
+    Typically, [base] should be a value between [3] and [16]. *)
+let init ?(base:int = default_base) (n:int) (f: int -> 'a) : 'a t =
+  (* assert (1 <= base && base < 32) -- array size should not exceed 2^32 (could be up to max_array_size) *)
+  let b = base in
+  let d = depth_of_base_and_length base n in
+  (*Printf.printf "base=%d, depth=%d\n" b d;*)
+  if d = 0 then begin
+     (* assert (0 <= n && n <= 1 lsl base) *)
+     Darray.init (2 + n) (fun i ->
+       if i >= 2 then f i
+       else if i = 0 then cast b
+       else (* if i = 1 then *) cast d)
+  end else begin
+    let rec aux : 'a. (int -> 'a) -> int -> int -> int -> 'a t =
+      fun (type a) (f: int -> a) (d:int) (offset:int) (n:int) : a t ->
+      (*Printf.printf "aux depth=%d offset=%d n=%d\n" d offset n;*)
+      assert (n > 0);
+      if d = 0 then begin
+        Darray.init n (fun i -> f (offset + i))
+      end else begin
+        let w = d * b in
+        let nb_subtrees = 1 + ((n - 1) lsr w) in (* rounding up of [n / 2^w] *)
+        let size_of_full_subtree = 1 lsl w in
+        (*Printf.printf "  size_of_full_subtree=%d nb_subtrees=%d\n" size_of_full_subtree nb_subtrees;*)
+        cast (Darray.init nb_subtrees (fun i ->
+           let size_of_subtree =
+             if i < nb_subtrees - 1
+               then size_of_full_subtree
+               else n - (nb_subtrees - 1) * size_of_full_subtree
+             in
+          aux (cast f) (d - 1) (offset + i lsl w) size_of_subtree))
+      end in
+    let t = aux f d 0 n in
+    Darray.of_array_with_transfer (Array.append (cast [| b; d |]) (Darray.to_array_with_transfer t))
+  end
+
+(** [iter f s] applies [f] to each of the elements from the sequence [s]. *)
+let iter (type a) (f : a -> unit) (s : a t) : unit =
+  let _b, d = params s in
+  let rec aux : 'a. ('a -> unit) -> ?offset:int -> int -> 'a t -> unit =
+   fun (type a) (f : a->unit) ?(offset : int = 0) (d:int) (s:a t) : unit ->
+    if d = 0 then begin
+       for i = offset to Darray.length s - 1 do
+         f (Darray.get s i)
+       done
+    end else begin
+       for i = offset to Darray.length s - 1 do
+         let subtree = Darray.get (cast s) i in
+         aux (cast f) (d-1) subtree
+       done
+    end
+    in
+  aux f ~offset:2 d s
+  (* note: if it wasn't for performance, the two for loops could be factorized with a if inside the body *)
+
+(** [print s] prints the sequence [s] for debugging purpose. *)
+let _print (type a) (f:a -> unit) (s : a t) : unit =
+  let pr = Printf.printf in
+  let b, d = params s in
+  pr "<base=%d, depth=%d>" b d;
+  let rec aux : 'a. ('a -> unit) -> ?offset:int -> int -> 'a t -> unit =
+   fun (type a) (f:a -> unit) ?(offset : int = 0) (d:int) (s:a t) : unit ->
+    pr "[";
+    let r = Darray.length s in
+    for i = offset to r - 1 do
+      if d = 0 then begin
+        f (Darray.get s i)
+      end else begin
+        aux (cast f) (d-1) (Darray.get (cast s) i)
+      end;
+      if i < r - 1
+        then pr "; ";
+    done;
+    pr "]";
+    in
+  aux f ~offset:2 d s;
+  pr "\n"
+
+
+(*
+module Test = struct
+
+let _test () =
+  let pr = Printf.printf in
+  let print = _print (fun x -> pr "%d" x) in
+
+ (*
+  for n = 0 to 20 do
+    pr "-----make %d:\n" n;
+    let t : int t = init ~base:1 n (fun i -> i)in
+    print t;
+    pr "----->length:%d\n" (length t);
+
+  done;
+   *)
+
+  pr "-----make one:\n";
+  let n = 30 in
+  let t : int t = init ~base:1 n (fun i -> i)in
+  print t;
+
+  pr "-----length:%d\n" (length t);
+
+  (*
+
+  for i = 0 to n-1 do
+    pr "-----get %d: \n" i;
+    pr "==> %d\n" (get t i);
+  done;
+  pr "\n";
+  *)
+  (*
+  pr "-----set one:\n";
+   print (set t (n/2) (-1))
+  *)
+
+  pr "-----set:\n";
+  let r = ref t in
+  for i = 0 to n-1 do
+    r := set !r i (-i);
+    print !r;
+  done;
+  pr "\n"
+(**)
+
+end
+
+*)
\ No newline at end of file
diff --git a/benchmark/src/pdarray_narytree.mli b/benchmark/src/pdarray_narytree.mli
new file mode 100644
index 0000000000000000000000000000000000000000..6007eba778bbe87886d95d00b5683e750eed599f
--- /dev/null
+++ b/benchmark/src/pdarray_narytree.mli
@@ -0,0 +1,24 @@
+
+(** A persistent, fixed-length sequence of values of type ['a]. *)
+type 'a t
+
+(** The number of elements. *)
+val length: 'a t -> int
+
+(** [get s i] returns the element at index [i].
+    Requires [0 <= i < length t]. *)
+val get: 'a t -> int -> 'a
+
+(** [set s i v] updates the element at index [i] with [v].
+    Requires [0 <= i < length t]. *)
+val set: 'a t -> int -> 'a  -> 'a t
+
+(** [init ?base n f] return a sequence of length [n], with [f i] as [i]-th element.
+    The optional argument [base] can be used to specify that arrays used to
+    represent every node in the tree should have length [2^base].
+    Typically, [base] should be a value between [3] and [8]. *)
+val init: ?base:int -> int -> (int -> 'a) -> 'a t
+
+(** [iter f s] applies [f] to each of the elements from the sequence [s]. *)
+val iter: ('a -> unit) -> 'a t -> unit
+
diff --git a/benchmark/src/reach/Main.ml b/benchmark/src/reach/Main.ml
index 7bb02d7c088f44b0852f55c6e1933789e269678c..a723a4418a908b32dad4f71ea003c25986887990 100644
--- a/benchmark/src/reach/Main.ml
+++ b/benchmark/src/reach/Main.ml
@@ -152,6 +152,18 @@ let benchmark =
           sum := !sum + x
         done;
         sink !sum
+   end else if seq = "PdarrayNarytree" then begin
+      let s = Pdarray_narytree.init n (fun i -> i) in
+      fun () ->
+        (* This is the timed section. *)
+        (* Calls to [get i] for each target index. *)
+        let sum = ref 0 in
+        for index = 0 to k - 1 do
+          let j = a.(index) in
+          let x = Pdarray_narytree.get s j in
+          sum := !sum + x
+        done;
+        sink !sum
    end else if seq = "ESek_get" then begin
       let s =
         Construction.ebuild n in
diff --git a/benchmark/src/set/Main.ml b/benchmark/src/set/Main.ml
index fb23ec50096d30f30dcd49f3ed96e038da18909e..435213fa5db972d866e374454b205e305cc60f51 100644
--- a/benchmark/src/set/Main.ml
+++ b/benchmark/src/set/Main.ml
@@ -82,6 +82,7 @@ let benchmark =
   and basis =
     k
   and run () =
+
    if seq = "Array" then begin
 
       let s = Array.init n (fun i -> i) in
@@ -96,6 +97,20 @@ let benchmark =
         (* TODO: do we need this sink to avoid all the loop being compiled away? *)
         sink (Array.get s 0)
 
+   end else if seq = "Fmarray" then begin (* warning: very slow *)
+
+      let s = ref (Array.init n (fun i -> i)) in
+      fun () ->
+        (* This is the timed section. *)
+        (* Calls to [Fmarray.replace i] for each target index. *)
+        for index = 0 to k - 1 do
+          let j = a.(index) in
+          let value = index in
+          s := Fmarray.replace j value !s
+        done;
+        (* TODO: do we need this sink to avoid all the loop being compiled away? *)
+        sink (Array.get !s 0)
+
    end else if seq = "FmlibRbtrees" then begin
 
       (* initialization function for rbtrees *)
@@ -129,9 +144,21 @@ let benchmark =
           Numerichunk.set s j value;
         done
 
-   end else if seq = "ParrayNarytree" then begin
+   end else if seq = "ParrayNarytree5" then begin (* LATER: factorize code for various chunksize *)
+
+      let s = ref (Parray_narytree.init ~base:5 n (fun i -> i)) in
+      fun () ->
+        (* This is the timed section. *)
+        (* Calls to [replace i] for each target index. *)
+        for index = 0 to k - 1 do
+          let j = a.(index) in
+          let value = index in
+          s := Parray_narytree.set !s j value;
+        done
+
+   end else if seq = "ParrayNarytree8" then begin
 
-      let s = ref (Parray_narytree.init n (fun i -> i)) in
+      let s = ref (Parray_narytree.init ~base:5 n (fun i -> i)) in
       fun () ->
         (* This is the timed section. *)
         (* Calls to [replace i] for each target index. *)
@@ -141,6 +168,54 @@ let benchmark =
           s := Parray_narytree.set !s j value;
         done
 
+   end else if seq = "PdarrayNarytree5" then begin (* LATER: factorize code for various chunksize *)
+
+      let s = ref (Pdarray_narytree.init ~base:5 n (fun i -> i)) in
+      fun () ->
+        (* This is the timed section. *)
+        (* Calls to [replace i] for each target index. *)
+        for index = 0 to k - 1 do
+          let j = a.(index) in
+          let value = index in
+          s := Pdarray_narytree.set !s j value;
+        done
+
+   end else if seq = "PdarrayNarytree7" then begin
+
+      let s = ref (Pdarray_narytree.init ~base:7 n (fun i -> i)) in
+      fun () ->
+        (* This is the timed section. *)
+        (* Calls to [replace i] for each target index. *)
+        for index = 0 to k - 1 do
+          let j = a.(index) in
+          let value = index in
+          s := Pdarray_narytree.set !s j value;
+        done
+
+   end else if seq = "PdarrayNarytree8" then begin
+
+      let s = ref (Pdarray_narytree.init ~base:8 n (fun i -> i)) in
+      fun () ->
+        (* This is the timed section. *)
+        (* Calls to [replace i] for each target index. *)
+        for index = 0 to k - 1 do
+          let j = a.(index) in
+          let value = index in
+          s := Pdarray_narytree.set !s j value;
+        done
+
+   end else if seq = "PdarrayNarytree10" then begin
+
+      let s = ref (Pdarray_narytree.init ~base:10 n (fun i -> i)) in
+      fun () ->
+        (* This is the timed section. *)
+        (* Calls to [replace i] for each target index. *)
+        for index = 0 to k - 1 do
+          let j = a.(index) in
+          let value = index in
+          s := Pdarray_narytree.set !s j value;
+        done
+
    end else if seq = "ESek" then begin
 
       let s = Construction.ebuild n in