From ae5ddf91511de1815810d1255efac39ce8988c25 Mon Sep 17 00:00:00 2001 From: Sylvain Dailler Date: Mon, 11 Mar 2019 17:19:35 +0100 Subject: [PATCH] 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 | 2 +- Makefile.in | 25 ++- bench/interp/float32.mlw | 74 +++++++++ bench/interp/float64.mlw | 74 +++++++++ bench/interp/real.mlw | 31 ++++ configure.in | 44 +++++ src/mlw/big_real.ml | 137 ++++++++++++++++ src/mlw/big_real.mli | 27 ++++ src/mlw/pinterp.ml | 314 ++++++++++++++++++++++++++++++------ src/mlw/pinterp.mli | 7 +- src/tools/why3execute.ml | 4 +- src/util/mlmpfr_dummy.ml | 59 +++++++ src/util/mlmpfr_real.ml | 6 + src/util/mlmpfr_wrapper.mli | 69 ++++++++ stdlib/ieee_float.mlw | 10 +- stdlib/real.mlw | 4 +- 16 files changed, 823 insertions(+), 64 deletions(-) create mode 100644 bench/interp/float32.mlw create mode 100644 bench/interp/float64.mlw create mode 100644 bench/interp/real.mlw create mode 100644 src/mlw/big_real.ml create mode 100644 src/mlw/big_real.mli create mode 100644 src/util/mlmpfr_dummy.ml create mode 100644 src/util/mlmpfr_real.ml create mode 100644 src/util/mlmpfr_wrapper.mli diff --git a/.merlin.in b/.merlin.in index 86c957d1b..4e7d57b47 100644 --- a/.merlin.in +++ b/.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@ diff --git a/Makefile.in b/Makefile.in index b9c4f96ed..35f0d2e79 100644 --- a/Makefile.in +++ b/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 \ diff --git a/bench/interp/float32.mlw b/bench/interp/float32.mlw new file mode 100644 index 000000000..3b45fc8c5 --- /dev/null +++ b/bench/interp/float32.mlw @@ -0,0 +1,74 @@ +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 diff --git a/bench/interp/float64.mlw b/bench/interp/float64.mlw new file mode 100644 index 000000000..572322f80 --- /dev/null +++ b/bench/interp/float64.mlw @@ -0,0 +1,74 @@ +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 diff --git a/bench/interp/real.mlw b/bench/interp/real.mlw new file mode 100644 index 000000000..0fb43921d --- /dev/null +++ b/bench/interp/real.mlw @@ -0,0 +1,31 @@ +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 diff --git a/configure.in b/configure.in index 86af41c72..d115c42bb 100644 --- a/configure.in +++ b/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" diff --git a/src/mlw/big_real.ml b/src/mlw/big_real.ml new file mode 100644 index 000000000..3fdd05287 --- /dev/null +++ b/src/mlw/big_real.ml @@ -0,0 +1,137 @@ +(* 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 diff --git a/src/mlw/big_real.mli b/src/mlw/big_real.mli new file mode 100644 index 000000000..337c71842 --- /dev/null +++ b/src/mlw/big_real.mli @@ -0,0 +1,27 @@ +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 diff --git a/src/mlw/pinterp.ml b/src/mlw/pinterp.ml index a869d807f..f414617f0 100644 --- a/src/mlw/pinterp.ml +++ b/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) + | Modeconst, [] -> + Vreal op + | _ -> constr ls l + end + with + | Big_real.Undetermined -> (* Cannot decide interval comparison *) + constr ls l + | Mlmpfr_wrapper.Not_Implemented -> + raise CannotCompute + | _ -> assert false + ) let rec default_value_of_type env ity = match ity.ity_node with @@ -274,46 +376,126 @@ let exec_array_set _ args = end | _ -> assert false +let builtin_mode _kn _its = () + +let builtin_float_type _kn _its = () +(** Description of modules *) let built_in_modules = - [ + let bool_module = ["bool"],"Bool", [], [ "True", eval_true ; "False", eval_false ; - ] ; - ["int"],"Int", [], - [ Ident.op_infix "+", eval_int_op BigInt.add; - (* defined as x+(-y) - Ident.op_infix "-", eval_int_op BigInt.sub; *) - Ident.op_infix "*", eval_int_op BigInt.mul; - Ident.op_prefix "-", eval_int_uop BigInt.minus; - Ident.op_infix "=", eval_int_rel BigInt.eq; - Ident.op_infix "<", eval_int_rel BigInt.lt; - Ident.op_infix "<=", eval_int_rel BigInt.le; - Ident.op_infix ">", eval_int_rel BigInt.gt; - Ident.op_infix ">=", eval_int_rel BigInt.ge; - ] ; - ["int"],"MinMax", [], - [ "min", eval_int_op BigInt.min; - "max", eval_int_op BigInt.max; - ] ; - ["int"],"ComputerDivision", [], - [ "div", eval_int_op BigInt.computer_div; - "mod", eval_int_op BigInt.computer_mod; - ] ; - ["int"],"EuclideanDivision", [], - [ "div", eval_int_op BigInt.euclidean_div; - "mod", eval_int_op BigInt.euclidean_mod; - ] ; - ["array"],"Array", + ] + in + let int_modules = + [ + ["int"],"Int", [], + [ Ident.op_infix "+", eval_int_op BigInt.add; + (* defined as x+(-y) + Ident.op_infix "-", eval_int_op BigInt.sub; *) + Ident.op_infix "*", eval_int_op BigInt.mul; + Ident.op_prefix "-", eval_int_uop BigInt.minus; + Ident.op_infix "=", eval_int_rel BigInt.eq; + Ident.op_infix "<", eval_int_rel BigInt.lt; + Ident.op_infix "<=", eval_int_rel BigInt.le; + Ident.op_infix ">", eval_int_rel BigInt.gt; + Ident.op_infix ">=", eval_int_rel BigInt.ge; + ] ; + ["int"],"MinMax", [], + [ "min", eval_int_op BigInt.min; + "max", eval_int_op BigInt.max; + ] ; + ["int"],"ComputerDivision", [], + [ "div", eval_int_op BigInt.computer_div; + "mod", eval_int_op BigInt.computer_mod; + ] ; + ["int"],"EuclideanDivision", [], + [ "div", eval_int_op BigInt.euclidean_div; + "mod", eval_int_op BigInt.euclidean_mod; + ] + ] + in + let array_module = + ["array"],"Array", ["array", builtin_array_type], ["make", exec_array_make ; "empty", exec_array_empty ; Ident.op_get "", exec_array_get ; "length", exec_array_length ; Ident.op_set "", exec_array_set ; - ] ; - ] + ] + in + + let mode_module = + ["ieee_float"], "RoundingMode", + ["mode", builtin_mode], + ["RNE", (fun _ _ -> Vfloat_mode To_Nearest); + "RNA", (fun _ _ -> Vfloat_mode Away_From_Zero); + "RTP", (fun _ _ -> Vfloat_mode Toward_Plus_Infinity); + "RTN", (fun _ _ -> Vfloat_mode Toward_Minus_Infinity); + "RTZ", (fun _ _ -> Vfloat_mode Toward_Zero); + ] + in + let float_modules tyb ~prec m = + ["ieee_float"], m, + ["t", builtin_float_type], + ["zeroF", (fun _rs _l -> + Vfloat (make_zero ~prec Positive)); + "add", eval_float tyb Mode2 (fun rnd -> add ~rnd ~prec); + "sub", eval_float tyb Mode2 (fun rnd -> sub ~rnd ~prec); + "mul", eval_float tyb Mode2 (fun rnd -> mul ~rnd ~prec); + "div", eval_float tyb Mode2 (fun rnd -> div ~rnd ~prec); + "abs", eval_float tyb Mode1 (fun rnd -> abs ~rnd ~prec); + "neg", eval_float tyb Mode1 (fun rnd -> neg ~rnd ~prec); + "fma", eval_float tyb Mode3 (fun rnd -> fma ~rnd ~prec); + "sqrt", eval_float tyb Mode1 (fun rnd -> sqrt ~rnd ~prec); + "roundToIntegral", eval_float tyb Mode1 (fun rnd -> rint ~rnd ~prec); +(* Intentionnally removed from programs + "min", eval_float_minmax min; + "max", eval_float_minmax max; +*) + "le", eval_float tyb Mode_rel lessequal_p; + "lt", eval_float tyb Mode_rel less_p; + "eq", eval_float tyb Mode_rel equal_p; + "is_zero", eval_float tyb Mode_rel1 zero_p; + "is_infinite", eval_float tyb Mode_rel1 inf_p; + "is_nan", eval_float tyb Mode_rel1 nan_p; + "is_positive", eval_float tyb Mode_rel1 (fun s -> signbit s = Positive); + "is_negative", eval_float tyb Mode_rel1 (fun s -> signbit s = Negative); + ] + in + + let real_module = + ["real"], "Real", [], + [Ident.op_infix "=", eval_real Mode_relr Big_real.eq; + Ident.op_infix "<", eval_real Mode_relr Big_real.lt; + Ident.op_infix "<=", eval_real Mode_relr Big_real.le; + Ident.op_prefix "-", eval_real Mode1r Big_real.neg; + Ident.op_infix "+", eval_real Mode2r Big_real.add; + Ident.op_infix "*", eval_real Mode2r Big_real.mul; + Ident.op_infix "/", eval_real Mode2r Big_real.div; + ] + in + let real_square_module = + ["real"], "Square", [], + ["sqrt", eval_real Mode1r Big_real.sqrt; + ] + in + let real_trigo_module = + ["real"], "Trigonometry", [], + ["pi", eval_real Modeconst Big_real.pi; + ] + in + bool_module :: (int_modules @ [array_module; + real_module; + real_square_module; + real_trigo_module; + mode_module; + float_modules 32 ~prec:24 "Float32"; + float_modules 64 ~prec:53 "Float64";]) + +exception CannotFind of (Env.pathname * string * string) let add_builtin_mo env (l,n,t,d) = let mo = Pmodule.read_module env l n in @@ -321,19 +503,24 @@ let add_builtin_mo env (l,n,t,d) = let kn = mo.Pmodule.mod_known in List.iter (fun (id,r) -> - let its = Pmodule.ns_find_its exp [id] in - r kn its) + let its = + try Pmodule.ns_find_its exp [id] with + | Not_found -> raise (CannotFind (l, n, id)) + in + r kn its) t; List.iter (fun (id,f) -> - let ps = Pmodule.ns_find_rs exp [id] in - Hrs.add builtin_progs ps f) + let ps = + try Pmodule.ns_find_rs exp [id] with + | Not_found -> raise (CannotFind (l, n, id)) + in + Hrs.add builtin_progs ps f) d let get_builtin_progs lib = List.iter (add_builtin_mo lib) built_in_modules - let get_pvs env pvs = let t = try @@ -477,6 +664,22 @@ let find_definition env rs = (* else look for a global function *) find_global_definition env.known rs +(* Assuming the real is given in pow2 and pow5 *) +let compute_fraction { Number.rv_sig = i; Number.rv_pow2 = p2; Number.rv_pow5 = p5 } = + let p2_val = BigInt.pow_int_pos_bigint 2 (BigInt.abs p2) in + let p5_val = BigInt.pow_int_pos_bigint 5 (BigInt.abs p5) in + let num = ref BigInt.one in + let den = ref BigInt.one in + num := BigInt.mul i !num; + if BigInt.ge p2 BigInt.zero then + num := BigInt.mul p2_val !num + else + den := BigInt.mul p2_val !den; + if BigInt.ge p5 BigInt.zero then + num := BigInt.mul p5_val !num + else + den := BigInt.mul p5_val !den; + (!num, !den) (* evaluate expressions *) let rec eval_expr env (e : expr) : result = @@ -492,7 +695,19 @@ let rec eval_expr env (e : expr) : result = Normal v with Not_found -> assert false (* Irred e ? *) end - | Econst c -> Normal (Vnum (big_int_of_const c)) + | Econst (Number.ConstInt c) -> + Normal (Vnum (big_int_of_const c)) + | Econst (Number.ConstReal r) -> + (* ConstReal can be float or real *) + let is_real ity = ity_equal ity ity_real in + if is_real e.e_ity then + let (p, q) = compute_fraction r.Number.rl_real in + let (sp, sq) = BigInt.to_string p, BigInt.to_string q in + try Normal (Vreal (real_from_fraction sp sq)) with + | Mlmpfr_wrapper.Not_Implemented -> raise CannotCompute + else + let s = Format.asprintf "%a" Number.print_constant (Number.ConstReal r) in + Normal (Vfloat (make_from_str s)) | Eexec (ce,cty) -> assert (cty.cty_args = []); assert (ce.c_cty.cty_args = []); @@ -741,9 +956,6 @@ and exec_call env rs args ity_result = rs.rs_name.Ident.id_string; raise CannotCompute - - - let eval_global_expr env km locals e = get_builtin_progs env; let env = { @@ -777,9 +989,16 @@ let eval_global_expr env km locals e = let res = eval_expr env e in res, global_env +let init_real real_param = + match real_param with + | None -> () + | Some (emin, emax, prec) -> + Big_real.init emin emax prec - -let eval_global_symbol env m fmt rs = +(* real_param is used in case of real computation for the user to give the + precision she wants to use *) +let eval_global_symbol ?real_param env m fmt rs = + init_real real_param; try match find_global_definition m.Pmodule.mod_known rs with | Function(locals,d) -> @@ -813,6 +1032,9 @@ let eval_global_symbol env m fmt rs = end end | _ -> assert false - with Not_found -> + with + | Not_found -> eprintf "Symbol '%s' has no definition.@." rs.rs_name.Ident.id_string; exit 1 + | CannotFind (l, n, s) -> + eprintf "Cannot find builtin symbol %s %s %s@." (List.fold_left (fun x y -> x^y) "" l) n s diff --git a/src/mlw/pinterp.mli b/src/mlw/pinterp.mli index 094e6ab7b..a9b9f2bba 100644 --- a/src/mlw/pinterp.mli +++ b/src/mlw/pinterp.mli @@ -33,4 +33,9 @@ val eval_global_expr: *) val eval_global_symbol: - Env.env -> Pmodule.pmodule -> Format.formatter -> Expr.rsymbol -> unit + ?real_param:(int*int*int) -> Env.env -> Pmodule.pmodule -> + Format.formatter -> Expr.rsymbol -> unit +(* [eval_global_symbol ?real_param env m fmt rs] evaluate global symbol [rs] of + module [m] in environment [env]. It prints output to [fmt]. + [real_param] is an optional argument used to give a precision on real + computation if necessary *) diff --git a/src/tools/why3execute.ml b/src/tools/why3execute.ml index cbd4e8619..1f9e811b3 100644 --- a/src/tools/why3execute.ml +++ b/src/tools/why3execute.ml @@ -78,8 +78,8 @@ let () = try Opt.iter do_input !opt_file with e when not (Debug.test_flag Debug.stack_trace) -> - eprintf "%a@." Exn_printer.exn_printer e; - exit 1 + eprintf "%a@." Exn_printer.exn_printer e; + exit 1 (* Local Variables: diff --git a/src/util/mlmpfr_dummy.ml b/src/util/mlmpfr_dummy.ml new file mode 100644 index 000000000..772daa065 --- /dev/null +++ b/src/util/mlmpfr_dummy.ml @@ -0,0 +1,59 @@ +(* Dummy implementation for mpfr that always raise Not_implemented *) + +exception Not_Implemented + +type mpfr_rnd_t = + | To_Nearest + | Toward_Zero + | Toward_Plus_Infinity + | Toward_Minus_Infinity + | Away_From_Zero + | Faithful + +type sign = Positive | Negative + +type mpfr_float = unit + +let set_emin _ = () +let set_emax _ = () +let set_default_prec _ = () + +(* Elementary operations *) +let add ?rnd:_ ?prec:_ _ _ = raise Not_Implemented +let neg ?rnd:_ ?prec:_ _ = raise Not_Implemented +let mul ?rnd:_ ?prec:_ _ _ = raise Not_Implemented +let div ?rnd:_ ?prec:_ _ _ = raise Not_Implemented +let sqrt ?rnd:_ ?prec:_ _ = raise Not_Implemented +let sub ?rnd:_ ?prec:_ _ _ = raise Not_Implemented +let abs ?rnd:_ ?prec:_ _ = raise Not_Implemented +let fma ?rnd:_ ?prec:_ _ _ _ = raise Not_Implemented +let rint ?rnd:_ ?prec:_ _ = raise Not_Implemented + +let min ?rnd:_ ?prec:_ _ _ = raise Not_Implemented +let max ?rnd:_ ?prec:_ _ _ = raise Not_Implemented + +let signbit _ = raise Not_Implemented + +let subnormalize ?rnd:_ _ = () (* Used outside of try *) + +let make_from_str ?prec:_ ?rnd:_ ?base:_ _ = () +let make_from_int ?prec:_ ?rnd:_ _ = () +let make_zero ?prec:_ _ = () + +let get_formatted_str ?rnd:_ ?base:_ ?size:_ _ = "Error: MPFR not configured" + +(* Comparison functions *) + +let greater_p _ _ = raise Not_Implemented +let greaterequal_p _ _ = raise Not_Implemented +let less_p _ _ = raise Not_Implemented +let lessequal_p _ _ = raise Not_Implemented +let equal_p _ _ = raise Not_Implemented +let lessgreater_p _ _ = raise Not_Implemented + +let zero_p _ = raise Not_Implemented +let nan_p _ = raise Not_Implemented +let inf_p _ = raise Not_Implemented + +(* Constants *) +let const_pi ?rnd:_ _ = () diff --git a/src/util/mlmpfr_real.ml b/src/util/mlmpfr_real.ml new file mode 100644 index 000000000..764f2d916 --- /dev/null +++ b/src/util/mlmpfr_real.ml @@ -0,0 +1,6 @@ +(* TODO This wrapper should eventually be removed ! *) + +(* Exception to be raised if mpfr is not installed *) +exception Not_Implemented + +include Mpfr diff --git a/src/util/mlmpfr_wrapper.mli b/src/util/mlmpfr_wrapper.mli new file mode 100644 index 000000000..28c6ecfaa --- /dev/null +++ b/src/util/mlmpfr_wrapper.mli @@ -0,0 +1,69 @@ +(* TODO This wrapper should eventually be removed ! + For documentation refer to the mlmpfr documentation. +*) + +(* Wrapper for mlMPFR: + MPFR installed -> implemented with mlMPFR + MPFR not installed -> dummy implementation +*) + +type mpfr_float + +(* Exception to be raised if mpfr is not installed *) +exception Not_Implemented + +val set_emin: int -> unit +val set_emax: int -> unit +val set_default_prec: int -> unit + +type mpfr_rnd_t = + | To_Nearest + | Toward_Zero + | Toward_Plus_Infinity + | Toward_Minus_Infinity + | Away_From_Zero + | Faithful + +type sign = Positive | Negative + +(* Elementary operations *) +val add: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float -> mpfr_float +val neg: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float +val mul: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float -> mpfr_float +val div: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float -> mpfr_float +val sqrt: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float +val sub: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float -> mpfr_float +val abs: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float +val fma: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float -> mpfr_float -> mpfr_float +val rint: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float + + +val min: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float -> mpfr_float +val max: ?rnd:mpfr_rnd_t -> ?prec:int -> mpfr_float -> mpfr_float -> mpfr_float + +val signbit: mpfr_float -> sign + +val subnormalize : ?rnd:mpfr_rnd_t -> mpfr_float -> mpfr_float + + +val make_from_str: ?prec:int -> ?rnd:mpfr_rnd_t -> ?base:int -> string -> mpfr_float +val make_from_int: ?prec:int -> ?rnd:mpfr_rnd_t -> int -> mpfr_float +val make_zero: ?prec:int -> sign -> mpfr_float + +val get_formatted_str : ?rnd:mpfr_rnd_t -> ?base:int -> ?size:int -> mpfr_float -> string + +(* Comparison functions *) + +val greater_p : mpfr_float -> mpfr_float -> bool +val greaterequal_p : mpfr_float -> mpfr_float -> bool +val less_p : mpfr_float -> mpfr_float -> bool +val lessequal_p : mpfr_float -> mpfr_float -> bool +val equal_p : mpfr_float -> mpfr_float -> bool +val lessgreater_p : mpfr_float -> mpfr_float -> bool + +val zero_p: mpfr_float -> bool +val nan_p : mpfr_float -> bool +val inf_p : mpfr_float -> bool + +(* Constants *) +val const_pi: ?rnd:mpfr_rnd_t -> int -> mpfr_float diff --git a/stdlib/ieee_float.mlw b/stdlib/ieee_float.mlw index eb5bed068..dab1fad10 100644 --- a/stdlib/ieee_float.mlw +++ b/stdlib/ieee_float.mlw @@ -127,11 +127,11 @@ module GenericFloat predicate is_normal t predicate is_subnormal t - predicate is_zero t - predicate is_infinite t - predicate is_nan t - predicate is_positive t - predicate is_negative t + val predicate is_zero t + val predicate is_infinite t + val predicate is_nan t + val predicate is_positive t + val predicate is_negative t (** helper predicate for zeros, normals and subnormals. not defined so that the axiomatisation below can use it without talking about diff --git a/stdlib/real.mlw b/stdlib/real.mlw index d3981fe9c..bcfa9504e 100644 --- a/stdlib/real.mlw +++ b/stdlib/real.mlw @@ -202,7 +202,7 @@ module Square function sqr (x : real) : real = x * x - function sqrt real : real + val function sqrt real : real axiom Sqrt_positive: forall x:real. x >= 0.0 -> sqrt x >= 0.0 @@ -335,7 +335,7 @@ module Trigonometry axiom Cos_0: cos 0.0 = 1.0 axiom Sin_0: sin 0.0 = 0.0 - constant pi : real + val constant pi : real axiom Pi_double_precision_bounds: 0x1.921fb54442d18p+1 < pi < 0x1.921fb54442d19p+1 -- 2.24.1