diff --git a/thys/IEEE_Floating_Point/IEEE.thy b/thys/IEEE_Floating_Point/IEEE.thy
--- a/thys/IEEE_Floating_Point/IEEE.thy
+++ b/thys/IEEE_Floating_Point/IEEE.thy
@@ -1,448 +1,457 @@
 (* Formalization of IEEE-754 Standard for binary floating-point arithmetic *)
-(* Author: Lei Yu, University of Cambridge *)
+(* Author: Lei Yu, University of Cambridge
+
+  Contrib: Peter Lammich: fixed wrong sign handling in fmadd
+
+*)
 
 section \<open>Specification of the IEEE standard\<close>
 
 theory IEEE
   imports
     "HOL-Library.Float"
     Word_Lib.Word_Lemmas
 begin
 
 typedef (overloaded) ('e::len, 'f::len) float = "UNIV::(1 word \<times> 'e word \<times> 'f word) set"
   by auto
 
 setup_lifting type_definition_float
 
 syntax "_float" :: "type \<Rightarrow> type \<Rightarrow> type" ("'(_, _') float")
 text \<open>parse \<open>('a, 'b) float\<close> as ('a::len, 'b::len) float.\<close>
 
 parse_translation \<open>
   let
     fun float t u = Syntax.const @{type_syntax float} $ t $ u;
     fun len_tr u =
       (case Term_Position.strip_positions u of
         v as Free (x, _) =>
           if Lexicon.is_tid x then
             (Syntax.const @{syntax_const "_ofsort"} $ v $
               Syntax.const @{class_syntax len})
           else u
       | _ => u)
     fun len_float_tr [t, u] =
       float (len_tr t) (len_tr u)
   in
     [(@{syntax_const "_float"}, K len_float_tr)]
   end
 \<close>
 
 subsection \<open>Derived parameters for floating point formats\<close>
 
 definition wordlength :: "('e, 'f) float itself \<Rightarrow> nat"
   where "wordlength x = LENGTH('e) + LENGTH('f) + 1"
 
 definition bias :: "('e, 'f) float itself \<Rightarrow> nat"
   where "bias x = 2^(LENGTH('e) - 1) - 1"
 
 definition emax :: "('e, 'f) float itself \<Rightarrow> nat"
   where "emax x = unat (- 1::'e word)"
 
 abbreviation fracwidth::"('e, 'f) float itself \<Rightarrow> nat" where
   "fracwidth _ \<equiv> LENGTH('f)"
 
 subsection \<open>Predicates for the four IEEE formats\<close>
 
 definition is_single :: "('e, 'f) float itself \<Rightarrow> bool"
   where "is_single x \<longleftrightarrow> LENGTH('e) = 8 \<and> wordlength x = 32"
 
 definition is_double :: "('e, 'f) float itself \<Rightarrow> bool"
   where "is_double x \<longleftrightarrow> LENGTH('e) = 11 \<and> wordlength x = 64"
 
 definition is_single_extended :: "('e, 'f) float itself \<Rightarrow> bool"
   where "is_single_extended x \<longleftrightarrow> LENGTH('e) \<ge> 11 \<and> wordlength x \<ge> 43"
 
 definition is_double_extended :: "('e, 'f) float itself \<Rightarrow> bool"
   where "is_double_extended x \<longleftrightarrow> LENGTH('e) \<ge> 15 \<and> wordlength x \<ge> 79"
 
 
 subsection \<open>Extractors for fields\<close>
 
 lift_definition sign::"('e, 'f) float \<Rightarrow> nat" is
   "\<lambda>(s::1 word, _::'e word, _::'f word). unat s" .
 
 lift_definition exponent::"('e, 'f) float \<Rightarrow> nat" is
   "\<lambda>(_, e::'e word, _). unat e" .
 
 lift_definition fraction::"('e, 'f) float \<Rightarrow> nat" is
   "\<lambda>(_, _, f::'f word). unat f" .
 
 abbreviation "real_of_word x \<equiv> real (unat x)"
 
 lift_definition valof :: "('e, 'f) float \<Rightarrow> real"
 is "\<lambda>(s, e, f).
     let x = (TYPE(('e, 'f) float)) in
     (if e = 0
      then (-1::real)^(unat s) * (2 / (2^bias x)) * (real_of_word f/2^(LENGTH('f)))
      else (-1::real)^(unat s) * ((2^(unat e)) / (2^bias x)) * (1 + real_of_word f/2^LENGTH('f)))"
   .
 
 subsection \<open>Partition of numbers into disjoint classes\<close>
 
 definition is_nan :: "('e, 'f) float \<Rightarrow> bool"
   where "is_nan a \<longleftrightarrow> exponent a = emax TYPE(('e, 'f)float) \<and> fraction a \<noteq> 0"
 
 definition is_infinity :: "('e, 'f) float \<Rightarrow> bool"
   where "is_infinity a \<longleftrightarrow> exponent a = emax TYPE(('e, 'f)float) \<and> fraction a = 0"
 
 definition is_normal :: "('e, 'f) float \<Rightarrow> bool"
   where "is_normal a \<longleftrightarrow> 0 < exponent a \<and> exponent a < emax TYPE(('e, 'f)float)"
 
 definition is_denormal :: "('e, 'f) float \<Rightarrow> bool"
   where "is_denormal a \<longleftrightarrow> exponent a = 0 \<and> fraction a \<noteq> 0"
 
 definition is_zero :: "('e, 'f) float \<Rightarrow> bool"
   where "is_zero a \<longleftrightarrow> exponent a = 0 \<and> fraction a = 0"
 
 definition is_finite :: "('e, 'f) float \<Rightarrow> bool"
   where "is_finite a \<longleftrightarrow> (is_normal a \<or> is_denormal a \<or> is_zero a)"
 
 
 subsection \<open>Special values\<close>
 
 lift_definition plus_infinity :: "('e, 'f) float" ("\<infinity>") is "(0, - 1, 0)" .
 
 lift_definition topfloat :: "('e, 'f) float" is "(0, - 2, 2^LENGTH('f) - 1)" .
 
 instantiation float::(len, len) zero begin
 
 lift_definition zero_float :: "('e, 'f) float" is "(0, 0, 0)" .
 
 instance proof qed
 
 end
 
 subsection \<open>Negation operation on floating point values\<close>
 
 instantiation float::(len, len) uminus begin
 lift_definition uminus_float :: "('e, 'f) float \<Rightarrow> ('e, 'f) float" is  "\<lambda>(s, e, f). (1 - s, e, f)" .
 instance proof qed
 end
 
 abbreviation (input) "minus_zero \<equiv> - (0::('e, 'f)float)"
 abbreviation (input) "minus_infinity \<equiv> - \<infinity>"
 abbreviation (input) "bottomfloat \<equiv> - topfloat"
 
 
 subsection \<open>Real number valuations\<close>
 
 text \<open>The largest value that can be represented in floating point format.\<close>
 definition largest :: "('e, 'f) float itself \<Rightarrow> real"
   where "largest x = (2^(emax x - 1) / 2^bias x) * (2 - 1/(2^fracwidth x))"
 
 text \<open>Threshold, used for checking overflow.\<close>
 definition threshold :: "('e, 'f) float itself \<Rightarrow> real"
   where "threshold x = (2^(emax x - 1) / 2^bias x) * (2 - 1/(2^(Suc(fracwidth x))))"
 
 text \<open>Unit of least precision.\<close>
 
 lift_definition one_lp::"('e ,'f) float \<Rightarrow> ('e ,'f) float" is "\<lambda>(s, e, f). (0, e::'e word, 1)" .
 lift_definition zero_lp::"('e ,'f) float \<Rightarrow> ('e ,'f) float" is "\<lambda>(s, e, f). (0, e::'e word, 0)" .
 
 definition ulp :: "('e, 'f) float \<Rightarrow> real" where "ulp a = valof (one_lp a) - valof (zero_lp a)"
 
 text \<open>Enumerated type for rounding modes.\<close>
 datatype roundmode = To_nearest | float_To_zero | To_pinfinity | To_ninfinity
 
 text \<open>Characterization of best approximation from a set of abstract values.\<close>
 definition "is_closest v s x a \<longleftrightarrow> a \<in> s \<and> (\<forall>b. b \<in> s \<longrightarrow> \<bar>v a - x\<bar> \<le> \<bar>v b - x\<bar>)"
 
 text \<open>Best approximation with a deciding preference for multiple possibilities.\<close>
 definition "closest v p s x =
   (SOME a. is_closest v s x a \<and> ((\<exists>b. is_closest v s x b \<and> p b) \<longrightarrow> p a))"
 
 
 subsection \<open>Rounding\<close>
 
 fun round :: "roundmode \<Rightarrow> real \<Rightarrow> ('e ,'f) float"
 where
   "round To_nearest y =
    (if y \<le> - threshold TYPE(('e ,'f) float) then minus_infinity
     else if y \<ge> threshold TYPE(('e ,'f) float) then plus_infinity
     else closest (valof) (\<lambda>a. even (fraction a)) {a. is_finite a} y)"
 | "round float_To_zero y =
    (if y < - largest TYPE(('e ,'f) float) then bottomfloat
     else if y > largest TYPE(('e ,'f) float) then topfloat
     else closest (valof) (\<lambda>a. True) {a. is_finite a \<and> \<bar>valof a\<bar> \<le> \<bar>y\<bar>} y)"
 | "round To_pinfinity y =
    (if y < - largest TYPE(('e ,'f) float) then bottomfloat
     else if y > largest TYPE(('e ,'f) float) then plus_infinity
     else closest (valof) (\<lambda>a. True) {a. is_finite a \<and> valof a \<ge> y} y)"
 | "round To_ninfinity y =
    (if y < - largest TYPE(('e ,'f) float) then minus_infinity
     else if y > largest TYPE(('e ,'f) float) then topfloat
     else closest (valof) (\<lambda>a. True) {a. is_finite a \<and> valof a \<le> y} y)"
 
 text \<open>Rounding to integer values in floating point format.\<close>
 
 definition is_integral :: "('e ,'f) float \<Rightarrow> bool"
   where "is_integral a \<longleftrightarrow> is_finite a \<and> (\<exists>n::nat. \<bar>valof a\<bar> = real n)"
 
 fun intround :: "roundmode \<Rightarrow> real \<Rightarrow> ('e ,'f) float"
 where
   "intround To_nearest y =
     (if y \<le> - threshold TYPE(('e ,'f) float) then minus_infinity
      else if y \<ge> threshold TYPE(('e ,'f) float) then plus_infinity
      else closest (valof) (\<lambda>a. (\<exists>n::nat. even n \<and> \<bar>valof a\<bar> = real n)) {a. is_integral a} y)"
 |"intround float_To_zero y =
     (if y < - largest TYPE(('e ,'f) float) then bottomfloat
      else if y > largest TYPE(('e ,'f) float) then topfloat
      else closest (valof) (\<lambda>x. True) {a. is_integral a \<and> \<bar>valof a\<bar> \<le> \<bar>y\<bar>} y)"
 |"intround To_pinfinity y =
     (if y < - largest TYPE(('e ,'f) float) then bottomfloat
      else if y > largest TYPE(('e ,'f) float) then plus_infinity
      else closest (valof) (\<lambda>x. True) {a. is_integral a \<and> valof a \<ge> y} y)"
 |"intround To_ninfinity y =
     (if y < - largest TYPE(('e ,'f) float) then minus_infinity
      else if y > largest TYPE(('e ,'f) float) then topfloat
      else closest (valof) (\<lambda>x. True) {a. is_integral a \<and> valof a \<ge> y} y)"
 
 text \<open>Round, choosing between -0.0 or +0.0\<close>
 
 definition float_round::"roundmode \<Rightarrow> bool \<Rightarrow> real \<Rightarrow> ('e, 'f) float"
   where "float_round mode toneg r =
       (let x = round mode r in
          if is_zero x
             then if toneg
                     then minus_zero
                  else 0
          else x)"
 
 text \<open>Non-standard of NaN.\<close>
 definition some_nan :: "('e ,'f) float"
   where "some_nan = (SOME a. is_nan a)"
 
 text \<open>Coercion for signs of zero results.\<close>
 definition zerosign :: "nat \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "zerosign s a =
     (if is_zero a then (if s = 0 then 0 else - 0) else a)"
 
 text \<open>Remainder operation.\<close>
 definition rem :: "real \<Rightarrow> real \<Rightarrow> real"
   where "rem x y =
     (let n = closest id (\<lambda>x. \<exists>n::nat. even n \<and> \<bar>x\<bar> = real n) {x. \<exists>n :: nat. \<bar>x\<bar> = real n} (x / y)
      in x - n * y)"
 
 definition frem :: "roundmode \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "frem m a b =
     (if is_nan a \<or> is_nan b \<or> is_infinity a \<or> is_zero b then some_nan
      else zerosign (sign a) (round m (rem (valof a) (valof b))))"
 
 
 subsection \<open>Definitions of the arithmetic operations\<close>
 
 definition fintrnd :: "roundmode \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "fintrnd m a =
     (if is_nan a then (some_nan)
      else if is_infinity a then a
      else zerosign (sign a) (intround m (valof a)))"
 
 definition fadd :: "roundmode \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "fadd m a b =
     (if is_nan a \<or> is_nan b \<or> (is_infinity a \<and> is_infinity b \<and> sign a \<noteq> sign b)
      then some_nan
      else if (is_infinity a) then a
      else if (is_infinity b) then b
      else
       zerosign
         (if is_zero a \<and> is_zero b \<and> sign a = sign b then sign a
          else if m = To_ninfinity then 1 else 0)
         (round m (valof a + valof b)))"
 
 definition fsub :: "roundmode \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "fsub m a b =
     (if is_nan a \<or> is_nan b \<or> (is_infinity a \<and> is_infinity b \<and> sign a = sign b)
      then some_nan
      else if is_infinity a then a
      else if is_infinity b then - b
      else
       zerosign
         (if is_zero a \<and> is_zero b \<and> sign a \<noteq> sign b then sign a
          else if m = To_ninfinity then 1 else 0)
         (round m (valof a - valof b)))"
 
 definition fmul :: "roundmode \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "fmul m a b =
     (if is_nan a \<or> is_nan b \<or> (is_zero a \<and> is_infinity b) \<or> (is_infinity a \<and> is_zero b)
      then some_nan
      else if is_infinity a \<or> is_infinity b
      then (if sign a = sign b then plus_infinity else minus_infinity)
      else zerosign (if sign a = sign b then 0 else 1 ) (round m (valof a * valof b)))"
 
 definition fdiv :: "roundmode \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "fdiv m a b =
     (if is_nan a \<or> is_nan b \<or> (is_zero a \<and> is_zero b) \<or> (is_infinity a \<and> is_infinity b)
      then some_nan
      else if is_infinity a \<or> is_zero b
      then (if sign a = sign b then plus_infinity else minus_infinity)
      else if is_infinity b
      then (if sign a = sign b then 0 else - 0)
      else zerosign (if sign a = sign b then 0 else 1) (round m (valof a / valof b)))"
 
 definition fsqrt :: "roundmode \<Rightarrow> ('e ,'f) float \<Rightarrow> ('e ,'f) float"
   where "fsqrt m a =
     (if is_nan a then some_nan
      else if is_zero a \<or> is_infinity a \<and> sign a = 0 then a
      else if sign a = 1 then some_nan
      else zerosign (sign a) (round m (sqrt (valof a))))"
 
 definition fmul_add :: "roundmode \<Rightarrow> ('t ,'w) float \<Rightarrow> ('t ,'w) float \<Rightarrow> ('t ,'w) float \<Rightarrow> ('t ,'w) float"
-  where "fmul_add mode x y z =
-      (let signP = if sign x = sign y then 0 else 1 in
-      let infP = is_infinity x  \<or> is_infinity y
-      in
-         if is_nan x \<or> is_nan y \<or> is_nan z then some_nan
-         else if is_infinity x \<and> is_zero y \<or>
-                 is_zero x \<and> is_infinity y \<or>
-                 is_infinity z \<and> infP \<and> signP \<noteq> sign z then
-            some_nan
-         else if is_infinity z \<and> (sign z = 0) \<or> infP \<and> (signP = 0)
-            then plus_infinity
-         else if is_infinity z \<and> (sign z = 1) \<or> infP \<and> (signP = 1)
-            then minus_infinity
-         else
-            let r1 = valof x * valof y;
-                r2 = valof z
-            in
-              float_round mode
-                (if (r1 = 0) \<and> (r2 = 0) \<and> (signP = sign z) then
-                   signP = 1
-                 else mode = To_ninfinity) (r1 + r2))"
+  where "fmul_add mode x y z = (let 
+    signP = if sign x = sign y then 0 else 1;
+    infP = is_infinity x  \<or> is_infinity y
+  in
+    if is_nan x \<or> is_nan y \<or> is_nan z then some_nan
+    else if is_infinity x \<and> is_zero y \<or>
+           is_zero x \<and> is_infinity y \<or>
+           is_infinity z \<and> infP \<and> signP \<noteq> sign z 
+      then some_nan
+    else if is_infinity z \<and> (sign z = 0) \<or> infP \<and> (signP = 0)
+      then plus_infinity
+    else if is_infinity z \<and> (sign z = 1) \<or> infP \<and> (signP = 1)
+      then minus_infinity
+    else let 
+      r1 = valof x * valof y;
+      r2 = valof z;
+      r = r1+r2
+    in 
+      if r=0 then ( \<comment> \<open>Exact Zero Case. Same sign rules as for add apply. \<close>
+        if r1=0 \<and> r2=0 \<and> signP=sign z then zerosign signP 0
+        else if mode = To_ninfinity then -0
+        else 0
+      ) else ( \<comment> \<open>Not exactly zero: Rounding has sign of exact value, even if rounded val is zero\<close>
+        zerosign (if r<0 then 1 else 0) (round mode r)
+      )
+  )"
 
 
 subsection \<open>Comparison operations\<close>
 
 datatype ccode = Gt | Lt | Eq | Und
 
 definition fcompare :: "('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> ccode"
   where "fcompare a b =
     (if is_nan a \<or> is_nan b then Und
      else if is_infinity a \<and> sign a = 1
      then (if is_infinity b \<and> sign b = 1 then Eq else Lt)
      else if is_infinity a \<and> sign a = 0
      then (if is_infinity b \<and> sign b = 0 then Eq else Gt)
      else if is_infinity b \<and> sign b = 1 then Gt
      else if is_infinity b \<and> sign b = 0 then Lt
      else if valof a < valof b then Lt
      else if valof a = valof b then Eq
      else Gt)"
 
 definition flt :: "('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> bool"
   where "flt a b \<longleftrightarrow> fcompare a b = Lt"
 
 definition fle :: "('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> bool"
   where "fle a b \<longleftrightarrow> fcompare a b = Lt \<or> fcompare a b = Eq"
 
 definition fgt :: "('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> bool"
   where "fgt a b \<longleftrightarrow> fcompare a b = Gt"
 
 definition fge :: "('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> bool"
   where "fge a b \<longleftrightarrow> fcompare a b = Gt \<or> fcompare a b = Eq"
 
 definition feq :: "('e ,'f) float \<Rightarrow> ('e ,'f) float \<Rightarrow> bool"
   where "feq a b \<longleftrightarrow> fcompare a b = Eq"
 
 
 section \<open>Specify float to be double  precision and round to even\<close>
 
 instantiation float :: (len, len) plus
 begin
 
 definition plus_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "a + b = fadd To_nearest a b"
 
 instance ..
 
 end
 
 instantiation float :: (len, len) minus
 begin
 
 definition minus_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "a - b = fsub To_nearest a b"
 
 instance ..
 
 end
 
 instantiation float :: (len, len) times
 begin
 
 definition times_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "a * b = fmul To_nearest a b"
 
 instance ..
 
 end
 
 instantiation float :: (len, len) one
 begin
 
 lift_definition one_float :: "('a, 'b) float" is "(0, 2^(LENGTH('a) - 1) - 1, 0)" .
 
 instance ..
 
 end
 
 instantiation float :: (len, len) inverse
 begin
 
 definition divide_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "a div b = fdiv To_nearest a b"
 
 definition inverse_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "inverse_float a = fdiv To_nearest 1 a"
 
 instance ..
 
 end
 
 definition float_rem :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "float_rem a b = frem To_nearest a b"
 
 definition float_sqrt :: "('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "float_sqrt a = fsqrt To_nearest a"
 
 definition ROUNDFLOAT ::"('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "ROUNDFLOAT a = fintrnd To_nearest a"
 
 
 instantiation float :: (len, len) ord
 begin
 
 definition less_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> bool"
   where "a < b \<longleftrightarrow> flt a b"
 
 definition less_eq_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> bool"
   where "a \<le> b \<longleftrightarrow> fle a b"
 
 instance ..
 
 end
 
 definition float_eq :: "('a, 'b) float \<Rightarrow> ('a, 'b) float \<Rightarrow> bool"  (infixl "\<doteq>" 70)
   where "float_eq a b = feq a b"
 
 instantiation float :: (len, len) abs
 begin
 
 definition abs_float :: "('a, 'b) float \<Rightarrow> ('a, 'b) float"
   where "abs_float a = (if sign a = 0 then a else - a)"
 
 instance ..
 end
 
 text \<open>The \<open>1 + \<epsilon>\<close> property.\<close>
 definition normalizes :: "_ itself \<Rightarrow> real \<Rightarrow> bool"
   where "normalizes float_format x =
     (1/ (2::real)^(bias float_format - 1) \<le> \<bar>x\<bar> \<and> \<bar>x\<bar> < threshold float_format)"
 
 end