diff --git a/Changes b/Changes
index 76f4a5c80..67c75b91d 100644
--- a/Changes
+++ b/Changes
@@ -177,6 +177,9 @@ Working version
 - GPR#2159, MPR#7874: annotate {String,Bytes}.equal as being [@@noalloc].
   (Pierre-Marie Pédrot, review by Nicolás Ojeda Bär)
 
+- GPR#1936: Add module Float.Array
+  (Damien Doligez, review by Xavier Clerc and Alain Frisch)
+
 ### Other libraries:
 
 - GPR#1061: Add ?follow parameter to Unix.link. This allows hardlinking
diff --git a/stdlib/Compflags b/stdlib/Compflags
index d22afa400..0f3138cd8 100755
--- a/stdlib/Compflags
+++ b/stdlib/Compflags
@@ -30,5 +30,6 @@ case $1 in
       echo ' -w Ae';;
   stdlib__scanf.cmx|stdlib__scanf.p.cmx) echo ' -inline 9';;
   *Labels.cm[ox]|*Labels.p.cmx) echo ' -nolabels -no-alias-deps';;
+  stdlib__float.cm[ox]|stdlib__float.p.cmx) echo ' -nolabels -no-alias-deps';;
   *) echo ' ';;
 esac
diff --git a/stdlib/Makefile b/stdlib/Makefile
index 5afd1bb34..eecbbd33e 100644
--- a/stdlib/Makefile
+++ b/stdlib/Makefile
@@ -44,7 +44,7 @@ OBJS=camlinternalFormatBasics.cmo stdlib.cmo $(OTHERS)
 OTHERS= $(P)pervasives.cmo $(P)seq.cmo $(P)option.cmo $(P)result.cmo \
   $(P)bool.cmo $(P)char.cmo $(P)uchar.cmo $(P)sys.cmo $(P)list.cmo \
   $(P)bytes.cmo $(P)string.cmo $(P)fun.cmo $(P)unit.cmo \
-  $(P)marshal.cmo $(P)obj.cmo $(P)float.cmo $(P)array.cmo \
+  $(P)marshal.cmo $(P)obj.cmo $(P)array.cmo $(P)float.cmo \
   $(P)int.cmo $(P)int32.cmo $(P)int64.cmo $(P)nativeint.cmo \
   $(P)lexing.cmo $(P)parsing.cmo \
   $(P)set.cmo $(P)map.cmo $(P)stack.cmo $(P)queue.cmo \
diff --git a/stdlib/array.mli b/stdlib/array.mli
index 77824b54a..6e3d40530 100644
--- a/stdlib/array.mli
+++ b/stdlib/array.mli
@@ -27,7 +27,7 @@ external get : 'a array -> int -> 'a = "%array_safe_get"
    The last element has number [Array.length a - 1].
    You can also write [a.(n)] instead of [Array.get a n].
 
-   Raise [Invalid_argument "index out of bounds"]
+   Raise [Invalid_argument]
    if [n] is outside the range 0 to [(Array.length a - 1)]. *)
 
 external set : 'a array -> int -> 'a -> unit = "%array_safe_set"
@@ -35,7 +35,7 @@ external set : 'a array -> int -> 'a -> unit = "%array_safe_set"
    element number [n] with [x].
    You can also write [a.(n) <- x] instead of [Array.set a n x].
 
-   Raise [Invalid_argument "index out of bounds"]
+   Raise [Invalid_argument]
    if [n] is outside the range 0 to [Array.length a - 1]. *)
 
 external make : int -> 'a -> 'a array = "caml_make_vect"
@@ -93,7 +93,10 @@ val create_matrix : int -> int -> 'a -> 'a array array
 
 val append : 'a array -> 'a array -> 'a array
 (** [Array.append v1 v2] returns a fresh array containing the
-   concatenation of the arrays [v1] and [v2]. *)
+   concatenation of the arrays [v1] and [v2].
+
+   Raise [Invalid_argument] if
+   [Array.length v1 + Array.length v2 > Sys.max_array_length]. *)
 
 val concat : 'a array list -> 'a array
 (** Same as {!Array.append}, but concatenates a list of arrays. *)
@@ -103,7 +106,7 @@ val sub : 'a array -> int -> int -> 'a array
    containing the elements number [start] to [start + len - 1]
    of array [a].
 
-   Raise [Invalid_argument "Array.sub"] if [start] and [len] do not
+   Raise [Invalid_argument] if [start] and [len] do not
    designate a valid subarray of [a]; that is, if
    [start < 0], or [len < 0], or [start + len > Array.length a]. *)
 
@@ -115,7 +118,7 @@ val fill : 'a array -> int -> int -> 'a -> unit
 (** [Array.fill a ofs len x] modifies the array [a] in place,
    storing [x] in elements number [ofs] to [ofs + len - 1].
 
-   Raise [Invalid_argument "Array.fill"] if [ofs] and [len] do not
+   Raise [Invalid_argument] if [ofs] and [len] do not
    designate a valid subarray of [a]. *)
 
 val blit : 'a array -> int -> 'a array -> int -> int -> unit
@@ -125,7 +128,7 @@ val blit : 'a array -> int -> 'a array -> int -> int -> unit
    [v1] and [v2] are the same array, and the source and
    destination chunks overlap.
 
-   Raise [Invalid_argument "Array.blit"] if [o1] and [len] do not
+   Raise [Invalid_argument] if [o1] and [len] do not
    designate a valid subarray of [v1], or if [o2] and [len] do not
    designate a valid subarray of [v2]. *)
 
@@ -134,7 +137,10 @@ val to_list : 'a array -> 'a list
 
 val of_list : 'a list -> 'a array
 (** [Array.of_list l] returns a fresh array containing the elements
-   of [l]. *)
+   of [l].
+
+   Raise [Invalid_argument] if the length of [l] is greater than
+   [Sys.max_array_length].*)
 
 
 (** {1 Iterators} *)
@@ -204,13 +210,14 @@ val exists : ('a -> bool) -> 'a array -> bool
     @since 4.03.0 *)
 
 val mem : 'a -> 'a array -> bool
-(** [mem a l] is true if and only if [a] is equal
-   to an element of [l].
-   @since 4.03.0 *)
+(** [mem a l] is true if and only if [a] is structurally equal
+    to an element of [l] (i.e. there is an [x] in [l] such that
+    [compare a x = 0]).
+    @since 4.03.0 *)
 
 val memq : 'a -> 'a array -> bool
 (** Same as {!Array.mem}, but uses physical equality instead of structural
-   equality to compare array elements.
+   equality to compare elements.
    @since 4.03.0 *)
 
 
@@ -223,8 +230,7 @@ val sort : ('a -> 'a -> int) -> 'a array -> unit
    compare as equal, a positive integer if the first is greater,
    and a negative integer if the first is smaller (see below for a
    complete specification).  For example, {!Stdlib.compare} is
-   a suitable comparison function, provided there are no floating-point
-   NaN values in the data.  After calling [Array.sort], the
+   a suitable comparison function.  After calling [Array.sort], the
    array is sorted in place in increasing order.
    [Array.sort] is guaranteed to run in constant heap space
    and (at most) logarithmic stack space.
@@ -234,7 +240,7 @@ val sort : ('a -> 'a -> int) -> 'a array -> unit
 
    Specification of the comparison function:
    Let [a] be the array and [cmp] the comparison function.  The following
-   must be true for all x, y, z in a :
+   must be true for all [x], [y], [z] in [a] :
 -   [cmp x y] > 0 if and only if [cmp y x] < 0
 -   if [cmp x y] >= 0 and [cmp y z] >= 0 then [cmp x z] >= 0
 
@@ -248,8 +254,8 @@ val stable_sort : ('a -> 'a -> int) -> 'a array -> unit
    elements that compare equal are kept in their original order) and
    not guaranteed to run in constant heap space.
 
-   The current implementation uses Merge Sort. It uses [n/2]
-   words of heap space, where [n] is the length of the array.
+   The current implementation uses Merge Sort. It uses a temporary
+   array of length [n/2], where [n] is the length of the array.
    It is usually faster than the current implementation of {!Array.sort}.
 *)
 
diff --git a/stdlib/float.ml b/stdlib/float.ml
index e5a4d5271..3145f1c66 100644
--- a/stdlib/float.ml
+++ b/stdlib/float.ml
@@ -148,11 +148,364 @@ external seeded_hash_param : int -> int -> int -> float -> int
 let hash x = seeded_hash_param 10 100 0 x
 
 module Array = struct
+
   type t = floatarray
-  external create : int -> t = "caml_floatarray_create"
+
   external length : t -> int = "%floatarray_length"
   external get : t -> int -> float = "%floatarray_safe_get"
   external set : t -> int -> float -> unit = "%floatarray_safe_set"
+  external create : int -> t = "caml_floatarray_create"
   external unsafe_get : t -> int -> float = "%floatarray_unsafe_get"
   external unsafe_set : t -> int -> float -> unit = "%floatarray_unsafe_set"
+
+  let unsafe_fill a ofs len v =
+    for i = ofs to ofs + len - 1 do unsafe_set a i v done
+
+  let unsafe_blit src sofs dst dofs len =
+    for i = 0 to len - 1 do
+      unsafe_set dst (dofs + i) (unsafe_get src (sofs + i))
+    done
+
+  let check a ofs len msg =
+    if ofs < 0 || len < 0 || ofs + len < 0 || ofs + len > length a then
+      invalid_arg msg
+
+  let make n v =
+    let result = create n in
+    unsafe_fill result 0 n v;
+    result
+
+  let init l f =
+    if l < 0 then invalid_arg "Float.Array.init"
+    else
+      let res = create l in
+      for i = 0 to l - 1 do
+        unsafe_set res i (f i)
+      done;
+      res
+
+  let append a1 a2 =
+    let l1 = length a1 in
+    let l2 = length a2 in
+    let result = create (l1 + l2) in
+    unsafe_blit a1 0 result 0 l1;
+    unsafe_blit a2 0 result l1 l2;
+    result
+
+  (* next 3 functions: modified copy of code from string.ml *)
+  let ensure_ge (x:int) y =
+    if x >= y then x else invalid_arg "Float.Array.concat"
+
+  let rec sum_lengths acc = function
+    | [] -> acc
+    | hd :: tl -> sum_lengths (ensure_ge (length hd + acc) acc) tl
+
+  let concat l =
+    let len = sum_lengths 0 l in
+    let result = create len in
+    let rec loop l i =
+      match l with
+      | [] -> assert (i = len)
+      | hd :: tl ->
+        let hlen = length hd in
+        unsafe_blit hd 0 result i hlen;
+        loop tl (i + hlen)
+    in
+    loop l 0;
+    result
+
+  let sub a ofs len =
+    check a ofs len "Float.Array.sub";
+    let result = create len in
+    unsafe_blit a ofs result 0 len;
+    result
+
+  let copy a =
+    let l = length a in
+    let result = create l in
+    unsafe_blit a 0 result 0 l;
+    result
+
+  let fill a ofs len v =
+    check a ofs len "Float.Array.fill";
+    unsafe_fill a ofs len v
+
+  let blit src sofs dst dofs len =
+    check src sofs len "Float.array.blit";
+    check dst dofs len "Float.array.blit";
+    unsafe_blit src sofs dst dofs len
+
+  let to_list a =
+    List.init (length a) (unsafe_get a)
+
+  let of_list l =
+    let result = create (List.length l) in
+    let rec fill i l =
+      match l with
+      | [] -> result
+      | h :: t -> unsafe_set result i h; fill (i + 1) t
+    in
+    fill 0 l
+
+  (* duplicated from array.ml *)
+  let iter f a =
+    for i = 0 to length a - 1 do f (unsafe_get a i) done
+
+  (* duplicated from array.ml *)
+  let iter2 f a b =
+    if length a <> length b then
+      invalid_arg "Float.Array.iter2: arrays must have the same length"
+    else
+      for i = 0 to length a - 1 do f (unsafe_get a i) (unsafe_get b i) done
+
+  let map f a =
+    let l = length a in
+    let r = create l in
+    for i = 0 to l - 1 do
+      unsafe_set r i (f (unsafe_get a i))
+    done;
+    r
+
+  let map2 f a b =
+    let la = length a in
+    let lb = length b in
+    if la <> lb then
+      invalid_arg "Float.Array.map2: arrays must have the same length"
+    else begin
+      let r = create la in
+      for i = 0 to la - 1 do
+        unsafe_set r i (f (unsafe_get a i) (unsafe_get b i))
+      done;
+      r
+    end
+
+  (* duplicated from array.ml *)
+  let iteri f a =
+    for i = 0 to length a - 1 do f i (unsafe_get a i) done
+
+  let mapi f a =
+    let l = length a in
+    let r = create l in
+    for i = 0 to l - 1 do
+      unsafe_set r i (f i (unsafe_get a i))
+    done;
+    r
+
+  (* duplicated from array.ml *)
+  let fold_left f x a =
+    let r = ref x in
+    for i = 0 to length a - 1 do
+      r := f !r (unsafe_get a i)
+    done;
+    !r
+
+  (* duplicated from array.ml *)
+  let fold_right f a x =
+    let r = ref x in
+    for i = length a - 1 downto 0 do
+      r := f (unsafe_get a i) !r
+    done;
+    !r
+
+  (* duplicated from array.ml *)
+  let exists p a =
+    let n = length a in
+    let rec loop i =
+      if i = n then false
+      else if p (unsafe_get a i) then true
+      else loop (i + 1) in
+    loop 0
+
+  (* duplicated from array.ml *)
+  let for_all p a =
+    let n = length a in
+    let rec loop i =
+      if i = n then true
+      else if p (unsafe_get a i) then loop (i + 1)
+      else false in
+    loop 0
+
+  (* duplicated from array.ml *)
+  let mem x a =
+    let n = length a in
+    let rec loop i =
+      if i = n then false
+      else if compare (unsafe_get a i) x = 0 then true
+      else loop (i + 1)
+    in
+    loop 0
+
+  (* mostly duplicated from array.ml, but slightly different *)
+  let mem_ieee x a =
+    let n = length a in
+    let rec loop i =
+      if i = n then false
+      else if x = (unsafe_get a i) then true
+      else loop (i + 1)
+    in
+    loop 0
+
+  (* duplicated from array.ml *)
+  exception Bottom of int
+  let sort cmp a =
+    let maxson l i =
+      let i31 = i+i+i+1 in
+      let x = ref i31 in
+      if i31+2 < l then begin
+        if cmp (get a i31) (get a (i31+1)) < 0 then x := i31+1;
+        if cmp (get a !x) (get a (i31+2)) < 0 then x := i31+2;
+        !x
+      end else
+        if i31+1 < l && cmp (get a i31) (get a (i31+1)) < 0
+        then i31+1
+        else if i31 < l then i31 else raise (Bottom i)
+    in
+    let rec trickledown l i e =
+      let j = maxson l i in
+      if cmp (get a j) e > 0 then begin
+        set a i (get a j);
+        trickledown l j e;
+      end else begin
+        set a i e;
+      end;
+    in
+    let trickle l i e = try trickledown l i e with Bottom i -> set a i e in
+    let rec bubbledown l i =
+      let j = maxson l i in
+      set a i (get a j);
+      bubbledown l j
+    in
+    let bubble l i = try bubbledown l i with Bottom i -> i in
+    let rec trickleup i e =
+      let father = (i - 1) / 3 in
+      assert (i <> father);
+      if cmp (get a father) e < 0 then begin
+        set a i (get a father);
+        if father > 0 then trickleup father e else set a 0 e;
+      end else begin
+        set a i e;
+      end;
+    in
+    let l = length a in
+    for i = (l + 1) / 3 - 1 downto 0 do trickle l i (get a i); done;
+    for i = l - 1 downto 2 do
+      let e = (get a i) in
+      set a i (get a 0);
+      trickleup (bubble i 0) e;
+    done;
+    if l > 1 then (let e = (get a 1) in set a 1 (get a 0); set a 0 e)
+
+  (* duplicated from array.ml, except for the call to [create] *)
+  let cutoff = 5
+  let stable_sort cmp a =
+    let merge src1ofs src1len src2 src2ofs src2len dst dstofs =
+      let src1r = src1ofs + src1len and src2r = src2ofs + src2len in
+      let rec loop i1 s1 i2 s2 d =
+        if cmp s1 s2 <= 0 then begin
+          set dst d s1;
+          let i1 = i1 + 1 in
+          if i1 < src1r then
+            loop i1 (get a i1) i2 s2 (d + 1)
+          else
+            blit src2 i2 dst (d + 1) (src2r - i2)
+        end else begin
+          set dst d s2;
+          let i2 = i2 + 1 in
+          if i2 < src2r then
+            loop i1 s1 i2 (get src2 i2) (d + 1)
+          else
+            blit a i1 dst (d + 1) (src1r - i1)
+        end
+      in loop src1ofs (get a src1ofs) src2ofs (get src2 src2ofs) dstofs;
+    in
+    let isortto srcofs dst dstofs len =
+      for i = 0 to len - 1 do
+        let e = (get a (srcofs + i)) in
+        let j = ref (dstofs + i - 1) in
+        while (!j >= dstofs && cmp (get dst !j) e > 0) do
+          set dst (!j + 1) (get dst !j);
+          decr j;
+        done;
+        set dst (!j + 1) e;
+      done;
+    in
+    let rec sortto srcofs dst dstofs len =
+      if len <= cutoff then isortto srcofs dst dstofs len else begin
+        let l1 = len / 2 in
+        let l2 = len - l1 in
+        sortto (srcofs + l1) dst (dstofs + l1) l2;
+        sortto srcofs a (srcofs + l2) l1;
+        merge (srcofs + l2) l1 dst (dstofs + l1) l2 dst dstofs;
+      end;
+    in
+    let l = length a in
+    if l <= cutoff then isortto 0 a 0 l else begin
+      let l1 = l / 2 in
+      let l2 = l - l1 in
+      let t = create l2 in
+      sortto l1 t 0 l2;
+      sortto 0 a l2 l1;
+      merge l2 l1 t 0 l2 a 0;
+    end
+
+  let fast_sort = stable_sort
+
+  (* duplicated from array.ml *)
+  let to_seq a =
+    let rec aux i () =
+      if i < length a
+      then
+        let x = unsafe_get a i in
+        Seq.Cons (x, aux (i+1))
+      else Seq.Nil
+    in
+    aux 0
+
+  (* duplicated from array.ml *)
+  let to_seqi a =
+    let rec aux i () =
+      if i < length a
+      then
+        let x = unsafe_get a i in
+        Seq.Cons ((i,x), aux (i+1))
+      else Seq.Nil
+    in
+    aux 0
+
+  (* mostly duplicated from array.ml *)
+  let of_rev_list l =
+    let len = List.length l in
+    let a = create len in
+    let rec fill i = function
+        [] -> a
+      | hd::tl -> unsafe_set a i hd; fill (i-1) tl
+    in
+    fill (len-1) l
+
+  (* duplicated from array.ml *)
+  let of_seq i =
+    let l = Seq.fold_left (fun acc x -> x::acc) [] i in
+    of_rev_list l
+
+
+  let map_to_array f a =
+    let l = length a in
+    if l = 0 then [| |] else begin
+      let r = Array.make l (f (unsafe_get a 0)) in
+      for i = 1 to l - 1 do
+        Array.unsafe_set r i (f (unsafe_get a i))
+      done;
+      r
+    end
+
+  let map_from_array f a =
+    let l = Array.length a in
+    let r = create l in
+    for i = 0 to l - 1 do
+      unsafe_set r i (f (Array.unsafe_get a i))
+    done;
+    r
+
 end
+
+module ArrayLabels = Array
diff --git a/stdlib/float.mli b/stdlib/float.mli
index 38386e025..70d6ab8a4 100644
--- a/stdlib/float.mli
+++ b/stdlib/float.mli
@@ -389,11 +389,270 @@ val hash: t -> int
 (** The hash function for floating-point numbers. *)
 
 module Array : sig
+
   type t = floatarray
-  external create : int -> t = "caml_floatarray_create"
-  external length : t -> int = "%floatarray_length"
-  external get : t -> int -> float = "%floatarray_safe_get"
-  external set : t -> int -> float -> unit = "%floatarray_safe_set"
+  (** The type of float arrays with packed representation.  @since 4.08.0 *)
+
+  val length : t -> int
+  (** Return the length (number of elements) of the given floatarray. *)
+
+  val get : t -> int -> float
+  (** [get a n] returns the element number [n] of floatarray [a].
+
+      Raise [Invalid_argument] if [n] is outside the range 0 to
+      [(length a - 1)]. *)
+
+  val set : t -> int -> float -> unit
+  (** [set a n x] modifies floatarray [a] in place, replacing element
+      number [n] with [x].
+
+      Raise [Invalid_argument] if [n] is outside the range 0 to
+      [(length a - 1)]. *)
+
+  val make : int -> float -> t
+  (** [make n x] returns a fresh floatarray of length [n], initialized with [x].
+
+      Raise [Invalid_argument] if [n < 0] or [n > Sys.max_floatarray_length]. *)
+
+  val create : int -> t
+  (** [create n] returns a fresh floatarray of length [n],
+      with uninitialized data.
+
+      Raise [Invalid_argument] if [n < 0] or [n > Sys.max_floatarray_length]. *)
+
+  val init : int -> (int -> float) -> t
+  (** [init n f] returns a fresh floatarray of length [n],
+     with element number [i] initialized to the result of [f i].
+     In other terms, [init n f] tabulates the results of [f]
+     applied to the integers [0] to [n-1].
+
+     Raise [Invalid_argument] if [n < 0] or [n > Sys.max_floatarray_length]. *)
+
+  val append : t -> t -> t
+  (** [append v1 v2] returns a fresh floatarray containing the
+     concatenation of the floatarrays [v1] and [v2].
+
+     Raise [Invalid_argument] if
+     [length v1 + length v2 > Sys.max_floatarray_length]. *)
+
+  val concat : t list -> t
+  (** Same as {!append}, but concatenates a list of floatarrays. *)
+
+  val sub : t -> int -> int -> t
+  (** [sub a start len] returns a fresh floatarray of length [len],
+     containing the elements number [start] to [start + len - 1]
+     of floatarray [a].
+
+     Raise [Invalid_argument] if [start] and [len] do not
+     designate a valid subarray of [a]; that is, if
+     [start < 0], or [len < 0], or [start + len > length a]. *)
+
+  val copy : t -> t
+  (** [copy a] returns a copy of [a], that is, a fresh floatarray
+      containing the same elements as [a]. *)
+
+  val fill : t -> int -> int -> float -> unit
+  (** [fill a ofs len x] modifies the floatarray [a] in place,
+     storing [x] in elements number [ofs] to [ofs + len - 1].
+
+     Raise [Invalid_argument] if [ofs] and [len] do not
+     designate a valid subarray of [a]. *)
+
+  val blit : t -> int -> t -> int -> int -> unit
+  (** [blit v1 o1 v2 o2 len] copies [len] elements
+      from floatarray [v1], starting at element number [o1], to floatarray [v2],
+      starting at element number [o2]. It works correctly even if
+      [v1] and [v2] are the same floatarray, and the source and
+      destination chunks overlap.
+
+      Raise [Invalid_argument] if [o1] and [len] do not
+      designate a valid subarray of [v1], or if [o2] and [len] do not
+      designate a valid subarray of [v2]. *)
+
+  val to_list : t -> float list
+  (** [to_list a] returns the list of all the elements of [a]. *)
+
+  val of_list : float list -> t
+  (** [of_list l] returns a fresh floatarray containing the elements
+      of [l].
+
+      Raise [Invalid_argument] if the length of [l] is greater than
+      [Sys.max_floatarray_length].*)
+
+  (** {2 Iterators} *)
+
+  val iter : (float -> unit) -> t -> unit
+  (** [iter f a] applies function [f] in turn to all
+      the elements of [a].  It is equivalent to
+      [f a.(0); f a.(1); ...; f a.(length a - 1); ()]. *)
+
+  val iteri : (int -> float -> unit) -> t -> unit
+  (** Same as {!iter}, but the
+      function is applied with the index of the element as first argument,
+      and the element itself as second argument. *)
+
+  val map : (float -> float) -> t -> t
+  (** [map f a] applies function [f] to all the elements of [a],
+      and builds a floatarray with the results returned by [f]. *)
+
+  val mapi : (int -> float -> float) -> t -> t
+  (** Same as {!map}, but the
+      function is applied to the index of the element as first argument,
+      and the element itself as second argument. *)
+
+  val fold_left : ('a -> float -> 'a) -> 'a -> t -> 'a
+  (** [fold_left f x a] computes
+      [f (... (f (f x a.(0)) a.(1)) ...) a.(n-1)],
+      where [n] is the length of the floatarray [a]. *)
+
+  val fold_right : (float -> 'a -> 'a) -> t -> 'a -> 'a
+  (** [fold_right f a x] computes
+      [f a.(0) (f a.(1) ( ... (f a.(n-1) x) ...))],
+      where [n] is the length of the floatarray [a]. *)
+
+  (** {2 Iterators on two arrays} *)
+
+  val iter2 : (float -> float -> unit) -> t -> t -> unit
+  (** [Array.iter2 f a b] applies function [f] to all the elements of [a]
+      and [b].
+      Raise [Invalid_argument] if the floatarrays are not the same size. *)
+
+  val map2 : (float -> float -> float) -> t -> t -> t
+  (** [map2 f a b] applies function [f] to all the elements of [a]
+      and [b], and builds a floatarray with the results returned by [f]:
+      [[| f a.(0) b.(0); ...; f a.(length a - 1) b.(length b - 1)|]].
+      Raise [Invalid_argument] if the floatarrays are not the same size. *)
+
+  (** {2 Array scanning} *)
+
+  val for_all : (float -> bool) -> t -> bool
+  (** [for_all p [|a1; ...; an|]] checks if all elements of the floatarray
+      satisfy the predicate [p]. That is, it returns
+      [(p a1) && (p a2) && ... && (p an)]. *)
+
+  val exists : (float -> bool) -> t -> bool
+  (** [exists p [|a1; ...; an|]] checks if at least one element of
+      the floatarray satisfies the predicate [p]. That is, it returns
+      [(p a1) || (p a2) || ... || (p an)]. *)
+
+  val mem : float -> t -> bool
+  (** [mem a l] is true if and only if there is an element of [l] that is
+      structurally equal to [a], i.e. there is an [x] in [l] such
+      that [compare a x = 0]. *)
+
+  val mem_ieee : float -> t -> bool
+  (** Same as {!mem}, but uses IEEE equality instead of structural equality. *)
+
+  (** {2 Sorting} *)
+
+  val sort : (float -> float -> int) -> t -> unit
+  (** Sort a floatarray in increasing order according to a comparison
+      function.  The comparison function must return 0 if its arguments
+      compare as equal, a positive integer if the first is greater,
+      and a negative integer if the first is smaller (see below for a
+      complete specification).  For example, {!Pervasives.compare} is
+      a suitable comparison function.  After calling [sort], the
+      array is sorted in place in increasing order.
+      [sort] is guaranteed to run in constant heap space
+      and (at most) logarithmic stack space.
+
+      The current implementation uses Heap Sort.  It runs in constant
+      stack space.
+
+      Specification of the comparison function:
+      Let [a] be the floatarray and [cmp] the comparison function. The following
+      must be true for all [x], [y], [z] in [a] :
+-      [cmp x y] > 0 if and only if [cmp y x] < 0
+-      if [cmp x y] >= 0 and [cmp y z] >= 0 then [cmp x z] >= 0
+
+      When [sort] returns, [a] contains the same elements as before,
+      reordered in such a way that for all i and j valid indices of [a] :
+-      [cmp a.(i) a.(j)] >= 0 if and only if i >= j
+  *)
+
+  val stable_sort : (float -> float -> int) -> t -> unit
+  (** Same as {!sort}, but the sorting algorithm is stable (i.e.
+       elements that compare equal are kept in their original order) and
+       not guaranteed to run in constant heap space.
+
+       The current implementation uses Merge Sort. It uses a temporary
+       floatarray of length [n/2], where [n] is the length of the floatarray.
+       It is usually faster than the current implementation of {!sort}. *)
+
+  val fast_sort : (float -> float -> int) -> t -> unit
+  (** Same as {!sort} or {!stable_sort}, whichever is faster
+      on typical input. *)
+
+  (** {2 Iterators} *)
+
+  val to_seq : t -> float Seq.t
+  (** Iterate on the floatarray, in increasing order. Modifications of the
+      floatarray during iteration will be reflected in the iterator. *)
+
+  val to_seqi : t -> (int * float) Seq.t
+  (** Iterate on the floatarray, in increasing order, yielding indices along
+      elements. Modifications of the floatarray during iteration will be
+      reflected in the iterator. *)
+
+  val of_seq : float Seq.t -> t
+  (** Create an array from the generator. *)
+
+
+  val map_to_array : (float -> 'a) -> t -> 'a array
+  (** [map_to_array f a] applies function [f] to all the elements of [a],
+      and builds an array with the results returned by [f]:
+      [[| f a.(0); f a.(1); ...; f a.(length a - 1) |]]. *)
+
+  val map_from_array : ('a -> float) -> 'a array -> t
+  (** [map_from_array f a] applies function [f] to all the elements of [a],
+      and builds a floatarray with the results returned by [f]. *)
+
+  (** {2 Undocumented functions} *)
+
+  (* These functions are for system use only. Do not call directly. *)
+  external unsafe_get : t -> int -> float = "%floatarray_unsafe_get"
+  external unsafe_set : t -> int -> float -> unit = "%floatarray_unsafe_set"
+end
+
+module ArrayLabels : sig
+
+  type t = floatarray
+  val length : t -> int
+  val get : t -> int -> float
+  val set : t -> int -> float -> unit
+  val make : int -> float -> t
+  val create : int -> t
+  val init : int -> f:(int -> float) -> t
+  val append : t -> t -> t
+  val concat : t list -> t
+  val sub : t -> pos:int -> len:int -> t
+  val copy : t -> t
+  val fill : t -> pos:int -> len:int -> float -> unit
+  val blit : src:t -> src_pos:int -> dst:t -> dst_pos:int -> len:int -> unit
+  val to_list : t -> float list
+  val of_list : float list -> t
+  val iter : f:(float -> unit) -> t -> unit
+  val iteri : f:(int -> float -> unit) -> t -> unit
+  val map : f:(float -> float) -> t -> t
+  val mapi : f:(int -> float -> float) -> t -> t
+  val fold_left : f:('a -> float -> 'a) -> init:'a -> t -> 'a
+  val fold_right : f:(float -> 'a -> 'a) -> t -> init:'a -> 'a
+  val iter2 : f:(float -> float -> unit) -> t -> t -> unit
+  val map2 : f:(float -> float -> float) -> t -> t -> t
+  val for_all : f:(float -> bool) -> t -> bool
+  val exists : f:(float -> bool) -> t -> bool
+  val mem : float -> set:t -> bool
+  val mem_ieee : float -> set:t -> bool
+  val sort : cmp:(float -> float -> int) -> t -> unit
+  val stable_sort : cmp:(float -> float -> int) -> t -> unit
+  val fast_sort : cmp:(float -> float -> int) -> t -> unit
+  val to_seq : t -> float Seq.t
+  val to_seqi : t -> (int * float) Seq.t
+  val of_seq : float Seq.t -> t
+  val map_to_array : f:(float -> 'a) -> t -> 'a array
+  val map_from_array : f:('a -> float) -> 'a array -> t
+
+  (* These functions are for system use only. Do not call directly. *)
   external unsafe_get : t -> int -> float = "%floatarray_unsafe_get"
   external unsafe_set : t -> int -> float -> unit = "%floatarray_unsafe_set"
 end
diff --git a/stdlib/sys.mli b/stdlib/sys.mli
index 228849501..370449ecf 100644
--- a/stdlib/sys.mli
+++ b/stdlib/sys.mli
@@ -145,9 +145,16 @@ val max_string_length : int
 (** Maximum length of strings and byte sequences. *)
 
 val max_array_length : int
-(** Maximum length of a normal array.  The maximum length of a float
-    array is [max_array_length/2] on 32-bit machines and
-    [max_array_length] on 64-bit machines. *)
+(** Maximum length of a normal array (i.e. any array whose elements are
+    not of type [float]). The maximum length of a [float array]
+    is [max_floatarray_length] if OCaml was configured with
+    [--enable-flat-float-array] and [max_array_length] if configured
+    with [--disable-flat-float-array]. *)
+
+val max_floatarray_length : int
+(** Maximum length of a floatarray. This is also the maximum length of
+    a [float array] when OCaml is configured with
+    [--enable-flat-float-array]. *)
 
 external runtime_variant : unit -> string = "caml_runtime_variant"
 (** Return the name of the runtime variant the program is running on.
diff --git a/stdlib/sys.mlp b/stdlib/sys.mlp
index 62c84fb18..5d997e8ad 100644
--- a/stdlib/sys.mlp
+++ b/stdlib/sys.mlp
@@ -45,6 +45,7 @@ let unix = unix ()
 let win32 = win32 ()
 let cygwin = cygwin ()
 let max_array_length = max_wosize ()
+let max_floatarray_length = max_array_length / (64 / word_size)
 let max_string_length = word_size / 8 * max_array_length - 1
 external runtime_variant : unit -> string = "caml_runtime_variant"
 external runtime_parameters : unit -> string = "caml_runtime_parameters"
diff --git a/testsuite/tests/lib-floatarray/floatarray.ml b/testsuite/tests/lib-floatarray/floatarray.ml
new file mode 100644
index 000000000..763b2c7b5
--- /dev/null
+++ b/testsuite/tests/lib-floatarray/floatarray.ml
@@ -0,0 +1,528 @@
+(* TEST
+*)
+
+open Printf
+
+(* This is the module type of [Float.Array] except type [t] is abstract. *)
+module type S = sig
+  type t
+  val length : t -> int
+  val get : t -> int -> float
+  val set : t -> int -> float -> unit
+  val make : int -> float -> t
+  val create : int -> t
+  val init : int -> (int -> float) -> t
+  val append : t -> t -> t
+  val concat : t list -> t
+  val sub : t -> int -> int -> t
+  val copy : t -> t
+  val fill : t -> int -> int -> float -> unit
+  val blit : t -> int -> t -> int -> int -> unit
+  val to_list : t -> float list
+  val of_list : float list -> t
+  val iter : (float -> unit) -> t -> unit
+  val iteri : (int -> float -> unit) -> t -> unit
+  val map : (float -> float) -> t -> t
+  val mapi : (int -> float -> float) -> t -> t
+  val fold_left : ('a -> float -> 'a) -> 'a -> t -> 'a
+  val fold_right : (float -> 'a -> 'a) -> t -> 'a -> 'a
+  val iter2 : (float -> float -> unit) -> t -> t -> unit
+  val map2 : (float -> float -> float) -> t -> t -> t
+  val for_all : (float -> bool) -> t -> bool
+  val exists : (float -> bool) -> t -> bool
+  val mem : float -> t -> bool
+  val mem_ieee : float -> t -> bool
+  val sort : (float -> float -> int) -> t -> unit
+  val stable_sort : (float -> float -> int) -> t -> unit
+  val fast_sort : (float -> float -> int) -> t -> unit
+  val to_seq : t -> float Seq.t
+  val to_seqi : t -> (int * float) Seq.t
+  val of_seq : float Seq.t -> t
+  val map_to_array : (float -> 'a) -> t -> 'a array
+  val map_from_array : ('a -> float) -> 'a array -> t
+  val unsafe_get : t -> int -> float
+  val unsafe_set : t -> int -> float -> unit
+end
+
+(* module [Array] specialized to [float] and with a few changes,
+   satisfies signature S *)
+module Float_array : S = struct
+  include Stdlib.Array
+  let create = create_float
+  let map_to_array f a = map f a
+  let map_from_array f a = map f a
+  let mem_ieee x a = exists ((=) x) a
+  type t = float array
+end
+
+module Test (A : S) : sig end = struct
+
+  (* auxiliary functions *)
+
+  let neg_zero = 1.0 /. neg_infinity in
+
+  let rec check_i_upto a i =
+    if i >= 0 then begin
+      assert (A.get a i = Float.of_int i);
+      check_i_upto a (i - 1);
+    end
+  in
+
+  let check_i a = check_i_upto a (A.length a - 1) in
+
+  let check_inval f arg =
+    match f arg with
+    | _ -> assert false
+    | exception (Invalid_argument _) -> ()
+    | exception _ -> assert false
+  in
+
+  (* [make] [set] [get] *)
+  let a = A.make 1000 1.0 in
+  for i = 0 to 499 do A.set a i (Float.of_int i) done;
+  let rec loop i =
+    if i >= 0 then begin
+      assert (A.get a i = (if i < 500 then Float.of_int i else 1.0));
+      loop (i - 1);
+    end
+  in loop 999;
+  check_inval (A.get a) (-1);
+  check_inval (A.get a) (1000);
+  check_inval (fun i -> A.set a i 1.0) (-1);
+  check_inval (fun i -> A.set a i 1.0) 1000;
+  check_inval A.create (-1);
+  check_inval A.create (Sys.max_floatarray_length + 1);
+  check_inval (fun i -> A.make i 1.0) (-1);
+  check_inval (fun i -> A.make i 1.0) (Sys.max_floatarray_length + 1);
+
+  (* [length] *)
+  let test_length l = assert (l = (A.length (A.create l))) in
+  test_length 0;
+  test_length 10;
+  test_length 25;
+  test_length 255;
+  test_length 256;
+  test_length 1000;
+  test_length 123456;
+
+  (* [init] *)
+  let a = A.init 1000 Float.of_int in
+  check_i a;
+  check_inval (fun i -> A.init i Float.of_int) (-1);
+  check_inval (fun i -> A.init i Float.of_int) (Sys.max_floatarray_length + 1);
+
+  (* [append] *)
+  let check m n =
+    let a = A.init m Float.of_int in
+    let b = A.init n (fun x -> Float.of_int (x + m)) in
+    let c = A.append a b in
+    assert (A.length c = (m + n));
+    check_i c;
+  in
+  check 0 0;
+  check 0 100;
+  check 1 100;
+  check 100 0;
+  check 100 1;
+  check 100 100;
+  check 1000 1000;
+  (* check_inval omitted *)
+
+  (* [concat] *)
+  let check l =
+    let f (len, acc) n =
+      (len + n, A.init n (fun i -> Float.of_int (len + i)) :: acc)
+    in
+    let (total, ll) = List.fold_left f (0, []) l in
+    let b = A.concat (List.rev ll) in
+    assert (A.length b = total);
+    check_i b;
+  in
+  check [0; 0; 0];
+  check [1; 10; 100];
+  check [10; 0];
+  check [0];
+  check [1000; 1000; 1000];
+  check [];
+  (* check_inval omitted *)
+
+  (* [sub] *)
+  let a = A.init 1000 (fun i -> Float.of_int (i - 100)) in
+  let b = A.sub a 100 200 in
+  check_i b;
+  assert (A.length b = 200);
+  let b = A.sub a 1000 0 in
+  check_i (A.sub a 1000 0);
+  assert  (A.length b = 0);
+  check_inval (A.sub a (-1)) 0;
+  check_inval (A.sub a 0) (-1);
+  check_inval (A.sub a 0) 1001;
+  check_inval (A.sub a 1000) 1;
+
+  (* [copy] *)
+  let check len =
+    let a = A.init len Float.of_int in
+    let b = A.copy a in
+    check_i b;
+    assert (A.length b = len);
+  in
+  check 0;
+  check 1;
+  check 128;
+  check 1023;
+
+  (* [blit] [fill] *)
+  let test_blit_fill data initval ofs len =
+    let a = A.of_list data in
+    let b = A.create (List.length data) in
+    A.blit a 0 b 0 (A.length b);
+    assert (a = b);
+    A.fill b ofs len initval;
+    let rec check i = function
+      | [] -> ()
+      | hd :: tl ->
+          assert (A.get b i = (if i >= ofs && i < ofs + len
+                               then initval else hd));
+          check (i + 1) tl;
+    in
+    check 0 data
+  in
+  test_blit_fill [1.0;2.0;5.0;8.123;-100.456;212e19] 3.1415 3 2;
+  let a = A.create 100 in
+  check_inval (A.fill a (-1) 0) 1.0;
+  check_inval (A.fill a 0 (-1)) 1.0;
+  check_inval (A.fill a 0 101) 1.0;
+  check_inval (A.fill a 100 1) 1.0;
+  check_inval (A.fill a 101 0) 1.0;
+  check_inval (A.blit a (-1) a 0) 0;
+  check_inval (A.blit a 0 a 0) (-1);
+  check_inval (A.blit a 0 a 0) 101;
+  check_inval (A.blit a 100 a 0) 1;
+  check_inval (A.blit a 101 a 0) 0;
+  check_inval (A.blit a 0 a (-1)) 0;
+  check_inval (A.blit a 0 a 100) 1;
+  check_inval (A.blit a 0 a 101) 0;
+
+  (* [to_list] [of_list] *)
+  let a = A.create 1000 in
+  assert (A.of_list (A.to_list a) = a);
+  let a = A.create 0 in
+  assert (A.of_list (A.to_list a) = a);
+  (* check_inval omitted *)
+
+  (* [iter] *)
+  let a = A.init 300 (Float.of_int) in
+  let r = ref 0.0 in
+  A.iter (fun x -> assert (x = !r); r := x +. 1.0) a;
+  A.iter (fun _ -> assert false) (A.create 0);
+  assert (!r = 300.0);
+
+  (* [iteri] *)
+  let a = A.init 300 Float.of_int in
+  let r = ref 0 in
+  let f i x =
+    assert (i = !r);
+    assert (x = Float.of_int i);
+    r := i + 1
+  in
+  A.iteri f a;
+  A.iteri (fun _ _ -> assert false) (A.create 0);
+  assert (!r = 300);
+
+  (* [map], test result and order of evaluation *)
+  let a = A.init 500 Float.of_int in
+  let r = ref 0.0 in
+  let f x =
+    assert (x = !r);
+    r := !r +. 1.0;
+    x -. 1.0
+  in
+  let b = A.map f a in
+  check_i (A.sub b 1 499);
+
+  (* [mapi], test result and order of evaluation *)
+  let a = A.init 500 Float.of_int in
+  let r = ref 0.0 in
+  let f i x =
+    assert (x = Float.of_int i);
+    assert (x = !r);
+    r := !r +. 1.0;
+    x -. 1.0
+  in
+  let b = A.mapi f a in
+  check_i (A.sub b 1 499);
+
+  (* [fold_left], test result and order of evaluation *)
+  let a = A.init 500 Float.of_int in
+  let f acc x =
+    assert (acc = x);
+    x +. 1.0
+  in
+  let acc = A.fold_left f 0.0 a in
+  assert (acc = 500.0);
+
+  (* [fold_right], test result and order of evaluation *)
+  let a = A.init 500 Float.of_int in
+  let f x acc =
+    assert (x = acc -. 1.0);
+    x
+  in
+  let acc = A.fold_right f a 500.0 in
+  assert (acc = 0.0);
+
+  (* [iter2], test result and order of evaluation *)
+  let a = A.init 123 Float.of_int in
+  let b = A.init 123 Float.of_int in
+  let r = ref 0.0 in
+  let f x y =
+    assert (x = !r);
+    assert (y = !r);
+    r := !r +. 1.0;
+  in
+  A.iter2 f a b;
+  let c = A.create 456 in
+  check_inval (A.iter2 (fun _ _ -> assert false) a) c;
+  check_inval (A.iter2 (fun _ _ -> assert false) c) a;
+
+  (* [map2], test result and order of evaluation *)
+  let a = A.init 456 Float.of_int in
+  let b = A.init 456 (fun i -> Float.of_int i /. 2.0) in
+  let r = ref 0.0 in
+  let f x y =
+    assert (x = !r);
+    assert (y = !r /. 2.0);
+    r := !r +. 1.0;
+    2.0 *. (x -. y)
+  in
+  let c = A.map2 f a b in
+  check_i c;
+  let d = A.create 455 in
+  check_inval (A.map2 (fun _ _ -> assert false) a) d;
+  check_inval (A.map2 (fun _ _ -> assert false) d) a;
+
+  (* [for_all], test result and order of evaluation *)
+  let a = A.init 777 Float.of_int in
+  let r = ref 0.0 in
+  let f x =
+    assert (x = !r);
+    r := x +. 1.0;
+    true
+  in
+  assert (A.for_all f a);
+  let f x = assert (x = 0.0); false in
+  assert (not (A.for_all f a));
+
+  (* [exists], test result and order of evaluation *)
+  let a = A.init 777 Float.of_int in
+  let r = ref 0.0 in
+  let f x =
+    assert (x = !r);
+    r := x +. 1.0;
+    false
+  in
+  assert (not (A.exists f a));
+  let f x = assert (x = 0.0); true in
+  assert (A.exists f a);
+
+  (* [mem] *)
+  let a = A.init 7777 Float.of_int in
+  assert (A.mem 0.0 a);
+  assert (A.mem 7776.0 a);
+  assert (not (A.mem (-1.0) a));
+  assert (not (A.mem 7777.0 a));
+  let check v =
+    A.set a 1000 v;
+    assert (A.mem v a);
+  in
+  List.iter check [infinity; neg_infinity; neg_zero; nan];
+
+  (* [mem_ieee] *)
+  let a = A.init 7777 Float.of_int in
+  assert (A.mem_ieee 0.0 a);
+  assert (A.mem_ieee 7776.0 a);
+  assert (not (A.mem_ieee (-1.0) a));
+  assert (not (A.mem_ieee 7777.0 a));
+  let check v =
+    A.set a 1000 v;
+    assert (A.mem_ieee v a);
+  in
+  List.iter check [infinity; neg_infinity; neg_zero];
+  A.set a 0 nan;
+  assert (not (A.mem_ieee nan a));
+
+  (* [sort] [fast_sort] [stable_sort] *)
+  let check_sort sort cmp a =
+    let rec check_sorted a i =
+      if i + 1 < A.length a then begin
+        assert (cmp (A.get a i) (A.get a (i + 1)) <= 0);
+        check_sorted a (i + 1);
+      end
+    in
+    let rec check_permutation a b i =
+      let p = Array.make (A.length a) true in
+      let rec find lo hi x =
+        assert (lo < hi);
+        if hi = lo + 1 then begin
+          assert (cmp (A.get a lo) x = 0);
+          assert (p.(lo));
+          p.(lo) <- false;
+        end else begin
+          let mid = (lo + hi) / 2 in
+          assert (lo < mid && mid < hi);
+          match cmp (A.get a (mid - 1)) x with
+          | 0 when p.(mid - 1) -> find lo mid x
+          | 0 -> find mid hi x
+          | c when c < 0 -> find mid hi x
+          | c when c > 0 -> find lo mid x
+          | _ -> assert false
+        end
+      in
+      A.iter (find 0 (A.length a)) b
+    in
+    let b = A.copy a in
+    sort cmp a;
+    check_sorted a 0;
+    check_permutation a b 0;
+  in
+  Random.init 123;
+  let rand_float _ =
+    match Random.int 1004 with
+    | 1000 -> nan
+    | 1001 -> infinity
+    | 1002 -> neg_infinity
+    | 1003 -> neg_zero
+    | n when n < 500 -> Random.float 1.0
+    | _ -> -. Random.float 1.0
+  in
+  let check s =
+    let a = A.init 5 Float.of_int in
+    check_sort s Stdlib.compare a; (* already sorted *)
+    check_sort s (fun x y -> Stdlib.compare y x) a; (* reverse-sorted *)
+
+    let a = A.of_list [nan; neg_infinity; neg_zero; 0.; infinity] in
+    check_sort s Stdlib.compare a; (* already sorted *)
+    check_sort s (fun x y -> Stdlib.compare y x) a; (* reverse-sorted *)
+
+    let a = A.init 50000 rand_float in
+    check_sort s Stdlib.compare a;
+    let a = A.make 1000 1.0 in
+    check_sort s Stdlib.compare a;
+    let a = A.append (A.make 1000 1.0) (A.make 1000 2.0) in
+    check_sort s Stdlib.compare a;
+  in
+  check A.sort;
+  check A.stable_sort;
+  check A.fast_sort;
+
+  (* [to_seq] *)
+  let check_seq a =
+    let r = ref 0 in
+    let f x =
+      assert (A.get a !r = x);
+      r := !r + 1;
+    in
+    let s = A.to_seq a in
+    Seq.iter f s;
+  in
+  check_seq (A.init 999 Float.of_int);
+  check_seq (A.create 0);
+
+  (* [to_seqi] *)
+  let check_seqi a =
+    let r = ref 0 in
+    let f (i, x) =
+      assert (i = !r);
+      assert (A.get a !r = x);
+      r := !r + 1;
+    in
+    let s = A.to_seqi a in
+    Seq.iter f s;
+  in
+  check_seqi (A.init 999 Float.of_int);
+  check_seqi (A.create 0);
+
+  (* [of_seq] *)
+  let r = ref 0 in
+  let rec f () =
+    if !r = 100 then Seq.Nil else begin
+      let res = Seq.Cons (Float.of_int !r, f) in
+      r := !r + 1;
+      res
+    end
+  in
+  let a = A.of_seq f in
+  assert (a = A.init 100 Float.of_int);
+  assert (A.of_seq Seq.empty = A.create 0);
+
+  (* [map_to_array] *)
+  let r = ref 0 in
+  let f x =
+    assert (x = Float.of_int !r);
+    r := !r + 1;
+    x *. 2.0
+  in
+  let a = A.init 876 Float.of_int in
+  let ar1 = A.map_to_array f a in
+  let ar2 = Array.init 876 (fun x -> Float.of_int (2 * x)) in
+  assert (ar1 = ar2);
+  let ar = A.map_to_array (fun _ -> assert false) (A.create 0) in
+  assert (ar = [| |]);
+
+  (* [map_from_array] *)
+  let r = ref 0 in
+  let f x =
+    assert (x = Float.of_int !r);
+    r := !r + 1;
+    x *. 2.0
+  in
+  let ar = Array.init 876 Float.of_int in
+  let a1 = A.map_from_array f ar in
+  let a2 = A.init 876 (fun x -> Float.of_int (2 * x)) in
+  assert (a1 = a2);
+  let a = A.map_from_array (fun _ -> assert false) [| |] in
+  assert (a = A.create 0);
+
+  (* comparisons *)
+  let normalize_comparison n =
+    if n = 0 then 0 else if n < 0 then -1 else 1
+  in
+  let check c l1 l2 =
+    assert (c = (normalize_comparison (compare (A.of_list l1) (A.of_list l2))))
+  in
+  check 0    [0.0; 0.25; -4.0; 3.141592654; nan]
+             [0.0; 0.25; -4.0; 3.141592654; nan];
+  check (-1) [0.0; 0.25; nan]
+             [0.0; 0.25; 3.14];
+  check (-1) [0.0; 0.25; -4.0]
+             [0.0; 0.25; 3.14159];
+  check 1    [0.0; 2.718; -4.0]
+             [0.0; 0.25; 3.14159];
+  check 1    [0.0; 2.718; -4.0]
+             [nan; 0.25; 3.14159];
+
+  (* [unsafe_get] [unsafe_set] *)
+  let a = A.create 3 in
+  for i = 0 to 2 do A.unsafe_set a i (float i) done;
+  for i = 0 to 2 do assert (A.unsafe_get a i = float i) done;
+
+  (* I/O *)
+  let test_structured_io value =
+    let (tmp, oc) =
+      Filename.open_temp_file ~mode:[Open_binary] "floatarray" ".data"
+    in
+    Marshal.to_channel oc value [];
+    close_out oc;
+    let ic = open_in_bin tmp in
+    let value' = Marshal.from_channel ic in
+    close_in ic;
+    Sys.remove tmp;
+    assert (compare value value' = 0)
+  in
+  let l = [0.; 0.25; -4.; 3.14159265; nan; infinity; neg_infinity; neg_zero] in
+  test_structured_io (A.of_list l);
+
+end
+
+(* We run the same tests on [Float.Array] and [Array]. *)
+module T1 = Test (Stdlib.Float.Array)
+module T2 = Test (Float_array)
diff --git a/testsuite/tests/lib-floatarray/ocamltests b/testsuite/tests/lib-floatarray/ocamltests
new file mode 100644
index 000000000..abf190238
--- /dev/null
+++ b/testsuite/tests/lib-floatarray/ocamltests
@@ -0,0 +1 @@
+floatarray.ml