# HG changeset patch # User hoelzl # Date 1233830715 -3600 # Node ID e15b74577368db01c16659df84a63779c0c52e51 # Parent c56a5571f60a616214d24244a0727694fc48715a Added new Float theory and moved old Library/Float.thy to ComputeFloat diff -r c56a5571f60a -r e15b74577368 src/HOL/IsaMakefile --- a/src/HOL/IsaMakefile Thu Feb 05 11:34:42 2009 +0100 +++ b/src/HOL/IsaMakefile Thu Feb 05 11:45:15 2009 +0100 @@ -335,7 +335,7 @@ Library/Mapping.thy Library/Numeral_Type.thy Library/Reflection.thy \ Library/Boolean_Algebra.thy Library/Countable.thy \ Library/RBT.thy Library/Univ_Poly.thy \ - Library/Enum.thy Library/Float.thy $(SRC)/Tools/float.ML $(SRC)/HOL/Tools/float_arith.ML \ + Library/Enum.thy Library/Float.thy $(SRC)/Tools/float.ML \ Library/reify_data.ML Library/reflection.ML @cd Library; $(ISABELLE_TOOL) usedir $(OUT)/HOL Library @@ -883,6 +883,7 @@ $(SRC)/Tools/Compute_Oracle/am_ghc.ML \ $(SRC)/Tools/Compute_Oracle/am_sml.ML \ $(SRC)/Tools/Compute_Oracle/compute.ML \ + Tools/ComputeFloat.thy Tools/float_arith.ML \ Matrix/Matrix.thy Matrix/SparseMatrix.thy Matrix/LP.thy \ Matrix/document/root.tex Matrix/ROOT.ML Matrix/cplex/Cplex.thy \ Matrix/cplex/CplexMatrixConverter.ML Matrix/cplex/Cplex_tools.ML \ diff -r c56a5571f60a -r e15b74577368 src/HOL/Library/Float.thy --- a/src/HOL/Library/Float.thy Thu Feb 05 11:34:42 2009 +0100 +++ b/src/HOL/Library/Float.thy Thu Feb 05 11:45:15 2009 +0100 @@ -1,30 +1,56 @@ -(* Title: HOL/Real/Float.thy - Author: Steven Obua -*) - -header {* Floating Point Representation of the Reals *} - +(* Title: HOL/Library/Float.thy + * Author: Steven Obua 2008 + * Johannes Hölzl, TU Muenchen 2008 / 2009 + *) theory Float imports Complex_Main -uses "~~/src/Tools/float.ML" ("~~/src/HOL/Tools/float_arith.ML") begin definition pow2 :: "int \ real" where - "pow2 a = (if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a)))))" + [simp]: "pow2 a = (if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a)))))" + +datatype float = Float int int + +fun Ifloat :: "float \ real" where +"Ifloat (Float a b) = real a * pow2 b" -definition - float :: "int * int \ real" where - "float x = real (fst x) * pow2 (snd x)" +instantiation float :: zero begin +definition zero_float where "0 = Float 0 0" +instance .. +end + +instantiation float :: one begin +definition one_float where "1 = Float 1 0" +instance .. +end + +instantiation float :: number begin +definition number_of_float where "number_of n = Float n 0" +instance .. +end -lemma pow2_0[simp]: "pow2 0 = 1" -by (simp add: pow2_def) +fun mantissa :: "float \ int" where +"mantissa (Float a b) = a" + +fun scale :: "float \ int" where +"scale (Float a b) = b" + +lemma Ifloat_neg_exp: "e < 0 \ Ifloat (Float m e) = real m * inverse (2^nat (-e))" by auto +lemma Ifloat_nge0_exp: "\ 0 \ e \ Ifloat (Float m e) = real m * inverse (2^nat (-e))" by auto +lemma Ifloat_ge0_exp: "0 \ e \ Ifloat (Float m e) = real m * (2^nat e)" by auto -lemma pow2_1[simp]: "pow2 1 = 2" -by (simp add: pow2_def) +lemma Float_num[simp]: shows + "Ifloat (Float 1 0) = 1" and "Ifloat (Float 1 1) = 2" and "Ifloat (Float 1 2) = 4" and + "Ifloat (Float 1 -1) = 1/2" and "Ifloat (Float 1 -2) = 1/4" and "Ifloat (Float 1 -3) = 1/8" and + "Ifloat (Float -1 0) = -1" and "Ifloat (Float (number_of n) 0) = number_of n" + by auto -lemma pow2_neg: "pow2 x = inverse (pow2 (-x))" -by (simp add: pow2_def) +lemma pow2_0[simp]: "pow2 0 = 1" by simp +lemma pow2_1[simp]: "pow2 1 = 2" by simp +lemma pow2_neg: "pow2 x = inverse (pow2 (-x))" by simp + +declare pow2_def[simp del] lemma pow2_add1: "pow2 (1 + a) = 2 * (pow2 a)" proof - @@ -96,281 +122,43 @@ qed qed -lemma "float (a, e) + float (b, e) = float (a + b, e)" -by (simp add: float_def algebra_simps) +lemma float_components[simp]: "Float (mantissa f) (scale f) = f" by (cases f, auto) + +lemma float_split: "\ a b. x = Float a b" by (cases x, auto) -definition - int_of_real :: "real \ int" where - "int_of_real x = (SOME y. real y = x)" +lemma float_split2: "(\ a b. x \ Float a b) = False" by (auto simp add: float_split) + +lemma float_zero[simp]: "Ifloat (Float 0 e) = 0" by simp + +lemma abs_div_2_less: "a \ 0 \ a \ -1 \ abs((a::int) div 2) < abs a" +by arith -definition - real_is_int :: "real \ bool" where - "real_is_int x = (EX (u::int). x = real u)" +function normfloat :: "float \ float" where +"normfloat (Float a b) = (if a \ 0 \ even a then normfloat (Float (a div 2) (b+1)) else if a=0 then Float 0 0 else Float a b)" +by pat_completeness auto +termination by (relation "measure (nat o abs o mantissa)") (auto intro: abs_div_2_less) +declare normfloat.simps[simp del] -lemma real_is_int_def2: "real_is_int x = (x = real (int_of_real x))" -by (auto simp add: real_is_int_def int_of_real_def) +theorem normfloat[symmetric, simp]: "Ifloat f = Ifloat (normfloat f)" +proof (induct f rule: normfloat.induct) + case (1 a b) + have real2: "2 = real (2::int)" + by auto + show ?case + apply (subst normfloat.simps) + apply (auto simp add: float_zero) + apply (subst 1[symmetric]) + apply (auto simp add: pow2_add even_def) + done +qed -lemma float_transfer: "real_is_int ((real a)*(pow2 c)) \ float (a, b) = float (int_of_real ((real a)*(pow2 c)), b - c)" -by (simp add: float_def real_is_int_def2 pow2_add[symmetric]) +lemma pow2_neq_zero[simp]: "pow2 x \ 0" + by (auto simp add: pow2_def) lemma pow2_int: "pow2 (int c) = 2^c" by (simp add: pow2_def) -lemma float_transfer_nat: "float (a, b) = float (a * 2^c, b - int c)" -by (simp add: float_def pow2_int[symmetric] pow2_add[symmetric]) - -lemma real_is_int_real[simp]: "real_is_int (real (x::int))" -by (auto simp add: real_is_int_def int_of_real_def) - -lemma int_of_real_real[simp]: "int_of_real (real x) = x" -by (simp add: int_of_real_def) - -lemma real_int_of_real[simp]: "real_is_int x \ real (int_of_real x) = x" -by (auto simp add: int_of_real_def real_is_int_def) - -lemma real_is_int_add_int_of_real: "real_is_int a \ real_is_int b \ (int_of_real (a+b)) = (int_of_real a) + (int_of_real b)" -by (auto simp add: int_of_real_def real_is_int_def) - -lemma real_is_int_add[simp]: "real_is_int a \ real_is_int b \ real_is_int (a+b)" -apply (subst real_is_int_def2) -apply (simp add: real_is_int_add_int_of_real real_int_of_real) -done - -lemma int_of_real_sub: "real_is_int a \ real_is_int b \ (int_of_real (a-b)) = (int_of_real a) - (int_of_real b)" -by (auto simp add: int_of_real_def real_is_int_def) - -lemma real_is_int_sub[simp]: "real_is_int a \ real_is_int b \ real_is_int (a-b)" -apply (subst real_is_int_def2) -apply (simp add: int_of_real_sub real_int_of_real) -done - -lemma real_is_int_rep: "real_is_int x \ ?! (a::int). real a = x" -by (auto simp add: real_is_int_def) - -lemma int_of_real_mult: - assumes "real_is_int a" "real_is_int b" - shows "(int_of_real (a*b)) = (int_of_real a) * (int_of_real b)" -proof - - from prems have a: "?! (a'::int). real a' = a" by (rule_tac real_is_int_rep, auto) - from prems have b: "?! (b'::int). real b' = b" by (rule_tac real_is_int_rep, auto) - from a obtain a'::int where a':"a = real a'" by auto - from b obtain b'::int where b':"b = real b'" by auto - have r: "real a' * real b' = real (a' * b')" by auto - show ?thesis - apply (simp add: a' b') - apply (subst r) - apply (simp only: int_of_real_real) - done -qed - -lemma real_is_int_mult[simp]: "real_is_int a \ real_is_int b \ real_is_int (a*b)" -apply (subst real_is_int_def2) -apply (simp add: int_of_real_mult) -done - -lemma real_is_int_0[simp]: "real_is_int (0::real)" -by (simp add: real_is_int_def int_of_real_def) - -lemma real_is_int_1[simp]: "real_is_int (1::real)" -proof - - have "real_is_int (1::real) = real_is_int(real (1::int))" by auto - also have "\ = True" by (simp only: real_is_int_real) - ultimately show ?thesis by auto -qed - -lemma real_is_int_n1: "real_is_int (-1::real)" -proof - - have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto - also have "\ = True" by (simp only: real_is_int_real) - ultimately show ?thesis by auto -qed - -lemma real_is_int_number_of[simp]: "real_is_int ((number_of \ int \ real) x)" -proof - - have neg1: "real_is_int (-1::real)" - proof - - have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto - also have "\ = True" by (simp only: real_is_int_real) - ultimately show ?thesis by auto - qed - - { - fix x :: int - have "real_is_int ((number_of \ int \ real) x)" - unfolding number_of_eq - apply (induct x) - apply (induct_tac n) - apply (simp) - apply (simp) - apply (induct_tac n) - apply (simp add: neg1) - proof - - fix n :: nat - assume rn: "(real_is_int (of_int (- (int (Suc n)))))" - have s: "-(int (Suc (Suc n))) = -1 + - (int (Suc n))" by simp - show "real_is_int (of_int (- (int (Suc (Suc n)))))" - apply (simp only: s of_int_add) - apply (rule real_is_int_add) - apply (simp add: neg1) - apply (simp only: rn) - done - qed - } - note Abs_Bin = this - { - fix x :: int - have "? u. x = u" - apply (rule exI[where x = "x"]) - apply (simp) - done - } - then obtain u::int where "x = u" by auto - with Abs_Bin show ?thesis by auto -qed - -lemma int_of_real_0[simp]: "int_of_real (0::real) = (0::int)" -by (simp add: int_of_real_def) - -lemma int_of_real_1[simp]: "int_of_real (1::real) = (1::int)" -proof - - have 1: "(1::real) = real (1::int)" by auto - show ?thesis by (simp only: 1 int_of_real_real) -qed - -lemma int_of_real_number_of[simp]: "int_of_real (number_of b) = number_of b" -proof - - have "real_is_int (number_of b)" by simp - then have uu: "?! u::int. number_of b = real u" by (auto simp add: real_is_int_rep) - then obtain u::int where u:"number_of b = real u" by auto - have "number_of b = real ((number_of b)::int)" - by (simp add: number_of_eq real_of_int_def) - have ub: "number_of b = real ((number_of b)::int)" - by (simp add: number_of_eq real_of_int_def) - from uu u ub have unb: "u = number_of b" - by blast - have "int_of_real (number_of b) = u" by (simp add: u) - with unb show ?thesis by simp -qed - -lemma float_transfer_even: "even a \ float (a, b) = float (a div 2, b+1)" - apply (subst float_transfer[where a="a" and b="b" and c="-1", simplified]) - apply (simp_all add: pow2_def even_def real_is_int_def algebra_simps) - apply (auto) -proof - - fix q::int - have a:"b - (-1\int) = (1\int) + b" by arith - show "(float (q, (b - (-1\int)))) = (float (q, ((1\int) + b)))" - by (simp add: a) -qed - -lemma int_div_zdiv: "int (a div b) = (int a) div (int b)" -by (rule zdiv_int) - -lemma int_mod_zmod: "int (a mod b) = (int a) mod (int b)" -by (rule zmod_int) - -lemma abs_div_2_less: "a \ 0 \ a \ -1 \ abs((a::int) div 2) < abs a" -by arith - -function norm_float :: "int \ int \ int \ int" where - "norm_float a b = (if a \ 0 \ even a then norm_float (a div 2) (b + 1) - else if a = 0 then (0, 0) else (a, b))" -by auto - -termination by (relation "measure (nat o abs o fst)") - (auto intro: abs_div_2_less) - -lemma norm_float: "float x = float (split norm_float x)" -proof - - { - fix a b :: int - have norm_float_pair: "float (a, b) = float (norm_float a b)" - proof (induct a b rule: norm_float.induct) - case (1 u v) - show ?case - proof cases - assume u: "u \ 0 \ even u" - with prems have ind: "float (u div 2, v + 1) = float (norm_float (u div 2) (v + 1))" by auto - with u have "float (u,v) = float (u div 2, v+1)" by (simp add: float_transfer_even) - then show ?thesis - apply (subst norm_float.simps) - apply (simp add: ind) - done - next - assume "~(u \ 0 \ even u)" - then show ?thesis - by (simp add: prems float_def) - qed - qed - } - note helper = this - have "? a b. x = (a,b)" by auto - then obtain a b where "x = (a, b)" by blast - then show ?thesis by (simp add: helper) -qed - -lemma float_add_l0: "float (0, e) + x = x" - by (simp add: float_def) - -lemma float_add_r0: "x + float (0, e) = x" - by (simp add: float_def) - -lemma float_add: - "float (a1, e1) + float (a2, e2) = - (if e1<=e2 then float (a1+a2*2^(nat(e2-e1)), e1) - else float (a1*2^(nat (e1-e2))+a2, e2))" - apply (simp add: float_def algebra_simps) - apply (auto simp add: pow2_int[symmetric] pow2_add[symmetric]) - done - -lemma float_add_assoc1: - "(x + float (y1, e1)) + float (y2, e2) = (float (y1, e1) + float (y2, e2)) + x" - by simp - -lemma float_add_assoc2: - "(float (y1, e1) + x) + float (y2, e2) = (float (y1, e1) + float (y2, e2)) + x" - by simp - -lemma float_add_assoc3: - "float (y1, e1) + (x + float (y2, e2)) = (float (y1, e1) + float (y2, e2)) + x" - by simp - -lemma float_add_assoc4: - "float (y1, e1) + (float (y2, e2) + x) = (float (y1, e1) + float (y2, e2)) + x" - by simp - -lemma float_mult_l0: "float (0, e) * x = float (0, 0)" - by (simp add: float_def) - -lemma float_mult_r0: "x * float (0, e) = float (0, 0)" - by (simp add: float_def) - -definition - lbound :: "real \ real" -where - "lbound x = min 0 x" - -definition - ubound :: "real \ real" -where - "ubound x = max 0 x" - -lemma lbound: "lbound x \ x" - by (simp add: lbound_def) - -lemma ubound: "x \ ubound x" - by (simp add: ubound_def) - -lemma float_mult: - "float (a1, e1) * float (a2, e2) = - (float (a1 * a2, e1 + e2))" - by (simp add: float_def pow2_add) - -lemma float_minus: - "- (float (a,b)) = float (-a, b)" - by (simp add: float_def) - -lemma zero_less_pow2: +lemma zero_less_pow2[simp]: "0 < pow2 x" proof - { @@ -387,182 +175,1258 @@ done qed +lemma normfloat_imp_odd_or_zero: "normfloat f = Float a b \ odd a \ (a = 0 \ b = 0)" +proof (induct f rule: normfloat.induct) + case (1 u v) + from 1 have ab: "normfloat (Float u v) = Float a b" by auto + { + assume eu: "even u" + assume z: "u \ 0" + have "normfloat (Float u v) = normfloat (Float (u div 2) (v + 1))" + apply (subst normfloat.simps) + by (simp add: eu z) + with ab have "normfloat (Float (u div 2) (v + 1)) = Float a b" by simp + with 1 eu z have ?case by auto + } + note case1 = this + { + assume "odd u \ u = 0" + then have ou: "\ (u \ 0 \ even u)" by auto + have "normfloat (Float u v) = (if u = 0 then Float 0 0 else Float u v)" + apply (subst normfloat.simps) + apply (simp add: ou) + done + with ab have "Float a b = (if u = 0 then Float 0 0 else Float u v)" by auto + then have ?case + apply (case_tac "u=0") + apply (auto) + by (insert ou, auto) + } + note case2 = this + show ?case + apply (case_tac "odd u \ u = 0") + apply (rule case2) + apply simp + apply (rule case1) + apply auto + done +qed + +lemma float_eq_odd_helper: + assumes odd: "odd a'" + and floateq: "Ifloat (Float a b) = Ifloat (Float a' b')" + shows "b \ b'" +proof - + { + assume bcmp: "b > b'" + from floateq have eq: "real a * pow2 b = real a' * pow2 b'" by simp + { + fix x y z :: real + assume "y \ 0" + then have "(x * inverse y = z) = (x = z * y)" + by auto + } + note inverse = this + have eq': "real a * (pow2 (b - b')) = real a'" + apply (subst diff_int_def) + apply (subst pow2_add) + apply (subst pow2_neg[where x = "-b'"]) + apply simp + apply (subst mult_assoc[symmetric]) + apply (subst inverse) + apply (simp_all add: eq) + done + have "\ z > 0. pow2 (b-b') = 2^z" + apply (rule exI[where x="nat (b - b')"]) + apply (auto) + apply (insert bcmp) + apply simp + apply (subst pow2_int[symmetric]) + apply auto + done + then obtain z where z: "z > 0 \ pow2 (b-b') = 2^z" by auto + with eq' have "real a * 2^z = real a'" + by auto + then have "real a * real ((2::int)^z) = real a'" + by auto + then have "real (a * 2^z) = real a'" + apply (subst real_of_int_mult) + apply simp + done + then have a'_rep: "a * 2^z = a'" by arith + then have "a' = a*2^z" by simp + with z have "even a'" by simp + with odd have False by auto + } + then show ?thesis by arith +qed + +lemma float_eq_odd: + assumes odd1: "odd a" + and odd2: "odd a'" + and floateq: "Ifloat (Float a b) = Ifloat (Float a' b')" + shows "a = a' \ b = b'" +proof - + from + float_eq_odd_helper[OF odd2 floateq] + float_eq_odd_helper[OF odd1 floateq[symmetric]] + have beq: "b = b'" by arith + with floateq show ?thesis by auto +qed + +theorem normfloat_unique: + assumes Ifloat_eq: "Ifloat f = Ifloat g" + shows "normfloat f = normfloat g" +proof - + from float_split[of "normfloat f"] obtain a b where normf:"normfloat f = Float a b" by auto + from float_split[of "normfloat g"] obtain a' b' where normg:"normfloat g = Float a' b'" by auto + have "Ifloat (normfloat f) = Ifloat (normfloat g)" + by (simp add: Ifloat_eq) + then have float_eq: "Ifloat (Float a b) = Ifloat (Float a' b')" + by (simp add: normf normg) + have ab: "odd a \ (a = 0 \ b = 0)" by (rule normfloat_imp_odd_or_zero[OF normf]) + have ab': "odd a' \ (a' = 0 \ b' = 0)" by (rule normfloat_imp_odd_or_zero[OF normg]) + { + assume odd: "odd a" + then have "a \ 0" by (simp add: even_def, arith) + with float_eq have "a' \ 0" by auto + with ab' have "odd a'" by simp + from odd this float_eq have "a = a' \ b = b'" by (rule float_eq_odd) + } + note odd_case = this + { + assume even: "even a" + with ab have a0: "a = 0" by simp + with float_eq have a0': "a' = 0" by auto + from a0 a0' ab ab' have "a = a' \ b = b'" by auto + } + note even_case = this + from odd_case even_case show ?thesis + apply (simp add: normf normg) + apply (case_tac "even a") + apply auto + done +qed + +instantiation float :: plus begin +fun plus_float where +[simp del]: "(Float a_m a_e) + (Float b_m b_e) = + (if a_e \ b_e then Float (a_m + b_m * 2^(nat(b_e - a_e))) a_e + else Float (a_m * 2^(nat (a_e - b_e)) + b_m) b_e)" +instance .. +end + +instantiation float :: uminus begin +fun uminus_float where [simp del]: "uminus_float (Float m e) = Float (-m) e" +instance .. +end + +instantiation float :: minus begin +fun minus_float where [simp del]: "(z::float) - w = z + (- w)" +instance .. +end + +instantiation float :: times begin +fun times_float where [simp del]: "(Float a_m a_e) * (Float b_m b_e) = Float (a_m * b_m) (a_e + b_e)" +instance .. +end + +fun float_pprt :: "float \ float" where +"float_pprt (Float a e) = (if 0 <= a then (Float a e) else 0)" + +fun float_nprt :: "float \ float" where +"float_nprt (Float a e) = (if 0 <= a then 0 else (Float a e))" + +instantiation float :: ord begin +definition le_float_def: "z \ w \ Ifloat z \ Ifloat w" +definition less_float_def: "z < w \ Ifloat z < Ifloat w" +instance .. +end + +lemma Ifloat_add[simp]: "Ifloat (a + b) = Ifloat a + Ifloat b" + by (cases a, cases b, simp add: algebra_simps plus_float.simps, + auto simp add: pow2_int[symmetric] pow2_add[symmetric]) + +lemma Ifloat_minus[simp]: "Ifloat (- a) = - Ifloat a" + by (cases a, simp add: uminus_float.simps) + +lemma Ifloat_sub[simp]: "Ifloat (a - b) = Ifloat a - Ifloat b" + by (cases a, cases b, simp add: minus_float.simps) + +lemma Ifloat_mult[simp]: "Ifloat (a*b) = Ifloat a * Ifloat b" + by (cases a, cases b, simp add: times_float.simps pow2_add) + +lemma Ifloat_0[simp]: "Ifloat 0 = 0" + by (auto simp add: zero_float_def float_zero) + +lemma Ifloat_1[simp]: "Ifloat 1 = 1" + by (auto simp add: one_float_def) + lemma zero_le_float: - "(0 <= float (a,b)) = (0 <= a)" - apply (auto simp add: float_def) - apply (auto simp add: zero_le_mult_iff zero_less_pow2) + "(0 <= Ifloat (Float a b)) = (0 <= a)" + apply auto + apply (auto simp add: zero_le_mult_iff) apply (insert zero_less_pow2[of b]) apply (simp_all) done lemma float_le_zero: - "(float (a,b) <= 0) = (a <= 0)" - apply (auto simp add: float_def) + "(Ifloat (Float a b) <= 0) = (a <= 0)" + apply auto apply (auto simp add: mult_le_0_iff) apply (insert zero_less_pow2[of b]) apply auto done -lemma float_abs: - "abs (float (a,b)) = (if 0 <= a then (float (a,b)) else (float (-a,b)))" - apply (auto simp add: abs_if) - apply (simp_all add: zero_le_float[symmetric, of a b] float_minus) +declare Ifloat.simps[simp del] + +lemma Ifloat_pprt[simp]: "Ifloat (float_pprt a) = pprt (Ifloat a)" + by (cases a, auto simp add: float_pprt.simps zero_le_float float_le_zero float_zero) + +lemma Ifloat_nprt[simp]: "Ifloat (float_nprt a) = nprt (Ifloat a)" + by (cases a, auto simp add: float_nprt.simps zero_le_float float_le_zero float_zero) + +instance float :: ab_semigroup_add +proof (intro_classes) + fix a b c :: float + show "a + b + c = a + (b + c)" + by (cases a, cases b, cases c, auto simp add: algebra_simps power_add[symmetric] plus_float.simps) +next + fix a b :: float + show "a + b = b + a" + by (cases a, cases b, simp add: plus_float.simps) +qed + +instance float :: comm_monoid_mult +proof (intro_classes) + fix a b c :: float + show "a * b * c = a * (b * c)" + by (cases a, cases b, cases c, simp add: times_float.simps) +next + fix a b :: float + show "a * b = b * a" + by (cases a, cases b, simp add: times_float.simps) +next + fix a :: float + show "1 * a = a" + by (cases a, simp add: times_float.simps one_float_def) +qed + +(* Floats do NOT form a cancel_semigroup_add: *) +lemma "0 + Float 0 1 = 0 + Float 0 2" + by (simp add: plus_float.simps zero_float_def) + +instance float :: comm_semiring +proof (intro_classes) + fix a b c :: float + show "(a + b) * c = a * c + b * c" + by (cases a, cases b, cases c, simp, simp add: plus_float.simps times_float.simps algebra_simps) +qed + +(* Floats do NOT form an order, because "(x < y) = (x <= y & x <> y)" does NOT hold *) + +instance float :: zero_neq_one +proof (intro_classes) + show "(0::float) \ 1" + by (simp add: zero_float_def one_float_def) +qed + +lemma float_le_simp: "((x::float) \ y) = (0 \ y - x)" + by (auto simp add: le_float_def) + +lemma float_less_simp: "((x::float) < y) = (0 < y - x)" + by (auto simp add: less_float_def) + +lemma Ifloat_min: "Ifloat (min x y) = min (Ifloat x) (Ifloat y)" unfolding min_def le_float_def by auto +lemma Ifloat_max: "Ifloat (max a b) = max (Ifloat a) (Ifloat b)" unfolding max_def le_float_def by auto + +instantiation float :: power begin +fun power_float where [simp del]: "(Float m e) ^ n = Float (m ^ n) (e * int n)" +instance .. +end + +instance float :: recpower +proof (intro_classes) + fix a :: float show "a ^ 0 = 1" by (cases a, auto simp add: power_float.simps one_float_def) +next + fix a :: float and n :: nat show "a ^ (Suc n) = a * a ^ n" + by (cases a, auto simp add: power_float.simps times_float.simps algebra_simps) +qed + +lemma float_power: "Ifloat (x ^ n) = (Ifloat x) ^ n" +proof (cases x) + case (Float m e) + + have "pow2 e ^ n = pow2 (e * int n)" + proof (cases "e >= 0") + case True hence e_nat: "e = int (nat e)" by auto + hence "pow2 e ^ n = (2 ^ nat e) ^ n" using pow2_int[of "nat e"] by auto + thus ?thesis unfolding power_mult[symmetric] unfolding pow2_int[symmetric] int_mult e_nat[symmetric] . + next + case False hence e_minus: "-e = int (nat (-e))" by auto + hence "pow2 (-e) ^ n = (2 ^ nat (-e)) ^ n" using pow2_int[of "nat (-e)"] by auto + hence "pow2 (-e) ^ n = pow2 ((-e) * int n)" unfolding power_mult[symmetric] unfolding pow2_int[symmetric] int_mult e_minus[symmetric] zmult_zminus . + thus ?thesis unfolding pow2_neg[of "-e"] pow2_neg[of "-e * int n"] unfolding zmult_zminus zminus_zminus nonzero_power_inverse[OF pow2_neq_zero, symmetric] + using nonzero_inverse_eq_imp_eq[OF _ pow2_neq_zero pow2_neq_zero] by auto + qed + thus ?thesis by (auto simp add: Float power_mult_distrib Ifloat.simps power_float.simps) +qed + +lemma zero_le_pow2[simp]: "0 \ pow2 s" + apply (subgoal_tac "0 < pow2 s") + apply (auto simp only:) + apply auto done -lemma float_zero: - "float (0, b) = 0" - by (simp add: float_def) - -lemma float_pprt: - "pprt (float (a, b)) = (if 0 <= a then (float (a,b)) else (float (0, b)))" - by (auto simp add: zero_le_float float_le_zero float_zero) - -lemma pprt_lbound: "pprt (lbound x) = float (0, 0)" - apply (simp add: float_def) - apply (rule pprt_eq_0) - apply (simp add: lbound_def) +lemma pow2_less_0_eq_False[simp]: "(pow2 s < 0) = False" + apply auto + apply (subgoal_tac "0 \ pow2 s") + apply simp + apply simp done -lemma nprt_ubound: "nprt (ubound x) = float (0, 0)" - apply (simp add: float_def) - apply (rule nprt_eq_0) - apply (simp add: ubound_def) +lemma pow2_le_0_eq_False[simp]: "(pow2 s \ 0) = False" + apply auto + apply (subgoal_tac "0 < pow2 s") + apply simp + apply simp done -lemma float_nprt: - "nprt (float (a, b)) = (if 0 <= a then (float (0,b)) else (float (a, b)))" - by (auto simp add: zero_le_float float_le_zero float_zero) - -lemma norm_0_1: "(0::_::number_ring) = Numeral0 & (1::_::number_ring) = Numeral1" +lemma float_pos_m_pos: "0 < Float m e \ 0 < m" + unfolding less_float_def Ifloat.simps Ifloat_0 zero_less_mult_iff by auto -lemma add_left_zero: "0 + a = (a::'a::comm_monoid_add)" - by simp +lemma float_pos_less1_e_neg: assumes "0 < Float m e" and "Float m e < 1" shows "e < 0" +proof - + have "0 < m" using float_pos_m_pos `0 < Float m e` by auto + hence "0 \ real m" and "1 \ real m" by auto + + show "e < 0" + proof (rule ccontr) + assume "\ e < 0" hence "0 \ e" by auto + hence "1 \ pow2 e" unfolding pow2_def by auto + from mult_mono[OF `1 \ real m` this `0 \ real m`] + have "1 \ Float m e" by (simp add: le_float_def Ifloat.simps) + thus False using `Float m e < 1` unfolding less_float_def le_float_def by auto + qed +qed + +lemma float_less1_mantissa_bound: assumes "0 < Float m e" "Float m e < 1" shows "m < 2^(nat (-e))" +proof - + have "e < 0" using float_pos_less1_e_neg assms by auto + have "\x. (0::real) < 2^x" by auto + have "real m < 2^(nat (-e))" using `Float m e < 1` + unfolding less_float_def Ifloat_neg_exp[OF `e < 0`] Ifloat_1 + real_mult_less_iff1[of _ _ 1, OF `0 < 2^(nat (-e))`, symmetric] + real_mult_assoc by auto + thus ?thesis unfolding real_of_int_less_iff[symmetric] by auto +qed + +function bitlen :: "int \ int" where +"bitlen 0 = 0" | +"bitlen -1 = 1" | +"0 < x \ bitlen x = 1 + (bitlen (x div 2))" | +"x < -1 \ bitlen x = 1 + (bitlen (x div 2))" + apply (case_tac "x = 0 \ x = -1 \ x < -1 \ x > 0") + apply auto + done +termination by (relation "measure (nat o abs)", auto) + +lemma bitlen_ge0: "0 \ bitlen x" by (induct x rule: bitlen.induct, auto) +lemma bitlen_ge1: "x \ 0 \ 1 \ bitlen x" by (induct x rule: bitlen.induct, auto simp add: bitlen_ge0) + +lemma bitlen_bounds': assumes "0 < x" shows "2^nat (bitlen x - 1) \ x \ x + 1 \ 2^nat (bitlen x)" (is "?P x") + using `0 < x` +proof (induct x rule: bitlen.induct) + fix x + assume "0 < x" and hyp: "0 < x div 2 \ ?P (x div 2)" hence "0 \ x" and "x \ 0" by auto + { fix x have "0 \ 1 + bitlen x" using bitlen_ge0[of x] by auto } note gt0_pls1 = this + + have "0 < (2::int)" by auto -lemma add_right_zero: "a + 0 = (a::'a::comm_monoid_add)" - by simp + show "?P x" + proof (cases "x = 1") + case True show "?P x" unfolding True by auto + next + case False hence "2 \ x" using `0 < x` `x \ 1` by auto + hence "2 div 2 \ x div 2" by (rule zdiv_mono1, auto) + hence "0 < x div 2" and "x div 2 \ 0" by auto + hence bitlen_s1_ge0: "0 \ bitlen (x div 2) - 1" using bitlen_ge1[OF `x div 2 \ 0`] by auto -lemma mult_left_one: "1 * a = (a::'a::semiring_1)" - by simp + { from hyp[OF `0 < x div 2`] + have "2 ^ nat (bitlen (x div 2) - 1) \ x div 2" by auto + hence "2 ^ nat (bitlen (x div 2) - 1) * 2 \ x div 2 * 2" by (rule mult_right_mono, auto) + also have "\ \ x" using `0 < x` by auto + finally have "2^nat (1 + bitlen (x div 2) - 1) \ x" unfolding power_Suc2[symmetric] Suc_nat_eq_nat_zadd1[OF bitlen_s1_ge0] by auto + } moreover + { have "x + 1 \ x - x mod 2 + 2" + proof - + have "x mod 2 < 2" using `0 < x` by auto + hence "x < x - x mod 2 + 2" unfolding algebra_simps by auto + thus ?thesis by auto + qed + also have "x - x mod 2 + 2 = (x div 2 + 1) * 2" unfolding algebra_simps using `0 < x` zdiv_zmod_equality2[of x 2 0] by auto + also have "\ \ 2^nat (bitlen (x div 2)) * 2" using hyp[OF `0 < x div 2`, THEN conjunct2] by (rule mult_right_mono, auto) + also have "\ = 2^(1 + nat (bitlen (x div 2)))" unfolding power_Suc2[symmetric] by auto + finally have "x + 1 \ 2^(1 + nat (bitlen (x div 2)))" . + } + ultimately show ?thesis + unfolding bitlen.simps(3)[OF `0 < x`] nat_add_distrib[OF zero_le_one bitlen_ge0] + unfolding add_commute nat_add_distrib[OF zero_le_one gt0_pls1] + by auto + qed +next + fix x :: int assume "x < -1" and "0 < x" hence False by auto + thus "?P x" by auto +qed auto + +lemma bitlen_bounds: assumes "0 < x" shows "2^nat (bitlen x - 1) \ x \ x < 2^nat (bitlen x)" + using bitlen_bounds'[OF `0 real m / 2^nat (bitlen m - 1)" and "real m / 2^nat (bitlen m - 1) < 2" +proof - + let ?B = "2^nat(bitlen m - 1)" + + have "?B \ m" using bitlen_bounds[OF `0 real m" unfolding real_of_int_le_iff[symmetric] by auto + thus "1 \ real m / ?B" by auto + + have "m \ 0" using assms by auto + have "0 \ bitlen m - 1" using bitlen_ge1[OF `m \ 0`] by auto -lemma mult_right_one: "a * 1 = (a::'a::semiring_1)" - by simp + have "m < 2^nat(bitlen m)" using bitlen_bounds[OF `0 = 2^nat(bitlen m - 1 + 1)" using bitlen_ge1[OF `m \ 0`] by auto + also have "\ = ?B * 2" unfolding nat_add_distrib[OF `0 \ bitlen m - 1` zero_le_one] by auto + finally have "real m < 2 * ?B" unfolding real_of_int_less_iff[symmetric] by auto + hence "real m / ?B < 2 * ?B / ?B" by (rule divide_strict_right_mono, auto) + thus "real m / ?B < 2" by auto +qed + +lemma float_gt1_scale: assumes "1 \ Float m e" + shows "0 \ e + (bitlen m - 1)" +proof (cases "0 \ e") + have "0 < Float m e" using assms unfolding less_float_def le_float_def by auto + hence "0 < m" using float_pos_m_pos by auto + hence "m \ 0" by auto + case True with bitlen_ge1[OF `m \ 0`] show ?thesis by auto +next + have "0 < Float m e" using assms unfolding less_float_def le_float_def by auto + hence "0 < m" using float_pos_m_pos by auto + hence "m \ 0" and "1 < (2::int)" by auto + case False let ?S = "2^(nat (-e))" + have "1 \ real m * inverse ?S" using assms unfolding le_float_def Ifloat_nge0_exp[OF False] by auto + hence "1 * ?S \ real m * inverse ?S * ?S" by (rule mult_right_mono, auto) + hence "?S \ real m" unfolding mult_assoc by auto + hence "?S \ m" unfolding real_of_int_le_iff[symmetric] by auto + from this bitlen_bounds[OF `0 < m`, THEN conjunct2] + have "nat (-e) < (nat (bitlen m))" unfolding power_strict_increasing_iff[OF `1 < 2`, symmetric] by (rule order_le_less_trans) + hence "-e < bitlen m" using False bitlen_ge0 by auto + thus ?thesis by auto +qed + +lemma normalized_float: assumes "m \ 0" shows "Ifloat (Float m (- (bitlen m - 1))) = real m / 2^nat (bitlen m - 1)" +proof (cases "- (bitlen m - 1) = 0") + case True show ?thesis unfolding Ifloat.simps pow2_def using True by auto +next + case False hence P: "\ 0 \ - (bitlen m - 1)" using bitlen_ge1[OF `m \ 0`] by auto + show ?thesis unfolding Ifloat_nge0_exp[OF P] real_divide_def by auto +qed + +lemma bitlen_Pls: "bitlen (Int.Pls) = Int.Pls" by (subst Pls_def, subst Pls_def, simp) + +lemma bitlen_Min: "bitlen (Int.Min) = Int.Bit1 Int.Pls" by (subst Min_def, simp add: Bit1_def) + +lemma bitlen_B0: "bitlen (Int.Bit0 b) = (if iszero b then Int.Pls else Int.succ (bitlen b))" + apply (auto simp add: iszero_def succ_def) + apply (simp add: Bit0_def Pls_def) + apply (subst Bit0_def) + apply simp + apply (subgoal_tac "0 < 2 * b \ 2 * b < -1") + apply auto + done -lemma int_pow_0: "(a::int)^(Numeral0) = 1" - by simp +lemma bitlen_B1: "bitlen (Int.Bit1 b) = (if iszero (Int.succ b) then Int.Bit1 Int.Pls else Int.succ (bitlen b))" +proof - + have h: "! x. (2*x + 1) div 2 = (x::int)" + by arith + show ?thesis + apply (auto simp add: iszero_def succ_def) + apply (subst Bit1_def)+ + apply simp + apply (subgoal_tac "2 * b + 1 = -1") + apply (simp only:) + apply simp_all + apply (subst Bit1_def) + apply simp + apply (subgoal_tac "0 < 2 * b + 1 \ 2 * b + 1 < -1") + apply (auto simp add: h) + done +qed + +lemma bitlen_number_of: "bitlen (number_of w) = number_of (bitlen w)" + by (simp add: number_of_is_id) -lemma int_pow_1: "(a::int)^(Numeral1) = a" - by simp +lemma [code]: "bitlen x = + (if x = 0 then 0 + else if x = -1 then 1 + else (1 + (bitlen (x div 2))))" + by (cases "x = 0 \ x = -1 \ 0 < x") auto + +definition lapprox_posrat :: "nat \ int \ int \ float" +where + "lapprox_posrat prec x y = + (let + l = nat (int prec + bitlen y - bitlen x) ; + d = (x * 2^l) div y + in normfloat (Float d (- (int l))))" + +lemma pow2_minus: "pow2 (-x) = inverse (pow2 x)" + unfolding pow2_neg[of "-x"] by auto + +lemma lapprox_posrat: + assumes x: "0 \ x" + and y: "0 < y" + shows "Ifloat (lapprox_posrat prec x y) \ real x / real y" +proof - + let ?l = "nat (int prec + bitlen y - bitlen x)" + + have "real (x * 2^?l div y) * inverse (2^?l) \ (real (x * 2^?l) / real y) * inverse (2^?l)" + by (rule mult_right_mono, fact real_of_int_div4, simp) + also have "\ \ (real x / real y) * 2^?l * inverse (2^?l)" by auto + finally have "real (x * 2^?l div y) * inverse (2^?l) \ real x / real y" unfolding real_mult_assoc by auto + thus ?thesis unfolding lapprox_posrat_def Let_def normfloat Ifloat.simps + unfolding pow2_minus pow2_int minus_minus . +qed -lemma zero_eq_Numeral0_nring: "(0::'a::number_ring) = Numeral0" - by simp +lemma real_of_int_div_mult: + fixes x y c :: int assumes "0 < y" and "0 < c" + shows "real (x div y) \ real (x * c div y) * inverse (real c)" +proof - + have "c * (x div y) + 0 \ c * x div y" unfolding zdiv_zmult1_eq[of c x y] + by (rule zadd_left_mono, + auto intro!: mult_nonneg_nonneg + simp add: pos_imp_zdiv_nonneg_iff[OF `0 < y`] `0 < c`[THEN less_imp_le] pos_mod_sign[OF `0 < y`]) + hence "real (x div y) * real c \ real (x * c div y)" + unfolding real_of_int_mult[symmetric] real_of_int_le_iff zmult_commute by auto + hence "real (x div y) * real c * inverse (real c) \ real (x * c div y) * inverse (real c)" + using `0 < c` by auto + thus ?thesis unfolding real_mult_assoc using `0 < c` by auto +qed + +lemma lapprox_posrat_bottom: assumes "0 < y" + shows "real (x div y) \ Ifloat (lapprox_posrat n x y)" +proof - + have pow: "\x. (0::int) < 2^x" by auto + show ?thesis + unfolding lapprox_posrat_def Let_def Ifloat_add normfloat Ifloat.simps pow2_minus pow2_int + using real_of_int_div_mult[OF `0 < y` pow] by auto +qed + +lemma lapprox_posrat_nonneg: assumes "0 \ x" and "0 < y" + shows "0 \ Ifloat (lapprox_posrat n x y)" +proof - + show ?thesis + unfolding lapprox_posrat_def Let_def Ifloat_add normfloat Ifloat.simps pow2_minus pow2_int + using pos_imp_zdiv_nonneg_iff[OF `0 < y`] assms by (auto intro!: mult_nonneg_nonneg) +qed + +definition rapprox_posrat :: "nat \ int \ int \ float" +where + "rapprox_posrat prec x y = (let + l = nat (int prec + bitlen y - bitlen x) ; + X = x * 2^l ; + d = X div y ; + m = X mod y + in normfloat (Float (d + (if m = 0 then 0 else 1)) (- (int l))))" -lemma one_eq_Numeral1_nring: "(1::'a::number_ring) = Numeral1" - by simp +lemma rapprox_posrat: + assumes x: "0 \ x" + and y: "0 < y" + shows "real x / real y \ Ifloat (rapprox_posrat prec x y)" +proof - + let ?l = "nat (int prec + bitlen y - bitlen x)" let ?X = "x * 2^?l" + show ?thesis + proof (cases "?X mod y = 0") + case True hence "y \ 0" and "y dvd ?X" using `0 < y` by auto + from real_of_int_div[OF this] + have "real (?X div y) * inverse (2 ^ ?l) = real ?X / real y * inverse (2 ^ ?l)" by auto + also have "\ = real x / real y * (2^?l * inverse (2^?l))" by auto + finally have "real (?X div y) * inverse (2^?l) = real x / real y" by auto + thus ?thesis unfolding rapprox_posrat_def Let_def normfloat if_P[OF True] + unfolding Ifloat.simps pow2_minus pow2_int minus_minus by auto + next + case False + have "0 \ real y" and "real y \ 0" using `0 < y` by auto + have "0 \ real y * 2^?l" by (rule mult_nonneg_nonneg, rule `0 \ real y`, auto) -lemma zero_eq_Numeral0_nat: "(0::nat) = Numeral0" - by simp + have "?X = y * (?X div y) + ?X mod y" by auto + also have "\ \ y * (?X div y) + y" by (rule add_mono, auto simp add: pos_mod_bound[OF `0 < y`, THEN less_imp_le]) + also have "\ = y * (?X div y + 1)" unfolding zadd_zmult_distrib2 by auto + finally have "real ?X \ real y * real (?X div y + 1)" unfolding real_of_int_le_iff real_of_int_mult[symmetric] . + hence "real ?X / (real y * 2^?l) \ real y * real (?X div y + 1) / (real y * 2^?l)" + by (rule divide_right_mono, simp only: `0 \ real y * 2^?l`) + also have "\ = real y * real (?X div y + 1) / real y / 2^?l" by auto + also have "\ = real (?X div y + 1) * inverse (2^?l)" unfolding nonzero_mult_divide_cancel_left[OF `real y \ 0`] + unfolding real_divide_def .. + finally show ?thesis unfolding rapprox_posrat_def Let_def normfloat Ifloat.simps if_not_P[OF False] + unfolding pow2_minus pow2_int minus_minus by auto + qed +qed + +lemma rapprox_posrat_le1: assumes "0 \ x" and "0 < y" and "x \ y" + shows "Ifloat (rapprox_posrat n x y) \ 1" +proof - + let ?l = "nat (int n + bitlen y - bitlen x)" let ?X = "x * 2^?l" + show ?thesis + proof (cases "?X mod y = 0") + case True hence "y \ 0" and "y dvd ?X" using `0 < y` by auto + from real_of_int_div[OF this] + have "real (?X div y) * inverse (2 ^ ?l) = real ?X / real y * inverse (2 ^ ?l)" by auto + also have "\ = real x / real y * (2^?l * inverse (2^?l))" by auto + finally have "real (?X div y) * inverse (2^?l) = real x / real y" by auto + also have "real x / real y \ 1" using `0 \ x` and `0 < y` and `x \ y` by auto + finally show ?thesis unfolding rapprox_posrat_def Let_def normfloat if_P[OF True] + unfolding Ifloat.simps pow2_minus pow2_int minus_minus by auto + next + case False + have "x \ y" + proof (rule ccontr) + assume "\ x \ y" hence "x = y" by auto + have "?X mod y = 0" unfolding `x = y` using zmod_zmult_self2 by auto + thus False using False by auto + qed + hence "x < y" using `x \ y` by auto + hence "real x / real y < 1" using `0 < y` and `0 \ x` by auto -lemma one_eq_Numeral1_nat: "(1::nat) = Numeral1" - by simp + from real_of_int_div4[of "?X" y] + have "real (?X div y) \ (real x / real y) * 2^?l" unfolding real_of_int_mult times_divide_eq_left real_of_int_power[symmetric] real_number_of . + also have "\ < 1 * 2^?l" using `real x / real y < 1` by (rule mult_strict_right_mono, auto) + finally have "?X div y < 2^?l" unfolding real_of_int_less_iff[of _ "2^?l", symmetric] by auto + hence "?X div y + 1 \ 2^?l" by auto + hence "real (?X div y + 1) * inverse (2^?l) \ 2^?l * inverse (2^?l)" + unfolding real_of_int_le_iff[of _ "2^?l", symmetric] real_of_int_power[symmetric] real_number_of + by (rule mult_right_mono, auto) + hence "real (?X div y + 1) * inverse (2^?l) \ 1" by auto + thus ?thesis unfolding rapprox_posrat_def Let_def normfloat Ifloat.simps if_not_P[OF False] + unfolding pow2_minus pow2_int minus_minus by auto + qed +qed -lemma zpower_Pls: "(z::int)^Numeral0 = Numeral1" - by simp +lemma zdiv_greater_zero: fixes a b :: int assumes "0 < a" and "a \ b" + shows "0 < b div a" +proof (rule ccontr) + have "0 \ b" using assms by auto + assume "\ 0 < b div a" hence "b div a = 0" using `0 \ b`[unfolded pos_imp_zdiv_nonneg_iff[OF `0 b` by auto +qed + +lemma rapprox_posrat_less1: assumes "0 \ x" and "0 < y" and "2 * x < y" and "0 < n" + shows "Ifloat (rapprox_posrat n x y) < 1" +proof (cases "x = 0") + case True thus ?thesis unfolding rapprox_posrat_def True Let_def normfloat Ifloat.simps by auto +next + case False hence "0 < x" using `0 \ x` by auto + hence "x < y" using assms by auto + + let ?l = "nat (int n + bitlen y - bitlen x)" let ?X = "x * 2^?l" + show ?thesis + proof (cases "?X mod y = 0") + case True hence "y \ 0" and "y dvd ?X" using `0 < y` by auto + from real_of_int_div[OF this] + have "real (?X div y) * inverse (2 ^ ?l) = real ?X / real y * inverse (2 ^ ?l)" by auto + also have "\ = real x / real y * (2^?l * inverse (2^?l))" by auto + finally have "real (?X div y) * inverse (2^?l) = real x / real y" by auto + also have "real x / real y < 1" using `0 \ x` and `0 < y` and `x < y` by auto + finally show ?thesis unfolding rapprox_posrat_def Let_def normfloat Ifloat.simps if_P[OF True] + unfolding pow2_minus pow2_int minus_minus by auto + next + case False + hence "(real x / real y) < 1 / 2" using `0 < y` and `0 \ x` `2 * x < y` by auto -lemma zpower_Min: "(z::int)^((-1)::nat) = Numeral1" + have "0 < ?X div y" + proof - + have "2^nat (bitlen x - 1) \ y" and "y < 2^nat (bitlen y)" + using bitlen_bounds[OF `0 < x`, THEN conjunct1] bitlen_bounds[OF `0 < y`, THEN conjunct2] `x < y` by auto + hence "(2::int)^nat (bitlen x - 1) < 2^nat (bitlen y)" by (rule order_le_less_trans) + hence "bitlen x \ bitlen y" by auto + hence len_less: "nat (bitlen x - 1) \ nat (int (n - 1) + bitlen y)" by auto + + have "x \ 0" and "y \ 0" using `0 < x` `0 < y` by auto + + have exp_eq: "nat (int (n - 1) + bitlen y) - nat (bitlen x - 1) = ?l" + using `bitlen x \ bitlen y` bitlen_ge1[OF `x \ 0`] bitlen_ge1[OF `y \ 0`] `0 < n` by auto + + have "y * 2^nat (bitlen x - 1) \ y * x" + using bitlen_bounds[OF `0 < x`, THEN conjunct1] `0 < y`[THEN less_imp_le] by (rule mult_left_mono) + also have "\ \ 2^nat (bitlen y) * x" using bitlen_bounds[OF `0 < y`, THEN conjunct2, THEN less_imp_le] `0 \ x` by (rule mult_right_mono) + also have "\ \ x * 2^nat (int (n - 1) + bitlen y)" unfolding mult_commute[of x] by (rule mult_right_mono, auto simp add: `0 \ x`) + finally have "real y * 2^nat (bitlen x - 1) * inverse (2^nat (bitlen x - 1)) \ real x * 2^nat (int (n - 1) + bitlen y) * inverse (2^nat (bitlen x - 1))" + unfolding real_of_int_le_iff[symmetric] by auto + hence "real y \ real x * (2^nat (int (n - 1) + bitlen y) / (2^nat (bitlen x - 1)))" + unfolding real_mult_assoc real_divide_def by auto + also have "\ = real x * (2^(nat (int (n - 1) + bitlen y) - nat (bitlen x - 1)))" using power_diff[of "2::real", OF _ len_less] by auto + finally have "y \ x * 2^?l" unfolding exp_eq unfolding real_of_int_le_iff[symmetric] by auto + thus ?thesis using zdiv_greater_zero[OF `0 < y`] by auto + qed + + from real_of_int_div4[of "?X" y] + have "real (?X div y) \ (real x / real y) * 2^?l" unfolding real_of_int_mult times_divide_eq_left real_of_int_power[symmetric] real_number_of . + also have "\ < 1/2 * 2^?l" using `real x / real y < 1/2` by (rule mult_strict_right_mono, auto) + finally have "?X div y * 2 < 2^?l" unfolding real_of_int_less_iff[of _ "2^?l", symmetric] by auto + hence "?X div y + 1 < 2^?l" using `0 < ?X div y` by auto + hence "real (?X div y + 1) * inverse (2^?l) < 2^?l * inverse (2^?l)" + unfolding real_of_int_less_iff[of _ "2^?l", symmetric] real_of_int_power[symmetric] real_number_of + by (rule mult_strict_right_mono, auto) + hence "real (?X div y + 1) * inverse (2^?l) < 1" by auto + thus ?thesis unfolding rapprox_posrat_def Let_def normfloat Ifloat.simps if_not_P[OF False] + unfolding pow2_minus pow2_int minus_minus by auto + qed +qed + +lemma approx_rat_pattern: fixes P and ps :: "nat * int * int" + assumes Y: "\y prec x. \y = 0; ps = (prec, x, 0)\ \ P" + and A: "\x y prec. \0 \ x; 0 < y; ps = (prec, x, y)\ \ P" + and B: "\x y prec. \x < 0; 0 < y; ps = (prec, x, y)\ \ P" + and C: "\x y prec. \x < 0; y < 0; ps = (prec, x, y)\ \ P" + and D: "\x y prec. \0 \ x; y < 0; ps = (prec, x, y)\ \ P" + shows P proof - - have 1:"((-1)::nat) = 0" - by simp - show ?thesis by (simp add: 1) + obtain prec x y where [simp]: "ps = (prec, x, y)" by (cases ps, auto) + from Y have "y = 0 \ P" by auto + moreover { assume "0 < y" have P proof (cases "0 \ x") case True with A and `0 < y` show P by auto next case False with B and `0 < y` show P by auto qed } + moreover { assume "y < 0" have P proof (cases "0 \ x") case True with D and `y < 0` show P by auto next case False with C and `y < 0` show P by auto qed } + ultimately show P by (cases "y = 0 \ 0 < y \ y < 0", auto) qed -lemma fst_cong: "a=a' \ fst (a,b) = fst (a',b)" - by simp +function lapprox_rat :: "nat \ int \ int \ float" +where + "y = 0 \ lapprox_rat prec x y = 0" +| "0 \ x \ 0 < y \ lapprox_rat prec x y = lapprox_posrat prec x y" +| "x < 0 \ 0 < y \ lapprox_rat prec x y = - (rapprox_posrat prec (-x) y)" +| "x < 0 \ y < 0 \ lapprox_rat prec x y = lapprox_posrat prec (-x) (-y)" +| "0 \ x \ y < 0 \ lapprox_rat prec x y = - (rapprox_posrat prec x (-y))" +apply simp_all by (rule approx_rat_pattern) +termination by lexicographic_order -lemma snd_cong: "b=b' \ snd (a,b) = snd (a,b')" - by simp +lemma compute_lapprox_rat[code]: + "lapprox_rat prec x y = (if y = 0 then 0 else if 0 \ x then (if 0 < y then lapprox_posrat prec x y else - (rapprox_posrat prec x (-y))) + else (if 0 < y then - (rapprox_posrat prec (-x) y) else lapprox_posrat prec (-x) (-y)))" + by auto + +lemma lapprox_rat: "Ifloat (lapprox_rat prec x y) \ real x / real y" +proof - + have h[rule_format]: "! a b b'. b' \ b \ a \ b' \ a \ (b::real)" by auto + show ?thesis + apply (case_tac "y = 0") + apply simp + apply (case_tac "0 \ x \ 0 < y") + apply (simp add: lapprox_posrat) + apply (case_tac "x < 0 \ 0 < y") + apply simp + apply (subst minus_le_iff) + apply (rule h[OF rapprox_posrat]) + apply (simp_all) + apply (case_tac "x < 0 \ y < 0") + apply simp + apply (rule h[OF _ lapprox_posrat]) + apply (simp_all) + apply (case_tac "0 \ x \ y < 0") + apply (simp) + apply (subst minus_le_iff) + apply (rule h[OF rapprox_posrat]) + apply simp_all + apply arith + done +qed -lemma lift_bool: "x \ x=True" - by simp +lemma lapprox_rat_bottom: assumes "0 \ x" and "0 < y" + shows "real (x div y) \ Ifloat (lapprox_rat n x y)" + unfolding lapprox_rat.simps(2)[OF assms] using lapprox_posrat_bottom[OF `0 int \ int \ float" +where + "y = 0 \ rapprox_rat prec x y = 0" +| "0 \ x \ 0 < y \ rapprox_rat prec x y = rapprox_posrat prec x y" +| "x < 0 \ 0 < y \ rapprox_rat prec x y = - (lapprox_posrat prec (-x) y)" +| "x < 0 \ y < 0 \ rapprox_rat prec x y = rapprox_posrat prec (-x) (-y)" +| "0 \ x \ y < 0 \ rapprox_rat prec x y = - (lapprox_posrat prec x (-y))" +apply simp_all by (rule approx_rat_pattern) +termination by lexicographic_order -lemma nlift_bool: "~x \ x=False" - by simp +lemma compute_rapprox_rat[code]: + "rapprox_rat prec x y = (if y = 0 then 0 else if 0 \ x then (if 0 < y then rapprox_posrat prec x y else - (lapprox_posrat prec x (-y))) else + (if 0 < y then - (lapprox_posrat prec (-x) y) else rapprox_posrat prec (-x) (-y)))" + by auto -lemma not_false_eq_true: "(~ False) = True" by simp +lemma rapprox_rat: "real x / real y \ Ifloat (rapprox_rat prec x y)" +proof - + have h[rule_format]: "! a b b'. b' \ b \ a \ b' \ a \ (b::real)" by auto + show ?thesis + apply (case_tac "y = 0") + apply simp + apply (case_tac "0 \ x \ 0 < y") + apply (simp add: rapprox_posrat) + apply (case_tac "x < 0 \ 0 < y") + apply simp + apply (subst le_minus_iff) + apply (rule h[OF _ lapprox_posrat]) + apply (simp_all) + apply (case_tac "x < 0 \ y < 0") + apply simp + apply (rule h[OF rapprox_posrat]) + apply (simp_all) + apply (case_tac "0 \ x \ y < 0") + apply (simp) + apply (subst le_minus_iff) + apply (rule h[OF _ lapprox_posrat]) + apply simp_all + apply arith + done +qed -lemma not_true_eq_false: "(~ True) = False" by simp +lemma rapprox_rat_le1: assumes "0 \ x" and "0 < y" and "x \ y" + shows "Ifloat (rapprox_rat n x y) \ 1" + unfolding rapprox_rat.simps(2)[OF `0 \ x` `0 < y`] using rapprox_posrat_le1[OF assms] . + +lemma rapprox_rat_neg: assumes "x < 0" and "0 < y" + shows "Ifloat (rapprox_rat n x y) \ 0" + unfolding rapprox_rat.simps(3)[OF assms] using lapprox_posrat_nonneg[of "-x" y n] assms by auto + +lemma rapprox_rat_nonneg_neg: assumes "0 \ x" and "y < 0" + shows "Ifloat (rapprox_rat n x y) \ 0" + unfolding rapprox_rat.simps(5)[OF assms] using lapprox_posrat_nonneg[of x "-y" n] assms by auto -lemmas binarith = - normalize_bin_simps - pred_bin_simps succ_bin_simps - add_bin_simps minus_bin_simps mult_bin_simps +lemma rapprox_rat_nonpos_pos: assumes "x \ 0" and "0 < y" + shows "Ifloat (rapprox_rat n x y) \ 0" +proof (cases "x = 0") + case True hence "0 \ x" by auto show ?thesis unfolding rapprox_rat.simps(2)[OF `0 \ x` `0 < y`] + unfolding True rapprox_posrat_def Let_def by auto +next + case False hence "x < 0" using assms by auto + show ?thesis using rapprox_rat_neg[OF `x < 0` `0 < y`] . +qed + +fun float_divl :: "nat \ float \ float \ float" +where + "float_divl prec (Float m1 s1) (Float m2 s2) = + (let + l = lapprox_rat prec m1 m2; + f = Float 1 (s1 - s2) + in + f * l)" -lemma int_eq_number_of_eq: - "(((number_of v)::int)=(number_of w)) = iszero ((number_of (v + uminus w))::int)" - by (rule eq_number_of_eq) +lemma float_divl: "Ifloat (float_divl prec x y) \ Ifloat x / Ifloat y" +proof - + from float_split[of x] obtain mx sx where x: "x = Float mx sx" by auto + from float_split[of y] obtain my sy where y: "y = Float my sy" by auto + have "real mx / real my \ (real mx * pow2 sx / (real my * pow2 sy)) / (pow2 (sx - sy))" + apply (case_tac "my = 0") + apply simp + apply (case_tac "my > 0") + apply (subst pos_le_divide_eq) + apply simp + apply (subst pos_le_divide_eq) + apply (simp add: mult_pos_pos) + apply simp + apply (subst pow2_add[symmetric]) + apply simp + apply (subgoal_tac "my < 0") + apply auto + apply (simp add: field_simps) + apply (subst pow2_add[symmetric]) + apply (simp add: field_simps) + done + then have "Ifloat (lapprox_rat prec mx my) \ (real mx * pow2 sx / (real my * pow2 sy)) / (pow2 (sx - sy))" + by (rule order_trans[OF lapprox_rat]) + then have "Ifloat (lapprox_rat prec mx my) * pow2 (sx - sy) \ real mx * pow2 sx / (real my * pow2 sy)" + apply (subst pos_le_divide_eq[symmetric]) + apply simp_all + done + then have "pow2 (sx - sy) * Ifloat (lapprox_rat prec mx my) \ real mx * pow2 sx / (real my * pow2 sy)" + by (simp add: algebra_simps) + then show ?thesis + by (simp add: x y Let_def Ifloat.simps) +qed -lemma int_iszero_number_of_Pls: "iszero (Numeral0::int)" - by (simp only: iszero_number_of_Pls) +lemma float_divl_lower_bound: assumes "0 \ x" and "0 < y" shows "0 \ float_divl prec x y" +proof (cases x, cases y) + fix xm xe ym ye :: int + assume x_eq: "x = Float xm xe" and y_eq: "y = Float ym ye" + have "0 \ xm" using `0 \ x`[unfolded x_eq le_float_def Ifloat.simps Ifloat_0 zero_le_mult_iff] by auto + have "0 < ym" using `0 < y`[unfolded y_eq less_float_def Ifloat.simps Ifloat_0 zero_less_mult_iff] by auto -lemma int_nonzero_number_of_Min: "~(iszero ((-1)::int))" - by simp + have "\n. 0 \ Ifloat (Float 1 n)" unfolding Ifloat.simps using zero_le_pow2 by auto + moreover have "0 \ Ifloat (lapprox_rat prec xm ym)" by (rule order_trans[OF _ lapprox_rat_bottom[OF `0 \ xm` `0 < ym`]], auto simp add: `0 \ xm` pos_imp_zdiv_nonneg_iff[OF `0 < ym`]) + ultimately show "0 \ float_divl prec x y" + unfolding x_eq y_eq float_divl.simps Let_def le_float_def Ifloat_0 by (auto intro!: mult_nonneg_nonneg) +qed + +lemma float_divl_pos_less1_bound: assumes "0 < x" and "x < 1" and "0 < prec" shows "1 \ float_divl prec 1 x" +proof (cases x) + case (Float m e) + from `0 < x` `x < 1` have "0 < m" "e < 0" using float_pos_m_pos float_pos_less1_e_neg unfolding Float by auto + let ?b = "nat (bitlen m)" and ?e = "nat (-e)" + have "1 \ m" and "m \ 0" using `0 < m` by auto + with bitlen_bounds[OF `0 < m`] have "m < 2^?b" and "(2::int) \ 2^?b" by auto + hence "1 \ bitlen m" using power_le_imp_le_exp[of "2::int" 1 ?b] by auto + hence pow_split: "nat (int prec + bitlen m - 1) = (prec - 1) + ?b" using `0 < prec` by auto + + have pow_not0: "\x. (2::real)^x \ 0" by auto -lemma int_iszero_number_of_Bit0: "iszero ((number_of (Int.Bit0 w))::int) = iszero ((number_of w)::int)" - by simp + from float_less1_mantissa_bound `0 < x` `x < 1` Float + have "m < 2^?e" by auto + with bitlen_bounds[OF `0 < m`, THEN conjunct1] + have "(2::int)^nat (bitlen m - 1) < 2^?e" by (rule order_le_less_trans) + from power_less_imp_less_exp[OF _ this] + have "bitlen m \ - e" by auto + hence "(2::real)^?b \ 2^?e" by auto + hence "(2::real)^?b * inverse (2^?b) \ 2^?e * inverse (2^?b)" by (rule mult_right_mono, auto) + hence "(1::real) \ 2^?e * inverse (2^?b)" by auto + also + let ?d = "real (2 ^ nat (int prec + bitlen m - 1) div m) * inverse (2 ^ nat (int prec + bitlen m - 1))" + { have "2^(prec - 1) * m \ 2^(prec - 1) * 2^?b" using `m < 2^?b`[THEN less_imp_le] by (rule mult_left_mono, auto) + also have "\ = 2 ^ nat (int prec + bitlen m - 1)" unfolding pow_split zpower_zadd_distrib by auto + finally have "2^(prec - 1) * m div m \ 2 ^ nat (int prec + bitlen m - 1) div m" using `0 < m` by (rule zdiv_mono1) + hence "2^(prec - 1) \ 2 ^ nat (int prec + bitlen m - 1) div m" unfolding zdiv_zmult_self1[OF `m \ 0`] . + hence "2^(prec - 1) * inverse (2 ^ nat (int prec + bitlen m - 1)) \ ?d" + unfolding real_of_int_le_iff[of "2^(prec - 1)", symmetric] by auto } + from mult_left_mono[OF this[unfolded pow_split power_add inverse_mult_distrib real_mult_assoc[symmetric] right_inverse[OF pow_not0] real_mult_1], of "2^?e"] + have "2^?e * inverse (2^?b) \ 2^?e * ?d" unfolding pow_split power_add by auto + finally have "1 \ 2^?e * ?d" . + + have e_nat: "0 - e = int (nat (-e))" using `e < 0` by auto + have "bitlen 1 = 1" using bitlen.simps by auto + + show ?thesis + unfolding one_float_def Float float_divl.simps Let_def lapprox_rat.simps(2)[OF zero_le_one `0 < m`] lapprox_posrat_def `bitlen 1 = 1` + unfolding le_float_def Ifloat_mult normfloat Ifloat.simps pow2_minus pow2_int e_nat + using `1 \ 2^?e * ?d` by (auto simp add: pow2_def) +qed -lemma int_iszero_number_of_Bit1: "\ iszero ((number_of (Int.Bit1 w))::int)" - by simp +fun float_divr :: "nat \