more tests, several bug fixes

bench/interp/real.mlw

@@ -7,9 +7,21 @@ module R
   exception BenchFailure
 
 
-  let bench0 () =
-    let x : real = 1.0 / 3.0 in
-    (x, 3.0 * x)
+  let test0 () =
+     let x:real = 42.0 in
+     let y:real = 1.5 in
+     let z:real = 0.1 in
+     (x,y,z)
+
+  let test1 () =
+    let x : real = 0.1 in
+    let y : real = 10.0 * x in
+    (x, y)
+
+  let test2 () =
+    let x : real = 3.0 in
+    let y : real = 1.0 / x in
+    (x, y, 3.0 * y)
 
   let bench1 ()
     (* Tries to calculate sqrt(2) *)
@@ -24,9 +36,12 @@ module R
     !i
 
   let bench2 ()
-(*    raises { BenchFailure -> false }*)
-    ensures { result = 4.0 } =
-    sqrt (16.0)
+    raises { BenchFailure -> false }
+    ensures { result = 4.0 }
+  =
+    let r = sqrt 16.0 in
+    if r <> 4.0 then raise BenchFailure;
+    r
 
   use real.Trigonometry
 

src/mlw/big_real.ml

@@ -7,13 +7,15 @@ type real = mpfr_float * mpfr_float
 (* computationally, a real is represented as an interval of two floating-point numbers.
    such an interval `[a;b]` represents the set of real numbers between `a` and `b` *)
 
-let init, set_exponents, get_prec =
+let init, set_exponents, get_prec, get_zero, get_one =
 (*
    By default, for approximating real numbers, let's use binary128 floats
 *)
   let emin = ref (-16493) in
   let emax = ref 16384 in
   let prec = ref 113 in
+  let zero = ref (make_zero ~prec:!prec Positive) in
+  let one = ref (make_from_int ~prec:!prec ~rnd:Toward_Minus_Infinity 1) in
   (fun emin_i emax_i prec_i ->
     emin := emin_i;
     emax := emax_i;
@@ -21,7 +23,16 @@ let init, set_exponents, get_prec =
   (fun () ->
     set_emin !emin;
     set_emax !emax),
-  (fun () -> !prec)
+  (fun () -> !prec),
+  (fun () -> !zero),
+  (fun () -> !one)
+
+let print_float fmt x =
+  Format.fprintf fmt "%s" (get_formatted_str ~base:10 x)
+
+let print_real fmt (x, y) =
+  Format.fprintf fmt "[%a, %a]" print_float x print_float y
+
 
 let add (xmin, xmax) (ymin, ymax) =
   (* Exponents can be changed if floats occur in the code. *)
@@ -58,17 +69,18 @@ let mul (xmin, xmax) (ymin, ymax) =
 let inv (xmin, xmax) =
   set_exponents ();
   let prec = get_prec () in
+  let zero = get_zero () in
   (* If 0 is inside the interval we cannot compute the expression *)
-  if signbit xmin <> signbit xmax then
+  if lessequal_p xmin zero && lessequal_p zero xmax then
     raise Undetermined
   else
-    let one = make_from_int ~prec ~rnd:Toward_Minus_Infinity 1 in
+    let one = get_one () 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 res_min = min ~prec ~rnd:Toward_Minus_Infinity 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
+    let res_max = max  ~prec ~rnd:Toward_Plus_Infinity inv1_max inv2_max in
     (res_min, res_max)
 
 let div x y =
@@ -77,7 +89,8 @@ let div x y =
 let sqrt (xmin, xmax) =
   set_exponents ();
   let prec = get_prec() in
-  if lessequal_p (make_zero ~prec Positive) xmin then
+  let zero = get_zero () in
+  if lessequal_p zero 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)
@@ -86,44 +99,36 @@ let sqrt (xmin, xmax) =
 
 let real_from_str s =
   let prec = get_prec () in
-  let n = make_from_str ~prec ~base:10 s in
-  (n, n)
+  let n1 = make_from_str ~prec ~base:10 ~rnd:Toward_Minus_Infinity s in
+  let n2 = make_from_str ~prec ~base:10 ~rnd:Toward_Plus_Infinity s in
+  (n1, n2)
 
 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
+let lt (x1,x2) (y1,y2) =
+  if less_p x2 y1 then true else
+    if lessequal_p y2 x1 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
+let le (x1,x2) (y1,y2) =
+  if lessequal_p x2 y1 then true else
+    if less_p y2 x1 then false else
       raise Undetermined
 
-let pi =
+let gt x y = lt y x
+
+let ge x y = le y x
+
+let pi () =
   let prec = get_prec () in
-  const_pi prec, const_pi prec
+  (const_pi ~rnd:Toward_Minus_Infinity prec,
+   const_pi ~rnd:Toward_Plus_Infinity prec)

src/mlw/big_real.mli

@@ -22,6 +22,8 @@ val sqrt: real -> real
 val eq: real -> real -> bool
 val lt: real -> real -> bool
 val le: real -> real -> bool
+val gt: real -> real -> bool
+val ge: real -> real -> bool
 
 (* Constants *)
-val pi: real
+val pi: unit -> real

src/mlw/pinterp.ml

@@ -484,7 +484,7 @@ let built_in_modules =
   in
   let real_trigo_module =
     ["real"], "Trigonometry", [],
-    ["pi", eval_real Modeconst Big_real.pi;
+    ["pi", eval_real Modeconst (Big_real.pi());
     ]
   in
   bool_module :: (int_modules @ [array_module;