Merge branch 'feature/kernel/arithmetic' into 'master'

Regroup all arithmetic modules in one place

See merge request frama-c/frama-c!4923

src/kernel_services/analysis/filter/linear_filter_test.ml

@@ -20,43 +20,6 @@
 (*                                                                        *)
 (**************************************************************************)
 
-module Rational = struct
-
-  type scalar = Q.t
-
-  module Type = struct
-    include Datatype.Serializable_undefined
-    type t = scalar
-    let name = "Linear.Filter.Test.Rational"
-    let reprs = [ Q.zero ]
-    let compare = Q.compare
-    let equal = Q.equal
-    let hash q = Z.hash (Q.num q) + 11 * Z.hash (Q.den q)
-  end
-
-  include Datatype.Make_with_collections (Type)
-
-  let pretty fmt n =
-    let ten = Z.of_int 10 in
-    let sign = if Q.sign n >= 0 then "" else "-" in
-    let num = Q.num n |> Z.abs and den = Q.den n |> Z.abs in
-    let finish n = Z.Compare.(n >= den || n = Z.zero) in
-    let rec f e n = if finish n then (n, e) else f (e + 1) Z.(n * ten) in
-    let num, exponent = f 0 num in
-    let default fmt n = Format.fprintf fmt "%1.7f" (Q.to_float n) in
-    if exponent > 0 then
-      let number = Q.make num den in
-      Format.fprintf fmt "%s%aE-%d" sign default number exponent
-    else default fmt n
-
-  include Q
-  let infinity = Q.inf
-  let ( = ) = Q.equal
-
-end
-
-
-
 module Filter = Linear_filter.Make (Rational)
 module Linear = Filter.Linear
 

src/libraries/arithmetic/field.ml

+(**************************************************************************)
+(*                                                                        *)
+(*  This file is part of Frama-C.                                         *)
+(*                                                                        *)
+(*  Copyright (C) 2007-2025                                               *)
+(*    CEA (Commissariat à l'énergie atomique et aux énergies              *)
+(*         alternatives)                                                  *)
+(*                                                                        *)
+(*  you can redistribute it and/or modify it under the terms of the GNU   *)
+(*  Lesser General Public License as published by the Free Software       *)
+(*  Foundation, version 2.1.                                              *)
+(*                                                                        *)
+(*  It is distributed in the hope that it will be useful,                 *)
+(*  but WITHOUT ANY WARRANTY; without even the implied warranty of        *)
+(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *)
+(*  GNU Lesser General Public License for more details.                   *)
+(*                                                                        *)
+(*  See the GNU Lesser General Public License version 2.1                 *)
+(*  for more details (enclosed in the file licenses/LGPLv2.1).            *)
+(*                                                                        *)
+(**************************************************************************)
+
+(** This module provides a generic signature describing mathematical fields,
+    i.e algebraic structure that behave like rationals, reals or complex
+    numbers. The signature also provides several functions that are not
+    directly part of fields definition, but are useful nonetheless, in
+    particular when using fields to model floating point computations.
+    @since Frama-C+dev *)
+
+(** This type is used to return bounds for mathematical functions that
+    cannot be exactly computed using rationals, like the square root
+    for instance. *)
+type 't bounds = { lower : 't ; upper : 't }
+
+module type S = sig
+
+  include Datatype.S_with_collections
+  type scalar = t
+
+  (** {2 Useful constants in the field.}
+
+      Two is provided for users that model floating point computations using
+      abstractions over the field. *)
+
+  val zero    : scalar
+  val one     : scalar
+  val two     : scalar
+  val pos_inf : scalar
+  val neg_inf : scalar
+
+  (** {2 Constructors interacting with standard OCaml types.} *)
+
+  val of_int    : int -> scalar
+  val of_float  : float -> scalar
+  val to_float  : scalar -> float
+  val of_string : string -> scalar
+
+  (** {2 Standard arithmetic unary operations.} *)
+
+  val neg : scalar -> scalar
+  val abs : scalar -> scalar
+
+  (** {2 Standard arithmetic binary operations.} *)
+
+  val ( + ) : scalar -> scalar -> scalar
+  val ( - ) : scalar -> scalar -> scalar
+  val ( * ) : scalar -> scalar -> scalar
+  val ( / ) : scalar -> scalar -> scalar
+
+  (** {2 Standard arithmetic relational operations.} *)
+
+  val ( =  ) : scalar -> scalar -> bool
+  val ( != ) : scalar -> scalar -> bool
+  val ( <= ) : scalar -> scalar -> bool
+  val ( <  ) : scalar -> scalar -> bool
+  val ( >= ) : scalar -> scalar -> bool
+  val ( >  ) : scalar -> scalar -> bool
+
+  (** {2 Other arithmetic operations. }*)
+
+  (** [pow2 e] computes the scalar [2 ^ e]. *)
+  val pow2 : int -> scalar
+
+  (** [log2 f] computes the integer [l] and [u] such as
+      [l = floor(log₂ f)] and [u = ceiling(log₂ f)]. *)
+  val log2 : scalar -> int bounds
+
+  (** [sqrt f] computes bounds of the square root of [f]. *)
+  val sqrt : scalar -> scalar bounds
+
+  val max : scalar -> scalar -> scalar
+  val min : scalar -> scalar -> scalar
+
+end

src/kernel_services/analysis/filter/finite.mli → src/libraries/arithmetic/finite.mli

@@ -22,11 +22,12 @@
 
 open Nat
 
+(** Encoding of finite set in OCaml type system. *)
 
-
-(* The type [n finite] encodes all finite sets of cardinal [n]. It is used by
-   the module Linear to represent accesses to vectors and matrices coefficients,
-   statically ensuring that no out of bounds access can be performed. *)
+(** The type [n finite] encodes all finite sets of cardinal [n]. It is used by
+    the module {!Linear} to represent accesses to vectors and matrices
+    coefficients, statically ensuring that no out of bounds access can be
+    performed. *)
 type 'n finite
 
 val first : 'n succ finite
@@ -34,15 +35,16 @@ val last  : 'n succ nat -> 'n succ finite
 val next  : 'n finite -> 'n succ finite
 val ( = ) : 'n finite -> 'n finite -> bool
 
-(* The call [of_int limit n] returns a finite value representing the n-nd
-   element of a finite set of cardinal limit. If n is not in the bounds, none is
-   returned. This function complexity is O(1). *)
+(** The call [of_int limit n] returns a finite value representing the n-nd
+    element of a finite set of cardinal limit. If n is not in the bounds, [None]
+    is returned. This function complexity is O(1). *)
 val of_int : 'n succ nat -> int -> 'n succ finite option
 
-(* The call [to_int n] returns an integer equal to n. This function complexity
-   is O(1). *)
+(** The call [to_int n] returns an integer equal to n. This function complexity
+    is O(1). *)
 val to_int : 'n finite -> int
 
-(* The call [for_each acc limit f] folds over each finite elements of a set of
-   cardinal limit, computing f at each step. The function complexity is O(n). *)
+(** The call [for_each f limit acc] folds over each finite elements of a set of
+    cardinal limit, computing f at each step.
+    The function complexity is O(n). *)
 val for_each : ('n finite -> 'a -> 'a) -> 'n nat -> 'a -> 'a

src/kernel_services/analysis/filter/linear.mli → src/libraries/arithmetic/linear.mli

@@ -25,19 +25,19 @@ open Finite
 
 
 
-(* Definition of a linear space over a field. Used by Linear_filter to represent
-   and compute linear filters invariants. *)
+(** Definition of a linear space over a field. Used by {!Linear_filter} to
+    represent and compute linear filters invariants. *)
 
 module Space (Field : Field.S) : sig
 
-  (* The type of scalars in the field 𝕂. *)
+  (** The type of scalars in the field 𝕂. *)
   type scalar = Field.scalar
 
-  (* The type of matrices in 𝕂ⁿˣᵐ *)
+  (** The type of matrices in 𝕂ⁿˣᵐ *)
   type ('n, 'm) matrix
 
-  (* Type representing a column vector of 𝕂ⁿ. One can use [transpose] if one
-     needs a row vector, for example when printing it. *)
+  (** Type representing a column vector of 𝕂ⁿ. One can use {!Matrix.transpose}
+      if one needs a row vector, for example when printing it. *)
   type 'n vector = ('n, zero succ) matrix
 
 
@@ -46,28 +46,28 @@ module Space (Field : Field.S) : sig
 
     val pretty : Format.formatter -> 'n vector -> unit
 
-    (* The call [zero n] returns the 0 vector in 𝕂ⁿ. *)
+    (** The call [zero n] returns the 0 vector in 𝕂ⁿ. *)
     val zero   : 'n succ nat -> 'n succ vector
 
-    (* The call [repeat x n] returns a vector in 𝕂ⁿ which each dimension
-       containing the scalar x. *)
+    (** The call [repeat x n] returns a vector in 𝕂ⁿ which each dimension
+        containing the scalar x. *)
     val repeat : scalar -> 'n succ nat -> 'n succ vector
 
-    (* The call [base i n] returns the i-th base vector in the orthonormal space
-       of 𝕂ⁿ. In other words, the returned vector contains zero except for the
-       i-th coordinate, which contains one. *)
+    (** The call [base i n] returns the i-th base vector in the orthonormal
+        space of 𝕂ⁿ. In other words, the returned vector contains zero except
+        for the i-th coordinate, which contains one. *)
     val base   : 'n succ finite -> 'n succ nat -> 'n succ vector
 
-    (* The call [set i x v] returns a new vector of the same linear space as
-       [v], and with the same coordinates, except for the i-th one, which is
-       set to the scalar [x]. *)
+    (** The call [set i x v] returns a new vector of the same linear space as
+        [v], and with the same coordinates, except for the i-th one, which is
+        set to the scalar [x]. *)
     val set    : 'n finite -> scalar -> 'n vector -> 'n vector
 
-    (* The call [size v] for [v] a vector of 𝕂ⁿ returns n. *)
+    (** The call [size v] for [v] a vector of 𝕂ⁿ returns n. *)
     val size   : 'n vector -> 'n nat
 
-    (* The call [norm v] computes the ∞-norm of [v], i.e the maximum of the
-       absolute values of its coordinates. *)
+    (** The call [norm v] computes the ∞-norm of [v], i.e the maximum of the
+        absolute values of its coordinates. *)
     val norm   : 'n vector -> scalar
 
   end
@@ -78,43 +78,43 @@ module Space (Field : Field.S) : sig
 
     val pretty : Format.formatter -> ('n, 'm) matrix -> unit
 
-    (* The call [id n] returns the identity matrix in 𝕂ⁿˣⁿ. *)
+    (** The call [id n] returns the identity matrix in 𝕂ⁿˣⁿ. *)
     val id : 'n succ nat -> ('n succ, 'n succ) matrix
 
-    (* The call [zero n m] returns the 0 matrix in 𝕂ⁿˣᵐ. *)
+    (** The call [zero n m] returns the 0 matrix in 𝕂ⁿˣᵐ. *)
     val zero : 'n succ nat -> 'm succ nat -> ('n succ, 'm succ) matrix
 
-    (* The call [get i j m] returns the coefficient of the i-th row and
-       the j-th column. *)
+    (** The call [get i j m] returns the coefficient of the i-th row and
+        the j-th column. *)
     val get : 'n finite -> 'm finite -> ('n, 'm) matrix -> scalar
 
-    (* The call [set i j x m] returns a new matrix of the same linear space as
-       [m], and with the same coefficients, except for the one of the i-th row
-       and the j-th column, which is set to the scalar [x]. *)
+    (** The call [set i j x m] returns a new matrix of the same linear space as
+        [m], and with the same coefficients, except for the one of the i-th row
+        and the j-th column, which is set to the scalar [x]. *)
     val set : 'n finite -> 'm finite -> scalar -> ('n, 'm) matrix -> ('n, 'm) matrix
 
-    (* The call [norm m] computes the ∞-norm of [m], i.e the maximum of the
-       absolute sums of the rows of [m]. *)
+    (** The call [norm m] computes the ∞-norm of [m], i.e the maximum of the
+        absolute sums of the rows of [m]. *)
     val norm : ('n, 'm) matrix -> scalar
 
-    (* The call [transpose m] for m in 𝕂ⁿˣᵐ returns a new matrix in 𝕂ᵐˣⁿ with
-       all the coefficients transposed. *)
+    (** The call [transpose m] for m in 𝕂ⁿˣᵐ returns a new matrix in 𝕂ᵐˣⁿ with
+        all the coefficients transposed. *)
     val transpose : ('n, 'm) matrix -> ('m, 'n) matrix
 
-    (* The call [dimensions m] for m in 𝕂ⁿˣᵐ returns the pair (n, m). *)
+    (** The call [dimensions m] for m in 𝕂ⁿˣᵐ returns the pair (n, m). *)
     val dimensions : ('m, 'n) matrix -> 'm nat * 'n nat
 
-    (* Matrices addition. The dimensions compatibility is statically ensured. *)
+    (** Matrices addition. The dimensions compatibility is statically ensured. *)
     val ( + ) : ('n, 'm) matrix -> ('n, 'm) matrix -> ('n, 'm) matrix
 
-    (* Matrices multiplication. The dimensions compatibility is statically
-       ensured. *)
+    (** Matrices multiplication. The dimensions compatibility is statically
+        ensured. *)
     val ( * ) : ('n, 'm) matrix -> ('m, 'p) matrix -> ('n, 'p) matrix
 
-    (* Matrix exponentiation. The call [power m] returns a memoized function.
-       When one needs to compute several exponentiations of the same matrix, one
-       should perform the call [power m] once and used the returned function
-       each times one needs it. *)
+    (** Matrix exponentiation. The call [power m] returns a memoized function.
+        When one needs to compute several exponentiations of the same matrix,
+        one should perform the call [power m] once and used the returned
+        function each times one needs it. *)
     val power : ('n, 'n) matrix -> (int -> ('n, 'n) matrix)
 
   end

src/kernel_services/analysis/filter/nat.mli → src/libraries/arithmetic/nat.mli

@@ -20,10 +20,11 @@
 (*                                                                        *)
 (**************************************************************************)
 
-(* Encoding of the Peano arithmetic in OCaml type system. A value of type
-   [n nat] contains [n] at the value and the type level, allowing to express
-   properties on objects sizes and accesses for example. It is used by the
-   module Linear to represent vectors and matrices dimensions. *)
+(** Encoding of the Peano arithmetic in OCaml type system. A value of type
+    [n nat] contains [n] at the value and the type level, allowing to express
+    properties on objects sizes and accesses for example. It is used by the
+    module {!Linear} to represent vectors and matrices dimensions. *)
+
 type zero = |
 type 'n succ = |
 type 'n nat
@@ -35,14 +36,14 @@ val one  : zero succ nat
 val succ : 'n nat -> 'n succ nat
 val prev : 'n succ nat -> 'n nat
 
-(* The call [to_int n] returns an integer equal to n. This function complexity
-   is O(1). *)
+(** The call [to_int n] returns an integer equal to n. This function complexity
+    is O(1). *)
 val to_int : 'n nat -> int
 
-(* Returns a positive or null natural. If the given parameter is stricly
-   negative then 0 is returned. This function complexity is O(n). *)
+(** Returns a positive or null natural. If the given parameter is stricly
+    negative then [None] is returned. This function complexity is O(1). *)
 val of_int : int -> positive_or_null option
 
-(* Returns a strictly positive natural. If the given parameter is less or equal
-   than zero, then 1 is returned. This function complexity is O(n). *)
+(** Returns a strictly positive natural. If the given parameter is less or equal
+    than zero, then [None] is returned. This function complexity is O(1). *)
 val of_strictly_positive_int : int -> strictly_positive option

src/libraries/arithmetic/rational.ml

+(**************************************************************************)
+(*                                                                        *)
+(*  This file is part of Frama-C.                                         *)
+(*                                                                        *)
+(*  Copyright (C) 2007-2025                                               *)
+(*    CEA (Commissariat à l'énergie atomique et aux énergies              *)
+(*         alternatives)                                                  *)
+(*                                                                        *)
+(*  you can redistribute it and/or modify it under the terms of the GNU   *)
+(*  Lesser General Public License as published by the Free Software       *)
+(*  Foundation, version 2.1.                                              *)
+(*                                                                        *)
+(*  It is distributed in the hope that it will be useful,                 *)
+(*  but WITHOUT ANY WARRANTY; without even the implied warranty of        *)
+(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *)
+(*  GNU Lesser General Public License for more details.                   *)
+(*                                                                        *)
+(*  See the GNU Lesser General Public License version 2.1                 *)
+(*  for more details (enclosed in the file licenses/LGPLv2.1).            *)
+(*                                                                        *)
+(**************************************************************************)
+
+include Datatype.Rational
+type scalar = t
+
+let zero = Q.zero
+let one = Q.one
+let two = Q.of_int 2
+let ten = Z.of_int 10
+let pos_inf = Q.inf
+let neg_inf = Q.minus_inf
+let of_float = Q.of_float
+let to_float = Q.to_float
+let of_int = Q.of_int
+
+let pow10 e = Z.(pow ten e) |> Q.of_bigint
+
+let split_significant_exponent s =
+  match String.split_on_char 'e' s with
+  | [] | _ :: _ :: _ :: _ -> assert false
+  | [ s ] -> s, 0
+  | [ m ; e ] -> m, int_of_string e
+
+let significant_of_string significant =
+  match String.split_on_char '.' significant with
+  | [] | _ :: _ :: _ :: _ -> assert false
+  | [ m ] -> Q.of_string m
+  | [ integer ; fractional ] ->
+    let shift = pow10 (String.length fractional) in
+    Q.(of_string (integer ^ fractional) / shift)
+
+let of_string str =
+  let str = String.lowercase_ascii str in
+  let length = String.length str in
+  let last = String.get str (length - 1) in
+  let str = if last == 'f' then String.sub str 0 (length - 1) else str in
+  let significant, e = split_significant_exponent str in
+  let significant = significant_of_string significant in
+  let shift = if e >= 0 then pow10 e else Q.inv (pow10 ~-e) in
+  Q.(significant * shift)
+
+let pow2 e =
+  Q.(mul_2exp one e)
+
+(* r = ⌊log₂ (n / d)⌋ = ⌊log₂ n - log₂ d⌋ which cannot be calculated directly.
+   However, we have the three following properties:
+   - ⌊x⌋ + ⌊y⌋ ≤ ⌊x + y⌋ ≤ ⌊x⌋ + ⌊y⌋ + 1 ;
+   - ⌊-x⌋ = -⌈x⌉ ;
+   - ⌈x⌉ = ⌊x⌋ + 1 ;
+     Thus, we deduce that ⌊log₂ n⌋ - ⌊log₂ d⌋ - 1 ≤ r ≤ ⌊log₂ n⌋ - ⌊log₂ d⌋. *)
+let log2 q =
+  if Q.(q <= zero) then raise (Invalid_argument (Q.to_string q)) ;
+  let num = Q.num q |> Z.log2 and den = Q.den q |> Z.log2 in
+  Field.{ lower = num - den - 1 ; upper = num - den }
+
+let neg = Q.neg
+let abs = Q.abs
+let min = Q.min
+let max = Q.max
+
+let sqrt q =
+  if Q.sign q = 1 then
+    let num = Q.num q and den = Q.den q in
+    let acceptable_delta = Q.inv (Q.of_bigint @@ Z.pow ten 32) in
+    let rec approx_starting_at_scaling scaling =
+      let lower = Z.(sqrt (num * den * scaling * scaling)) in
+      let upper = Z.(lower + one) in
+      let denominator = Z.(den * scaling) in
+      let lower = Q.(make lower denominator) in
+      let upper = Q.(make upper denominator) in
+      let delta = Q.(upper - lower) in
+      if Q.(delta <= acceptable_delta) then Field.{ lower ; upper }
+      else approx_starting_at_scaling Z.(scaling * scaling)
+    in approx_starting_at_scaling ten
+  else { lower = neg_inf ; upper = pos_inf }
+
+let ( + ) = Q.( + )
+let ( - ) = Q.( - )
+let ( * ) = Q.( * )
+let ( / ) = Q.( / )
+
+let ( =  ) l r = Q.equal l r
+let ( != ) l r = not (l = r)
+let ( <= ) = Q.( <= )
+let ( <  ) = Q.( <  )
+let ( >= ) = Q.( >= )
+let ( >  ) = Q.( >  )

src/kernel_services/analysis/filter/field.ml → src/libraries/arithmetic/rational.mli

@@ -20,33 +20,7 @@
 (*                                                                        *)
 (**************************************************************************)
 
-module type S = sig
+(** Implementation of the {!Field} signature based on rational numbers.
+    @since Frama-C+dev *)
 
-  type scalar
-  include Datatype.S_with_collections with type t = scalar
-
-  val zero     : scalar
-  val one      : scalar
-  val infinity : scalar
-
-  val of_int : int -> scalar
-  val of_float : float -> scalar
-  val to_float : scalar -> float
-
-  val neg : scalar -> scalar
-  val abs : scalar -> scalar
-  val max : scalar -> scalar -> scalar
-  val min : scalar -> scalar -> scalar
-
-  val ( =  ) : scalar -> scalar -> bool
-  val ( <= ) : scalar -> scalar -> bool
-  val ( <  ) : scalar -> scalar -> bool
-  val ( >= ) : scalar -> scalar -> bool
-  val ( >  ) : scalar -> scalar -> bool
-
-  val ( + ) : scalar -> scalar -> scalar
-  val ( - ) : scalar -> scalar -> scalar
-  val ( * ) : scalar -> scalar -> scalar
-  val ( / ) : scalar -> scalar -> scalar
-
-end
+include Field.S with type t = Q.t

src/libraries/datatype/datatype.ml

@@ -1601,6 +1601,23 @@ module Filepath = struct
 end
 let filepath = Filepath.ty
 
+module Rational =
+  Make_with_collections
+    (struct
+      type t = Q.t
+      let name = "Datatype.Rational"
+      let reprs = [ Q.zero ; Q.one ]
+      let structural_descr = Structural_descr.t_abstract
+      let equal = Q.equal
+      let compare = Q.compare
+      let copy = identity
+      let rehash = identity
+      let mem_project = never_any_project
+      let pretty fmt q = Format.pp_print_float fmt (Q.to_float q)
+      let hash q = Stdlib.Hashtbl.hash (Z.hash (Q.num q), Z.hash (Q.den q))
+    end)
+let rational = Rational.ty
+
 (* ****************************************************************************)
 (** {3 Triple} *)
 (* ****************************************************************************)

src/libraries/datatype/datatype.mli

@@ -377,6 +377,9 @@ val formatter: Format.formatter Type.t
 module Integer: S_with_collections with type t = Integer.t
 val integer: Integer.t Type.t
 
+module Rational: S_with_collections with type t = Q.t
+val rational : Rational.t Type.t
+
 (** Type-safe strings representing normalized filepaths.
     See module {!Filepath.Normalized}.
     @since 18.0-Argon *)

tests/float/oracle/linear_filter_test.res.oracle

 [kernel] Parsing linear_filter_test.i (no preprocessing)
 [kernel] Circle :
-  * 0 : [-31.3829787 .. 31.3829787]
-  * 1 : [-31.3829787 .. 31.3829787]
+  * 0 : [-31.3829787234 .. 31.3829787234]
+  * 1 : [-31.3829787234 .. 31.3829787234]
 [kernel] Simple :
-  * 0 : [5.7267253E-1 .. 2.4273275]
-  * 1 : [1.5726725 .. 3.4273275]
+  * 0 : [0.572672526718 .. 2.42732747328]
+  * 1 : [1.57267252672 .. 3.42732747328]