docs: Add refman+stdlib documentation

10 files changed, 242 insertions, 26 deletions

diff --git a/doc/sphinx/language/coq-library.rst b/doc/sphinx/language/coq-library.rst
index ac75240cfb..cad5e4e67e 100644
--- a/doc/sphinx/language/coq-library.rst
+++ b/doc/sphinx/language/coq-library.rst
@@ -756,6 +756,7 @@ subdirectories:
   * **Sets** : Sets (classical, constructive, finite, infinite, power set, etc.) 
   * **FSets** : Specification and implementations of finite sets and finite maps (by lists and by AVL trees)
   * **Reals** : Axiomatization of real numbers (classical, basic functions, integer part, fractional part, limit, derivative, Cauchy series, power series and results,...)
+  * **Floats** : Machine implementation of floating-point arithmetic (for the binary64 format)
   * **Relations** : Relations (definitions and basic results)
   * **Sorting** : Sorted list (basic definitions and heapsort correctness)
   * **Strings** : 8-bits characters and strings
@@ -768,7 +769,7 @@ are directly accessible with the command ``Require`` (see
 Section :ref:`compiled-files`).
 
 The different modules of the |Coq| standard library are documented
-online at http://coq.inria.fr/stdlib.
+online at https://coq.inria.fr/stdlib.
 
 Peano’s arithmetic (nat)
 ~~~~~~~~~~~~~~~~~~~~~~~~
@@ -988,6 +989,106 @@ Notation    Interpretation  Precedence  Associativity
 ``_ :: _``  ``cons``        60          right
 ==========  ==============  ==========  =============
 
+.. _floats_library:
+
+Floats library
+~~~~~~~~~~~~~~
+
+The library of primitive floating-point arithmetic can be loaded by
+requiring module ``Floats``:
+
+.. coqtop:: in
+
+  Require Import Floats.
+
+It exports the module ``PrimFloat`` that provides a primitive type
+named ``float``, defined in the kernel (see section :ref:`primitive-floats`),
+as well as two variant types ``float_comparison`` and ``float_class``:
+
+
+.. coqtop:: all
+
+  Print float.
+  Print float_comparison.
+  Print float_class.
+
+It then defines the primitive operators below, using the processor
+floating-point operators for binary64 in rounding-to-nearest even:
+
+* ``abs``
+* ``opp``
+* ``sub``
+* ``add``
+* ``mul``
+* ``div``
+* ``sqrt``
+* ``compare`` : compare two floats and return a ``float_comparison``
+* ``classify`` : analyze a float and return a ``float_class``
+* ``of_int63`` : round a primitive integer and convert it into a float
+* ``normfr_mantissa`` : take a float in ``[0.5; 1.0)`` and return its mantissa
+* ``frshiftexp`` : convert a float to fractional part in ``[0.5; 1.0)`` and integer part
+* ``ldshiftexp`` : multiply a float by an integral power of ``2``
+* ``next_up`` : return the next float towards positive infinity
+* ``next_down`` : return the next float towards negative infinity
+
+For special floating-point values, the following constants are also
+defined:
+
+* ``zero``
+* ``neg_zero``
+* ``one``
+* ``two``
+* ``infinity``
+* ``neg_infinity``
+* ``nan`` : Not a Number (assumed to be unique: the "payload" of NaNs is ignored)
+
+The following table shows the notations available when opening scope
+``float_scope``.
+
+===========  ==============
+Notation     Interpretation
+===========  ==============
+``- _``      ``opp``
+``_ - _``    ``sub``
+``_ + _``    ``add``
+``_ * _``    ``mul``
+``_ / _``    ``div``
+``_ == _``   ``eqb``
+``_ < _``    ``ltb``
+``_ <= _``   ``leb``
+``_ ?= _``   ``compare``
+===========  ==============
+
+Floating-point constants are parsed and pretty-printed as (17-digit)
+decimal constants. This ensures that the composition
+:math:`\text{parse} \circ \text{print}` amounts to the identity.
+
+.. example::
+
+  .. coqtop:: all
+
+    Open Scope float_scope.
+    Eval compute in 1 + 0.5.
+    Eval compute in 1 / 0.
+    Eval compute in 1 / -0.
+    Eval compute in 0 / 0.
+    Eval compute in 0 ?= -0.
+    Eval compute in nan ?= nan.
+    Eval compute in next_down (-1).
+
+The primitive operators are specified with respect to their Gallina
+counterpart, using the variant type ``spec_float``, and the injection
+``Prim2SF``:
+
+.. coqtop:: all
+
+  Print spec_float.
+  Check Prim2SF.
+  Check mul_spec.
+
+For more details on the available definitions and lemmas, see the
+online documentation of the ``Floats`` library.
+
 .. _userscontributions:
 
 Users’ contributions
diff --git a/doc/sphinx/language/gallina-extensions.rst b/doc/sphinx/language/gallina-extensions.rst
index f50cf9340c..3f122f122b 100644
--- a/doc/sphinx/language/gallina-extensions.rst
+++ b/doc/sphinx/language/gallina-extensions.rst
@@ -2411,7 +2411,7 @@ by means of the interactive proof engine.
 .. _primitive-integers:
 
 Primitive Integers
---------------------------------
+------------------
 
 The language of terms features 63-bit machine integers as values. The type of
 such a value is *axiomatized*; it is declared through the following sentence
@@ -2464,6 +2464,48 @@ wrapped into the :g:`Uint63.of_int` (resp. :g:`Uint63.of_int64`) constructor on
 64-bit (resp. 32-bit) platforms. Currently, this cannot be customized (see the
 function :g:`Uint63.compile` from the kernel).
 
+.. _primitive-floats:
+
+Primitive Floats
+----------------
+
+The language of terms features Binary64 floating-point numbers as values.
+The type of such a value is *axiomatized*; it is declared through the
+following sentence (excerpt from the :g:`PrimFloat` module):
+
+.. coqdoc::
+
+   Primitive float := #float64_type.
+
+This type is equipped with a few operators, that must be similarly declared.
+For instance, the product of two primitive floats can be computed using the
+:g:`PrimFloat.mul` function, declared and specified as follows:
+
+.. coqdoc::
+
+   Primitive mul := #float64_mul.
+   Notation "x * y" := (mul x y) : float_scope.
+
+   Axiom mul_spec : forall x y, Prim2SF (x * y)%float = SF64mul (Prim2SF x) (Prim2SF y).
+
+where :g:`Prim2SF` is defined in the :g:`FloatOps` module.
+
+The set of such operators is described in section :ref:`floats_library`.
+
+These primitive declarations are regular axioms. As such, they must be trusted, and are listed by the
+:g:`Print Assumptions` command.
+
+The reduction machines (:tacn:`vm_compute`, :tacn:`native_compute`) implement
+dedicated, efficient rules to reduce the applications of these primitive
+operations, using the floating-point processor operators that are assumed
+to comply with the IEEE 754 standard for floating-point arithmetic.
+
+These primitives, when extracted to OCaml (see :ref:`extraction`), are mapped to
+types and functions of a :g:`Float64` module. Said module is not produced by
+extraction. Instead, it has to be provided by the user (if they want to compile
+or execute the extracted code). For instance, an implementation of this module
+can be taken from the kernel of Coq.
+
 Bidirectionality hints
 ----------------------
 
diff --git a/doc/stdlib/index-list.html.template b/doc/stdlib/index-list.html.template
index f0ada745e7..eedd8a3d61 100644
--- a/doc/stdlib/index-list.html.template
+++ b/doc/stdlib/index-list.html.template
@@ -328,6 +328,19 @@ through the <tt>Require Import</tt> command.</p>
     theories/Numbers/Integer/Binary/ZBinary.v
     theories/Numbers/Integer/NatPairs/ZNatPairs.v
     </dd>
+
+    <dt> <b>&nbsp;&nbsp;Floats</b>:
+       Floating-point arithmetic
+    </dt>
+    <dd>
+    theories/Floats/FloatClass.v
+    theories/Floats/PrimFloat.v
+    theories/Floats/SpecFloat.v
+    theories/Floats/FloatOps.v
+    theories/Floats/FloatAxioms.v
+    theories/Floats/FloatLemmas.v
+    (theories/Floats/Floats.v)
+    </dd>
     </dl>
   </dd>
 
diff --git a/theories/Floats/FloatAxioms.v b/theories/Floats/FloatAxioms.v
index d094032805..dc637e50a6 100644
--- a/theories/Floats/FloatAxioms.v
+++ b/theories/Floats/FloatAxioms.v
@@ -1,5 +1,7 @@
 Require Import ZArith Int63 SpecFloat PrimFloat FloatOps.
 
+(** * Properties of the primitive operators for the Binary64 format *)
+
 Notation valid_binary := (valid_binary prec emax).
 
 Definition SF64classify := SFclassify prec.
diff --git a/theories/Floats/FloatClass.v b/theories/Floats/FloatClass.v
index 3241035aef..627cb648f9 100644
--- a/theories/Floats/FloatClass.v
+++ b/theories/Floats/FloatClass.v
@@ -1,2 +1,2 @@
-Inductive float_class : Set :=
+Variant float_class : Set :=
   | PNormal | NNormal | PSubn | NSubn | PZero | NZero | PInf | NInf | NaN.
diff --git a/theories/Floats/FloatLemmas.v b/theories/Floats/FloatLemmas.v
index 7b3714ab77..4e1f14610d 100644
--- a/theories/Floats/FloatLemmas.v
+++ b/theories/Floats/FloatLemmas.v
@@ -1,6 +1,8 @@
 Require Import ZArith Int63 SpecFloat PrimFloat FloatOps FloatAxioms.
 Require Import Psatz.
 
+(** * Support results involving frexp and ldexp *)
+
 Lemma shift_value : [|shift|]%int63 = (2*emax + prec)%Z.
   reflexivity.
 Qed.
diff --git a/theories/Floats/FloatOps.v b/theories/Floats/FloatOps.v
index 6cc7cb0568..5ffbfc7215 100644
--- a/theories/Floats/FloatOps.v
+++ b/theories/Floats/FloatOps.v
@@ -1,6 +1,7 @@
 Require Import ZArith Int63 SpecFloat PrimFloat.
 
-(* Properties of the Binary64 IEEE 754 format *)
+(** * Derived operations and mapping between primitive [float]s and [spec_float]s *)
+
 Definition prec := 53%Z.
 Definition emax := 1024%Z.
 Notation emin := (emin prec emax).
@@ -15,6 +16,9 @@ Definition ldexp f e :=
 
 Definition ulp f := ldexp one (fexp prec emax (snd (frexp f))).
 
+(** [Prim2SF] is an injective function that will be useful to express
+the properties of the implemented Binary64 format (see [FloatAxioms]).
+*)
 Definition Prim2SF f :=
   if is_nan f then S754_nan
   else if is_zero f then S754_zero (get_sign f)
diff --git a/theories/Floats/Floats.v b/theories/Floats/Floats.v
index 53dd73f6d0..700c69b99d 100644
--- a/theories/Floats/Floats.v
+++ b/theories/Floats/Floats.v
@@ -1,6 +1,17 @@
+(** The Floats library is split in 6 theories:
+- FloatClass: define the [float_class] inductive
+- PrimFloat: define the floating-point values and operators as kernel primitives
+- SpecFloat: specify the floating-point operators with binary integers
+- FloatOps: define conversion functions between [spec_float] and [float]
+- FloatAxioms: state properties of the primitive operators w.r.t. [spec_float]
+- FloatLemmas: prove a few results involving frexp and ldexp
+
+For a brief overview of the Floats library,
+see {{https://coq.inria.fr/distrib/current/refman/language/coq-library.html#floats-library}} *)
+
 Require Export FloatClass.
-Require Export SpecFloat.
 Require Export PrimFloat.
+Require Export SpecFloat.
 Require Export FloatOps.
 Require Export FloatAxioms.
 Require Export FloatLemmas.
diff --git a/theories/Floats/PrimFloat.v b/theories/Floats/PrimFloat.v
index 24e4ac2da2..880252c2b9 100644
--- a/theories/Floats/PrimFloat.v
+++ b/theories/Floats/PrimFloat.v
@@ -1,23 +1,38 @@
 Require Import Int63 FloatClass.
 
-Inductive float_comparison : Set := FEq | FLt | FGt | FNotComparable.
+(** * Definition of the interface for primitive floating-point arithmetic
+
+This interface provides processor operators for the Binary64 format of the
+IEEE 754-2008 standard. *)
+
+(** ** Type definition for the co-domain of [compare] *)
+Variant float_comparison : Set := FEq | FLt | FGt | FNotComparable.
 
 Register float_comparison as kernel.ind_f_cmp.
+
 Register float_class as kernel.ind_f_class.
 
+(** ** The main type *)
+(** [float]: primitive type for Binary64 floating-point numbers. *)
 Primitive float := #float64_type.
 
+(** ** Syntax support *)
 Declare Scope float_scope.
 Delimit Scope float_scope with float.
 Bind Scope float_scope with float.
 
 Declare ML Module "float_syntax_plugin".
 
-Primitive opp := #float64_opp.
-Notation "- x" := (opp x) : float_scope.
+(** ** Floating-point operators *)
+Primitive classify := #float64_classify.
 
 Primitive abs := #float64_abs.
 
+Primitive sqrt := #float64_sqrt.
+
+Primitive opp := #float64_opp.
+Notation "- x" := (opp x) : float_scope.
+
 Primitive eqb := #float64_eq.
 Notation "x == y" := (eqb x y) (at level 70, no associativity) : float_scope.
 
@@ -30,43 +45,52 @@ Notation "x <= y" := (leb x y) (at level 70, no associativity) : float_scope.
 Primitive compare := #float64_compare.
 Notation "x ?= y" := (compare x y) (at level 70, no associativity) : float_scope.
 
-Primitive classify := #float64_classify.
-
 Primitive mul := #float64_mul.
 Notation "x * y" := (mul x y) : float_scope.
 
 Primitive add := #float64_add.
 Notation "x + y" := (add x y) : float_scope.
+
 Primitive sub := #float64_sub.
 Notation "x - y" := (sub x y) : float_scope.
 
 Primitive div := #float64_div.
 Notation "x / y" := (div x y) : float_scope.
 
-Primitive sqrt := #float64_sqrt.
+(** ** Conversions *)
 
-(* Convert a primitive integer into a float value.
- * The value is rounded if necessary. *)
+(** [of_int63]: convert a primitive integer into a float value.
+    The value is rounded if need be. *)
 Primitive of_int63 := #float64_of_int63.
 
-(* If the input is a float value with an absolute value
- * inside [0.5; 1.) then return its mantissa as a primitive integer.
- * The mantissa will be a 53-bit integer with its most significant bit set to 1.
- * Else return zero.
- * The sign bit is always ignored. *)
+(** Specification of [normfr_mantissa]:
+- If the input is a float value with an absolute value inside $[0.5, 1.)$#[0.5, 1.)#;
+- Then return its mantissa as a primitive integer.
+  The mantissa will be a 53-bit integer with its most significant bit set to 1;
+- Else return zero.
+
+The sign bit is always ignored. *)
 Primitive normfr_mantissa := #float64_normfr_mantissa.
 
-(* Exponent manipulation functions *)
-Definition shift := (2101)%int63. (* = 2*emax + prec *)
+(** ** Exponent manipulation functions *)
+Definition shift := 2101%int63. (** [= 2*emax + prec] *)
+
+(** [frshiftexp]: convert a float to fractional part in $[0.5, 1.)$#[0.5, 1.)#
+and integer part. *)
 Primitive frshiftexp := #float64_frshiftexp.
+
+(** [ldshiftexp]: multiply a float by an integral power of 2. *)
 Primitive ldshiftexp := #float64_ldshiftexp.
 
+(** ** Predecesor/Successor functions *)
+
+(** [next_up]: return the next float towards positive infinity. *)
 Primitive next_up := #float64_next_up.
-Primitive next_down := #float64_next_down.
 
-Local Open Scope float_scope.
+(** [next_down]: return the next float towards negative infinity. *)
+Primitive next_down := #float64_next_down.
 
-(* Special values, needed for pretty-printing *)
+(** ** Special values (needed for pretty-printing) *)
 Definition infinity := Eval compute in div (of_int63 1) (of_int63 0).
 Definition neg_infinity := Eval compute in opp infinity.
 Definition nan := Eval compute in div (of_int63 0) (of_int63 0).
@@ -75,12 +99,13 @@ Register infinity as num.float.infinity.
 Register neg_infinity as num.float.neg_infinity.
 Register nan as num.float.nan.
 
-(* Other special values *)
+(** ** Other special values *)
 Definition one := Eval compute in (of_int63 1).
 Definition zero := Eval compute in (of_int63 0).
-Definition neg_zero := Eval compute in (-zero).
+Definition neg_zero := Eval compute in (-zero)%float.
 Definition two := Eval compute in (of_int63 2).
 
+(** ** Predicates and helper functions *)
 Definition is_nan f := negb (f == f)%float.
 
 Definition is_zero f := (f == zero)%float. (* note: 0 == -0 with floats *)
@@ -89,6 +114,7 @@ Definition is_infinity f := (abs f == infinity)%float.
 
 Definition is_finite (x : float) := negb (is_nan x || is_infinity x).
 
+(** [get_sign]: return [true] for [-] sign, [false] for [+] sign. *)
 Definition get_sign f :=
   let f := if is_zero f then (one / f)%float else f in
   (f < zero)%float.
diff --git a/theories/Floats/SpecFloat.v b/theories/Floats/SpecFloat.v
index 8708782e39..fd0aa5e075 100644
--- a/theories/Floats/SpecFloat.v
+++ b/theories/Floats/SpecFloat.v
@@ -1,11 +1,25 @@
 Require Import ZArith FloatClass.
 
+(** * Specification of floating-point arithmetic
+
+This specification is mostly borrowed from the [IEEE754.Binary] module
+of the Flocq library (see {{http://flocq.gforge.inria.fr/}}) *)
+
+(** ** Inductive specification of floating-point numbers
+
+Similar to [Flocq.IEEE754.Binary.full_float], but with no NaN payload. *)
 Variant spec_float :=
   | S754_zero (s : bool)
   | S754_infinity (s : bool)
   | S754_nan
   | S754_finite (s : bool) (m : positive) (e : Z).
 
+(** ** Parameterized definitions
+
+[prec] is the number of bits of the mantissa including the implicit one;
+[emax] is the exponent of the infinities.
+
+For instance, Binary64 is defined by [prec = 53] and [emax = 1024]. *)
 Section FloatOps.
   Variable prec emax : Z.
 
@@ -132,7 +146,8 @@ Section FloatOps.
       end.
   End Rounding.
 
-  (* Define operations *)
+  (** ** Define operations *)
+
   Definition SFopp x :=
     match x with
     | S754_nan => S754_nan