why3execute: Add execution for floats and reals (depends from mpfr)

This adds an optional dependency on MPFR to run why3execute for floats. It
is also used for reals (represented as intervals of floats).
This commit does the following changes:
- update the configure/Makefile to allow MPFR dependency
- adds a MPFR wrapper so that why3execute can still be compiled if MPFR is
  not installed. In this case, an exception is raised when executing on
  real/float.
- some updates are made to the standard library so that real/float
  functions are part of the program world (and can be executed).
- pinterp changes to make elementary functions from real and float
  executable.
- add some tests under bench/interp for manual testing of this feature

Note that the correctness of the results given for reals depends on the
precision. A too low precision may give unexpected results.

.merlin.in

@@ -41,4 +41,4 @@ B plugins/tptp
 B plugins/python
 B lib/why3
 
-PKG str unix num dynlink @ZIPLIB@ @LABLGTKPKG@ @META_OCAMLGRAPH@ @JSOFOCAMLPKG@
+PKG str unix num dynlink @ZIPLIB@ @LABLGTKPKG@ @META_OCAMLGRAPH@ @JSOFOCAMLPKG@ @MLMPFR@

Makefile.in

@@ -87,7 +87,7 @@ EMACS = @EMACS@
 #PSVIEWER  = @PSVIEWER@
 #PDFVIEWER = @PDFVIEWER@
 
-INCLUDES = @WHY3INCLUDE@ @ZIPINCLUDE@ @MENHIRINCLUDE@ @NUMINCLUDE@
+INCLUDES = @WHY3INCLUDE@ @ZIPINCLUDE@ @MENHIRINCLUDE@ @NUMINCLUDE@ @MLMPFRINCLUDE@
 
 # warnings are enabled and non fatal by default, except:
 # - disabled:
@@ -113,8 +113,8 @@ WARNINGS = A-4-9-41-44-45-50-52@5@8@48
 OFLAGS = -w $(WARNINGS) -safe-string -keep-locs -bin-annot -dtypes -g $(INCLUDES)
 BFLAGS = -w $(WARNINGS) -safe-string -keep-locs -bin-annot -dtypes -g $(INCLUDES)
 
-OLINKFLAGS = -linkall $(EXTCMXA)
-BLINKFLAGS = -linkall $(EXTCMA)
+OLINKFLAGS = -linkall @MLMPFR_LINK@ $(EXTCMXA)
+BLINKFLAGS = -linkall @MLMPFR_LINK@ $(EXTCMA)
 
 ifeq (@enable_profiling@,yes)
 OFLAGS += -g -p
@@ -126,7 +126,7 @@ CMIHACK = -intf-suffix .cmi
 # external libraries common to all binaries
 
 EXTOBJS = menhirLib
-EXTLIBS = str unix nums dynlink @ZIPLIB@ @WHY3LIB@
+EXTLIBS = str unix nums dynlink @ZIPLIB@ @MLMPFR@ @WHY3LIB@
 
 EXTCMA	= $(addsuffix .cmo,$(EXTOBJS)) $(addsuffix .cma,$(EXTLIBS))
 EXTCMXA = $(addsuffix .cmx,$(EXTOBJS)) $(addsuffix .cmxa,$(EXTLIBS))
@@ -187,7 +187,7 @@ GENERATED =
 LIBGENERATED = src/util/config.ml \
 	       src/util/rc.ml src/util/lexlib.ml \
 	       src/util/json_parser.mli src/util/json_parser.ml \
-	       src/util/json_lexer.ml \
+	       src/util/json_lexer.ml src/util/mlmpfr_wrapper.ml \
 	       src/parser/lexer.ml \
 	       src/parser/parser.mli src/parser/parser.ml \
 	       src/driver/driver_parser.mli src/driver/driver_parser.ml \
@@ -197,7 +197,7 @@ LIBGENERATED = src/util/config.ml \
 	       src/session/compress.ml src/session/xml.ml \
 	       src/session/strategy_parser.ml
 
-LIB_UTIL = config bigInt util opt lists strings \
+LIB_UTIL = config bigInt mlmpfr_wrapper util opt lists strings \
 	   pp extmap extset exthtbl weakhtbl \
 	   hashcons wstdlib exn_printer \
 	   json_base json_parser json_lexer \
@@ -213,7 +213,7 @@ LIB_DRIVER = prove_client call_provers driver_ast driver_parser driver_lexer dri
 		collect_data_model parse_smtv2_model_lexer parse_smtv2_model \
 		parse_smtv2_model
 
-LIB_MLW = ity expr pdecl eval_match typeinv vc pmodule dexpr \
+LIB_MLW = ity expr pdecl eval_match typeinv vc pmodule dexpr big_real \
           pinterp mltree compile mlinterp pdriver cprinter ml_printer \
 	  ocaml_printer cakeml_printer
 
@@ -269,6 +269,16 @@ $(LIBCMX): OFLAGS += -for-pack Why3
 
 $(LIBDEP): $(LIBGENERATED)
 
+# Mlmpfr
+
+ifeq (@found_mlmpfr@,yes)
+src/util/mlmpfr_wrapper.ml: config.status src/util/mlmpfr_real.ml
+	cp src/util/mlmpfr_real.ml $@
+else
+src/util/mlmpfr_wrapper.ml: config.status src/util/mlmpfr_dummy.ml
+	cp src/util/mlmpfr_dummy.ml $@
+endif
+
 # Ocamlzip
 
 ifeq (@enable_zip@,yes)
@@ -2016,6 +2026,7 @@ MODULESTODOC = \
 	util/extmap util/extset util/exthtbl \
 	util/weakhtbl util/wstdlib util/rc util/debug \
 	util/loc util/pp util/bigInt util/number \
+	util/mlmpfr_wrapper \
 	core/ident core/ty core/term core/decl core/coercion core/theory \
 	core/env core/task core/trans core/pretty core/printer \
 	parser/typing \

bench/interp/float32.mlw

+module N
+
+  use ieee_float.Float32
+
+  let f (x: t)
+  =
+    add RNE x (1.0: t)
+
+  let g () =
+    add RNE (24.0: t) (1.0: t)
+
+  exception BenchFailure
+
+  let bench1 () raises { BenchFailure -> false } =
+    let x = f (1.0: t) in
+    if not (eq x (2.0: t)) then raise BenchFailure;
+    x
+
+  let bench2 () raises { BenchFailure -> false } =
+    let o = (0x1p-23: t) in
+    let x = sqrt RNE ((1.0: t) .+ o .- (1.0: t) .- o) in (* Expected result: 0.0 *)
+    if not (eq x (0.0: t)) then raise BenchFailure;
+    let o = (0x1p-24: t) in
+    let x = sqrt RNE ((1.0: t) .+ o .- (1.0: t) .- o) in (* Expected result: Nan *)
+    if eq x x then raise BenchFailure;
+    x
+
+  let bench3 ()
+    raises { BenchFailure -> false }
+    ensures { eq result (0.0: t)} =
+    let o = (0x1p-129: t) in
+    let v = o .* o in
+    if not (eq v (0.0:t)) then raise BenchFailure;
+    v
+
+  let bench4 ()
+    raises { BenchFailure -> false }
+    ensures { is_infinite result } =
+    let o = (0x1p+64: t) in
+    let res = o .* o in
+    if not (is_infinite res) then raise BenchFailure;
+    res
+
+  let bench5 ()
+    raises { BenchFailure -> false }
+    ensures { eq result (1.0:t) } =
+    let o = add RNE (1.0:t) (0x1p-24:t) in
+    if not (eq o (1.0:t)) then raise BenchFailure;
+    o
+
+  let bench6 ()
+    raises { BenchFailure -> false }
+    ensures { not (eq result (1.0:t)) } =
+    let o = add RNE (1.0:t) (0x1p-23:t) in
+    if eq o (1.0:t) then raise BenchFailure;
+    o
+
+  let bench7 ()
+    raises { BenchFailure -> false }
+    ensures { eq result (0.0: t)} =
+    let o = (0x1p-149: t) in
+    let o2 = o ./ (2.0: t) in
+    if not (eq o2 (0.0:t)) then raise BenchFailure;
+    o2
+
+  let bench8 ()
+    raises { BenchFailure -> false }
+    ensures { not (eq result (0.0: t))} =
+    let o = (0x1p-148: t) in
+    let o2 = o ./ (2.0: t) in
+    if (eq o2 (0.0:t)) then raise BenchFailure;
+    o2
+
+end

bench/interp/float64.mlw

+module N
+
+  use ieee_float.Float64
+
+  let f (x: t)
+  =
+    add RNE x (1.0: t)
+
+  let g () =
+    add RNE (24.0: t) (1.0: t)
+
+  exception BenchFailure
+
+  let bench1 () raises { BenchFailure -> false } =
+    let x = f (1.0: t) in
+    if not (eq x (2.0: t)) then raise BenchFailure;
+    x
+
+  let bench2 () raises { BenchFailure -> false } =
+    let o = (0x1p-52: t) in
+    let x = sqrt RNE ((1.0: t) .+ o .- (1.0: t) .- o) in (* Expected result: 0.0 *)
+    if not (eq x (0.0: t)) then raise BenchFailure;
+    let o = (0x1p-53: t) in
+    let x = sqrt RNE ((1.0: t) .+ o .- (1.0: t) .- o) in (* Expected result: Nan *)
+    if eq x x then raise BenchFailure;
+    x
+
+  let bench3 ()
+    raises { BenchFailure -> false }
+    ensures { eq result (0.0: t)} =
+    let o = (0x1p-1023: t) in
+    let v = o .* o in
+    if not (eq v (0.0:t)) then raise BenchFailure;
+    v
+
+  let bench4 ()
+    raises { BenchFailure -> false }
+    ensures { is_infinite result } =
+    let o = (0x1p+1023: t) in
+    let res = o .* o in
+    if not (is_infinite res) then raise BenchFailure;
+    res
+
+  let bench5 ()
+    raises { BenchFailure -> false }
+    ensures { eq result (1.0:t) } =
+    let o = add RNE (1.0:t) (0x1p-53:t) in
+    if not (eq o (1.0:t)) then raise BenchFailure;
+    o
+
+  let bench6 ()
+    raises { BenchFailure -> false }
+    ensures { not (eq result (1.0:t)) } =
+    let o = add RNE (1.0:t) (0x1p-52:t) in
+    if eq o (1.0:t) then raise BenchFailure;
+    o
+
+  let bench7 ()
+    raises { BenchFailure -> false }
+    ensures { eq result (0.0: t)} =
+    let o = (0x1p-1074: t) in
+    let o2 = o ./ (2.0: t) in
+    if not (eq o2 (0.0:t)) then raise BenchFailure;
+    o2
+
+  let bench8 ()
+    raises { BenchFailure -> false }
+    ensures { not (eq result (0.0: t))} =
+    let o = (0x1p-1073: t) in
+    let o2 = o ./ (2.0: t) in
+    if (eq o2 (0.0:t)) then raise BenchFailure;
+    o2
+
+end

bench/interp/real.mlw

+module R
+
+  use real.Real
+  use ref.Ref
+  use real.Square
+
+  exception BenchFailure
+
+  let bench1 ()
+    (* Tries to calculate sqrt(2) *)
+    diverges
+    raises { BenchFailure -> false }
+    ensures { result = 10.0 } =
+    let incr = 0.002 (*(1.0 / 10000.0)*) in
+    let i : ref real = ref 0.0 in
+    while (!i * !i <= 2.0) do
+      i := !i + incr;
+    done;
+    !i
+
+  let bench2 ()
+(*    raises { BenchFailure -> false }*)
+    ensures { result = 4.0 } =
+    sqrt (16.0)
+
+  use real.Trigonometry
+
+  let bench3 () =
+    sqrt (pi)
+
+end

configure.in

@@ -688,6 +688,45 @@ elif test "$enable_js_of_ocaml" != yes; then
    reason_web_ide=" (Javascript support not available)"
 fi
 
+# Mlmpfr
+MLMPFR_LINK=
+if test "$enable_js_of_ocaml" != no; then
+   found_mlmpfr=no
+   reason_mlmpfr=" (Cannot allow mlmpfr with js_of_ocaml) "
+else
+   if test "$USEOCAMLFIND" = yes; then
+      DIR=$(ocamlfind query mlmpfr)
+      if test -n "$DIR"; then
+         echo "ocamlfind found mlmpfr in $DIR"
+         # Test that mpfr version is higher than 4.0.0 (because of
+         # Faithful constructor incompatibility).
+         MPFR_VER=$(grep -q "4.0.0" $DIR/META 2> /dev/null)
+         EXIT_CODE=$?
+         if test $EXIT_CODE = 0; then
+            MLMPFRINCLUDE="-I $DIR"
+            MLMPFR=mlmpfr
+            found_mlmpfr=yes
+            MLMPFR_LINK="-cclib -lmpfr"
+            AC_CHECK_FILE($DIR/mlmpfr.cma,,found_mlmpfr=yes)
+         else
+            reason_mlmpfr=" (Version of mlmpfr is not 4.0.0) "
+            found_mlmpfr=no
+         fi
+      else
+         reason_mlmpfr=" (Could not find mlmpfr) "
+         found_mlmpfr=no
+      fi
+   else
+      reason_mlmpfr=" (Should use ocamlfind to use mlmpfr) "
+      found_mlmpfr=no
+   fi
+fi
+if test "$found_mlmpfr" = no; then
+   MLMPFRINCLUDE=
+   MLMPFR=
+fi
+
+
 # Coq
 
 enable_coq_support=yes
@@ -937,6 +976,10 @@ AC_SUBST(COMPILERLIBS)
 AC_SUBST(COMPILERLIBSPKG)
 
 AC_SUBST(NUMINCLUDE)
+AC_SUBST(MLMPFRINCLUDE)
+AC_SUBST(MLMPFR)
+AC_SUBST(MLMPFR_LINK)
+AC_SUBST(found_mlmpfr)
 
 AC_SUBST(enable_zarith)
 AC_SUBST(BIGINTINCLUDE)
@@ -1032,6 +1075,7 @@ echo "    Native compilation      : $enable_native_code"
 echo "    Profiling               : $enable_profiling"
 echo "    Compiler plugins        : $enable_compiler_plugins$reason_compiler_plugins"
 echo "    Javascript support      : $enable_js_of_ocaml$reason_js_of_ocaml"
+echo "    Mpfr support            : $found_mlmpfr$reason_mlmpfr"
 echo "Components"
 echo "    Why3 library            : $enable_why3_lib"
 echo "    GTK IDE                 : $enable_ide$reason_ide"

src/mlw/big_real.ml

+(* Experimental module under development *)
+
+exception Undetermined
+
+open Mlmpfr_wrapper
+type real = mpfr_float * mpfr_float
+
+(* TODO in all the following getting to plus_infinity or to zero on both side of
+   the interval should always be seen as an error. The precision should be
+   increased in those cases.
+   For example, we cannot simplify multiplication of [0; 0] with ]+inf; +inf[
+   because both may mean that the underlying floats overflow/underflow and the
+   result can be far from the solutions we get. Example: divide a number by 10 a
+   lot and then multiply the result by 10 the same number of time + 1. You would
+   obtain 0 whereas you would like to obtain number * 10.
+*)
+
+(* TODO precision cannot be changed once launched the first time. So we need
+   to init it once. *)
+let init, set_exponents, get_prec =
+  let emin = ref (-148) in
+  let emax = ref 128 in
+  let prec = ref 24 in
+  (fun emin_i emax_i prec_i ->
+    emin := emin_i;
+    emax := emax_i;
+    prec := prec_i),
+  (fun () ->
+    set_emin !emin;
+    set_emax !emax),
+  (fun () -> !prec)
+
+let add (xmin, xmax) (ymin, ymax) =
+  (* Exponents can be changed if floats occur in the code. *)
+  set_exponents ();
+  let prec = get_prec () in
+  let res_min = add ~prec ~rnd:Toward_Minus_Infinity xmin ymin in
+  let res_max = add ~prec ~rnd:Toward_Plus_Infinity xmax ymax in
+  (res_min, res_max)
+
+let neg (xmin, xmax) =
+  set_exponents ();
+  let prec = get_prec () in
+  let res_min = neg ~prec ~rnd:Toward_Minus_Infinity xmax in
+  let res_max = neg ~prec ~rnd:Toward_Plus_Infinity xmin in
+  (res_min, res_max)
+
+let mul (xmin, xmax) (ymin, ymax) =
+  set_exponents ();
+  let prec = get_prec () in
+  let min = min ~prec ~rnd:Toward_Minus_Infinity in
+  let max = max ~prec ~rnd:Toward_Plus_Infinity in
+  let mul1_min = mul ~prec ~rnd:Toward_Minus_Infinity xmin ymin in
+  let mul2_min = mul ~prec ~rnd:Toward_Minus_Infinity xmin ymax in
+  let mul3_min = mul ~prec ~rnd:Toward_Minus_Infinity xmax ymin in
+  let mul4_min = mul ~prec ~rnd:Toward_Minus_Infinity xmax ymax in
+  let res_min = List.fold_left min mul1_min [mul2_min; mul3_min; mul4_min] in
+  let mul1_max = mul ~prec ~rnd:Toward_Plus_Infinity xmin ymin in
+  let mul2_max = mul ~prec ~rnd:Toward_Plus_Infinity xmin ymax in
+  let mul3_max = mul ~prec ~rnd:Toward_Plus_Infinity xmax ymin in
+  let mul4_max = mul ~prec ~rnd:Toward_Plus_Infinity xmax ymax in
+  let res_max = List.fold_left max mul1_max [mul2_max; mul3_max; mul4_max] in
+  (res_min, res_max)
+
+let inv (xmin, xmax) =
+  set_exponents ();
+  let prec = get_prec () in
+  (* If 0 is inside the interval we cannot compute the expression. The precision
+     should be increased. *)
+  if signbit xmin <> signbit xmax then
+    raise Undetermined
+  else
+    let one = make_from_int ~prec ~rnd:Toward_Minus_Infinity 1 in
+    let inv1_min = div ~prec ~rnd:Toward_Minus_Infinity one xmin in
+    let inv2_min = div ~prec ~rnd:Toward_Minus_Infinity one xmax in
+    let res_min = min inv1_min inv2_min in
+    let inv1_max = div ~prec ~rnd:Toward_Plus_Infinity one xmin in
+    let inv2_max = div ~prec ~rnd:Toward_Plus_Infinity one xmax in
+    let res_max = max inv1_max inv2_max in
+    (res_min, res_max)
+
+let div x y =
+  mul x (inv y)
+
+let sqrt (xmin, xmax) =
+  set_exponents ();
+  let prec = get_prec() in
+  if lessequal_p (make_zero ~prec Positive) xmin then
+    let res_min = sqrt ~rnd:Toward_Minus_Infinity ~prec xmin in
+    let res_max = sqrt ~rnd:Toward_Plus_Infinity ~prec xmax in
+    (res_min, res_max)
+  else
+    raise Undetermined
+
+let real_from_str s =
+  let prec = get_prec () in
+  let n = make_from_str ~prec ~base:10 s in
+  (n, n)
+
+let real_from_fraction p q =
+  let p = real_from_str p in
+  let q = real_from_str q in
+  div p q
+
+let print_real fmt r =
+  let (x, y) = r in
+  let x = get_formatted_str ~base:10 x in
+  let y = get_formatted_str ~base:10 y in
+  Format.fprintf fmt "[%s, %s]" x y
+
+let eq (xmin, xmax) (ymin, ymax) =
+  if (equal_p xmin xmax) && (equal_p ymin ymax) then
+    equal_p xmin ymin
+  else
+    raise Undetermined
+
+let lt x y =
+  if less_p (snd x) (fst y) then
+    true
+  else
+    if less_p (fst x) (snd y) then
+      false
+    else
+      raise Undetermined
+
+let le x y =
+  if lessequal_p (snd x) (fst y) then
+    true
+  else
+    if less_p (snd y) (fst x) then
+      false
+    else
+      raise Undetermined
+
+let pi =
+  let prec = get_prec () in
+  const_pi prec, const_pi prec

src/mlw/big_real.mli

+type real
+
+(* Exception raised when intervals do not allow decision of a comparison
+   function. For example, equality on non-singleton intervals. *)
+exception Undetermined
+
+val init: int -> int -> int -> unit
+
+val print_real: Format.formatter -> real -> unit
+
+val real_from_str: string -> real
+val real_from_fraction: string -> string -> real
+
+(* Real operations *)
+val neg: real -> real
+val add: real -> real -> real
+val mul: real -> real -> real
+val div: real -> real -> real
+val sqrt: real -> real
+
+(* Real comparisons *)
+val eq: real -> real -> bool
+val lt: real -> real -> bool
+val le: real -> real -> bool
+
+(* Constants *)
+val pi: real

src/mlw/pinterp.ml

@@ -20,13 +20,20 @@ let debug = Debug.register_info_flag "trace_exec"
 
 open Ity
 open Expr
-
+open Big_real
+open Mlmpfr_wrapper
 
 (* module Nummap = Map.Make(BigInt) *)
 
+type float_mode = mpfr_rnd_t
+type big_float = mpfr_float
+
 type value =
   | Vconstr of rsymbol * field list
   | Vnum of BigInt.t
+  | Vreal of real
+  | Vfloat_mode of float_mode
+  | Vfloat of big_float
   | Vbool of bool
   | Vvoid
   | Varray of value array
@@ -40,6 +47,15 @@ let field_get f =
 let constr rs vl =
   Vconstr(rs,List.map (fun v -> Fimmutable v) vl)
 
+let mode_to_string m =
+  match m with
+  | To_Nearest            -> "RNE"
+  | Away_From_Zero        -> "RNA"
+  | Toward_Plus_Infinity  -> "RTP"
+  | Toward_Minus_Infinity -> "RTN"
+  | Toward_Zero           -> "RTZ"
+  | Faithful              -> assert false
+
 let rec print_value fmt v =
   match v with
   | Vnum n ->
@@ -49,6 +65,16 @@ let rec print_value fmt v =
       fprintf fmt "(%s)" (BigInt.to_string n)
   | Vbool b ->
     fprintf fmt "%b" b
+  | Vreal r ->
+    print_real fmt r
+  | Vfloat f ->
+      (* Getting "@" is intentional in mlmpfr library for bases higher than 10.
+         So, we keep this notation. *)
+      let hexadecimal = get_formatted_str ~base:16 f in
+      let decimal = get_formatted_str ~base:10 f in
+      fprintf fmt "%s (%s)" decimal hexadecimal
+  | Vfloat_mode m ->
+      fprintf fmt "%s" (mode_to_string m)
   | Vvoid ->
     fprintf fmt "()"
   | Varray a ->
@@ -134,10 +160,8 @@ let rec matching env (v:value) p =
 
 exception NotNum
 
-let big_int_of_const c =
-  match c with
-    | Number.ConstInt i -> i.Number.il_int
-    | _ -> raise NotNum
+let big_int_of_const i =
+  i.Number.il_int
 
 let big_int_of_value v =
   match v with
@@ -178,6 +202,84 @@ let eval_int_rel op ls l =
     end
   | _ -> constr ls l
 
+(* This initialize Mpfr for float32 behavior *)
+let initialize_float32 () =
+  set_default_prec 24;
+  set_emin (-148);
+  set_emax 128
+
+(* This initialize Mpfr for float64 behavior *)
+let initialize_float64 () =
+  set_default_prec 53;
+  set_emin (-1073);
+  set_emax 1024
+
+type 'a float_arity =
+  | Mode1: (float_mode -> big_float -> big_float) float_arity (* Unary op *)
+  | Mode2: (float_mode -> big_float -> big_float -> big_float) float_arity (* binary op *)
+  | Mode3: (float_mode -> big_float -> big_float -> big_float -> big_float) float_arity (* ternary op *)
+  | Mode_rel: (big_float -> big_float -> bool) float_arity (* binary predicates *)
+  | Mode_rel1: (big_float -> bool) float_arity (* Unary predicate *)
+
+let use_fmode (fmode: int) =
+  match fmode with
+  | 32 -> initialize_float32 ()
+  | 64 -> initialize_float64 ()
+  | _ -> raise CannotCompute
+
+let eval_float:
+  type a. int -> a float_arity -> a -> Expr.rsymbol -> value list -> value =
+  (fun fmode ty op ls l ->
+     (* Simulate IEEE behavior *)
+     let subnormalize = subnormalize ~rnd:To_Nearest in
+     (* Set the exponent depending on Float type that are used: 32 or 64 *)
+     use_fmode fmode;
+     try
+       begin match ty, l with
+         | Mode1, [Vfloat_mode mode; Vfloat f] ->
+             Vfloat (subnormalize (op mode f))
+         | Mode2, [Vfloat_mode mode; Vfloat f1; Vfloat f2] ->
+             Vfloat (subnormalize (op mode f1 f2))
+         | Mode3, [Vfloat_mode mode; Vfloat f1; Vfloat f2; Vfloat f3] ->
+             Vfloat (subnormalize (op mode f1 f2 f3))
+         | Mode_rel, [Vfloat f1; Vfloat f2] ->
+             Vbool (op f1 f2)
+         | Mode_rel1, [Vfloat f] ->
+             Vbool (op f)
+         | _ -> constr ls l
+       end
+     with
+     | Mlmpfr_wrapper.Not_Implemented -> raise CannotCompute
+     | _ -> assert false
+  )
+
+type 'a real_arity =
+  | Modeconst: real real_arity
+  | Mode1r: (real -> real) real_arity
+  | Mode2r: (real -> real -> real) real_arity
+  | Mode_relr: (real -> real -> bool) real_arity
+
+let eval_real : type a. a real_arity -> a -> Expr.rsymbol -> value list -> value =
+  (fun ty op ls l ->
+     try
+       begin match ty, l with
+         | Mode1r, [Vreal r] ->
+             Vreal (op r)
+         | Mode2r, [Vreal r1; Vreal r2] ->
+             Vreal (op r1 r2)
+         | Mode_relr, [Vreal r1; Vreal r2] ->
+             Vbool (op r1 r2)