diff options
| author | Pierre Roux | 2019-04-04 00:14:47 +0200 |
|---|---|---|
| committer | Pierre Roux | 2019-11-01 10:21:03 +0100 |
| commit | dca0135a263717b3a1a1d7c4f054f039dc08109e (patch) | |
| tree | 3f2ab5ea79e084a9c72e1376d5399ad4a62cb771 /kernel/float64.ml | |
| parent | 3e0db1b645a8653c62b8b5a4978e6d8fbbe9a9cc (diff) | |
Make primitive float work on x86_32
Flag -fexcess-precision=standard is not enough on x86_32
where -msse2 -mfpmath=sse is required (-msse is not enough)
to avoid double rounding issues in the VM.
Most floating-point operation are now implemented in C because OCaml
is suffering double rounding issues on x86_32 with 80 bits extended
precision registers used for floating-point values, causing double
rounding making floating-point arithmetic incorrect with respect to
its specification.
Add a runtime test for double roundings.
Diffstat (limited to 'kernel/float64.ml')
| -rw-r--r-- | kernel/float64.ml | 52 |
1 files changed, 24 insertions, 28 deletions
diff --git a/kernel/float64.ml b/kernel/float64.ml index 5160d32892..55ad472ea9 100644 --- a/kernel/float64.ml +++ b/kernel/float64.ml @@ -70,20 +70,20 @@ let classify x = | FP_nan -> NaN [@@ocaml.inline always] -let mul x y = x *. y -[@@ocaml.inline always] +external mul : float -> float -> float = "coq_fmul_byte" "coq_fmul" +[@@unboxed] [@@noalloc] -let add x y = x +. y -[@@ocaml.inline always] +external add : float -> float -> float = "coq_fadd_byte" "coq_fadd" +[@@unboxed] [@@noalloc] -let sub x y = x -. y -[@@ocaml.inline always] +external sub : float -> float -> float = "coq_fsub_byte" "coq_fsub" +[@@unboxed] [@@noalloc] -let div x y = x /. y -[@@ocaml.inline always] +external div : float -> float -> float = "coq_fdiv_byte" "coq_fdiv" +[@@unboxed] [@@noalloc] -let sqrt x = sqrt x -[@@ocaml.inline always] +external sqrt : float -> float = "coq_fsqrt_byte" "coq_fsqrt" +[@@unboxed] [@@noalloc] let of_int63 x = Uint63.to_float x [@@ocaml.inline always] @@ -111,25 +111,11 @@ let frshiftexp f = let ldshiftexp f e = ldexp f (snd (Uint63.to_int2 e) - eshift) [@@ocaml.inline always] -let eta_float = ldexp 1. (-1074) (* smallest positive float (subnormal) *) +external next_up : float -> float = "coq_next_up_byte" "coq_next_up" +[@@unboxed] [@@noalloc] -let next_up f = - match classify_float f with - | FP_nan -> f - | FP_infinite -> if 0. < f then f else -.max_float - | FP_zero | FP_subnormal -> - let f = f +. eta_float in - if f = 0. then -0. else f (* or next_down may return -0. *) - | FP_normal -> - let f, e = frexp f in - if 0. < f || f <> -0.5 || e = -1021 then - ldexp (f +. epsilon_float /. 2.) e - else - ldexp (-0.5 +. epsilon_float /. 4.) e -[@@ocaml.inline always] - -let next_down f = -.(next_up (-.f)) -[@@ocaml.inline always] +external next_down : float -> float = "coq_next_down_byte" "coq_next_down" +[@@unboxed] [@@noalloc] let equal f1 f2 = match classify_float f1 with @@ -153,3 +139,13 @@ let total_compare f1 f2 = let is_float64 t = Obj.tag t = Obj.double_tag [@@ocaml.inline always] + +(*** Test at runtime that no harmful double rounding seems to + be performed with an intermediate 80 bits representation (x87). *) +let () = + let b = ldexp 1. 53 in + let s = add 1. (ldexp 1. (-52)) in + if add b s <= b || add b 1. <> b then + failwith "Detected double rounding due to the use of intermediate \ + 80 bits floating-point representation. Use of Float is \ + thus unsafe." |