src/HOL/Tools/float_arith.ML
changeset 28952 15a4b2cf8c34
parent 25300 bc3a1c964704
child 29804 e15b74577368
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/HOL/Tools/float_arith.ML	Wed Dec 03 15:58:44 2008 +0100
     1.3 @@ -0,0 +1,222 @@
     1.4 +(*  Title:  HOL/Real/float_arith.ML
     1.5 +    ID:     $Id$
     1.6 +    Author: Steven Obua
     1.7 +*)
     1.8 +
     1.9 +signature FLOAT_ARITH =
    1.10 +sig
    1.11 +  exception Destruct_floatstr of string
    1.12 +  val destruct_floatstr: (char -> bool) -> (char -> bool) -> string -> bool * string * string * bool * string
    1.13 +
    1.14 +  exception Floating_point of string
    1.15 +  val approx_dec_by_bin: int -> Float.float -> Float.float * Float.float
    1.16 +  val approx_decstr_by_bin: int -> string -> Float.float * Float.float
    1.17 +
    1.18 +  val mk_float: Float.float -> term
    1.19 +  val dest_float: term -> Float.float
    1.20 +
    1.21 +  val approx_float: int -> (Float.float * Float.float -> Float.float * Float.float)
    1.22 +    -> string -> term * term
    1.23 +end;
    1.24 +
    1.25 +structure FloatArith : FLOAT_ARITH =
    1.26 +struct
    1.27 +
    1.28 +exception Destruct_floatstr of string;
    1.29 +
    1.30 +fun destruct_floatstr isDigit isExp number =
    1.31 +  let
    1.32 +    val numlist = filter (not o Char.isSpace) (String.explode number)
    1.33 +
    1.34 +    fun countsigns ((#"+")::cs) = countsigns cs
    1.35 +      | countsigns ((#"-")::cs) =
    1.36 +      let
    1.37 +        val (positive, rest) = countsigns cs
    1.38 +      in
    1.39 +        (not positive, rest)
    1.40 +      end
    1.41 +      | countsigns cs = (true, cs)
    1.42 +
    1.43 +    fun readdigits [] = ([], [])
    1.44 +      | readdigits (q as c::cs) =
    1.45 +      if (isDigit c) then
    1.46 +        let
    1.47 +          val (digits, rest) = readdigits cs
    1.48 +        in
    1.49 +          (c::digits, rest)
    1.50 +        end
    1.51 +      else
    1.52 +        ([], q)
    1.53 +
    1.54 +    fun readfromexp_helper cs =
    1.55 +      let
    1.56 +        val (positive, rest) = countsigns cs
    1.57 +        val (digits, rest') = readdigits rest
    1.58 +      in
    1.59 +        case rest' of
    1.60 +          [] => (positive, digits)
    1.61 +          | _ => raise (Destruct_floatstr number)
    1.62 +      end
    1.63 +
    1.64 +    fun readfromexp [] = (true, [])
    1.65 +      | readfromexp (c::cs) =
    1.66 +      if isExp c then
    1.67 +        readfromexp_helper cs
    1.68 +      else
    1.69 +        raise (Destruct_floatstr number)
    1.70 +
    1.71 +    fun readfromdot [] = ([], readfromexp [])
    1.72 +      | readfromdot ((#".")::cs) =
    1.73 +      let
    1.74 +        val (digits, rest) = readdigits cs
    1.75 +        val exp = readfromexp rest
    1.76 +      in
    1.77 +        (digits, exp)
    1.78 +      end
    1.79 +      | readfromdot cs = readfromdot ((#".")::cs)
    1.80 +
    1.81 +    val (positive, numlist) = countsigns numlist
    1.82 +    val (digits1, numlist) = readdigits numlist
    1.83 +     val (digits2, exp) = readfromdot numlist
    1.84 +  in
    1.85 +    (positive, String.implode digits1, String.implode digits2, fst exp, String.implode (snd exp))
    1.86 +  end
    1.87 +
    1.88 +exception Floating_point of string;
    1.89 +
    1.90 +val ln2_10 = Math.ln 10.0 / Math.ln 2.0;
    1.91 +fun exp5 x = Integer.pow x 5;
    1.92 +fun exp10 x = Integer.pow x 10;
    1.93 +fun exp2 x = Integer.pow x 2;
    1.94 +
    1.95 +fun find_most_significant q r =
    1.96 +  let
    1.97 +    fun int2real i =
    1.98 +      case (Real.fromString o string_of_int) i of
    1.99 +        SOME r => r
   1.100 +        | NONE => raise (Floating_point "int2real")
   1.101 +    fun subtract (q, r) (q', r') =
   1.102 +      if r <= r' then
   1.103 +        (q - q' * exp10 (r' - r), r)
   1.104 +      else
   1.105 +        (q * exp10 (r - r') - q', r')
   1.106 +    fun bin2dec d =
   1.107 +      if 0 <= d then
   1.108 +        (exp2 d, 0)
   1.109 +      else
   1.110 +        (exp5 (~ d), d)
   1.111 +
   1.112 +    val L = Real.floor (int2real (IntInf.log2 q) + int2real r * ln2_10)
   1.113 +    val L1 = L + 1
   1.114 +
   1.115 +    val (q1, r1) = subtract (q, r) (bin2dec L1) 
   1.116 +  in
   1.117 +    if 0 <= q1 then
   1.118 +      let
   1.119 +        val (q2, r2) = subtract (q, r) (bin2dec (L1 + 1))
   1.120 +      in
   1.121 +        if 0 <= q2 then
   1.122 +          raise (Floating_point "find_most_significant")
   1.123 +        else
   1.124 +          (L1, (q1, r1))
   1.125 +      end
   1.126 +    else
   1.127 +      let
   1.128 +        val (q0, r0) = subtract (q, r) (bin2dec L)
   1.129 +      in
   1.130 +        if 0 <= q0 then
   1.131 +          (L, (q0, r0))
   1.132 +        else
   1.133 +          raise (Floating_point "find_most_significant")
   1.134 +      end
   1.135 +  end
   1.136 +
   1.137 +fun approx_dec_by_bin n (q,r) =
   1.138 +  let
   1.139 +    fun addseq acc d' [] = acc
   1.140 +      | addseq acc d' (d::ds) = addseq (acc + exp2 (d - d')) d' ds
   1.141 +
   1.142 +    fun seq2bin [] = (0, 0)
   1.143 +      | seq2bin (d::ds) = (addseq 0 d ds + 1, d)
   1.144 +
   1.145 +    fun approx d_seq d0 precision (q,r) =
   1.146 +      if q = 0 then
   1.147 +        let val x = seq2bin d_seq in
   1.148 +          (x, x)
   1.149 +        end
   1.150 +      else
   1.151 +        let
   1.152 +          val (d, (q', r')) = find_most_significant q r
   1.153 +        in
   1.154 +          if precision < d0 - d then
   1.155 +            let
   1.156 +              val d' = d0 - precision
   1.157 +              val x1 = seq2bin (d_seq)
   1.158 +              val x2 = (fst x1 * exp2 (snd x1 - d') + 1,  d') (* = seq2bin (d'::d_seq) *)
   1.159 +            in
   1.160 +              (x1, x2)
   1.161 +            end
   1.162 +          else
   1.163 +            approx (d::d_seq) d0 precision (q', r')
   1.164 +        end
   1.165 +
   1.166 +    fun approx_start precision (q, r) =
   1.167 +      if q = 0 then
   1.168 +        ((0, 0), (0, 0))
   1.169 +      else
   1.170 +        let
   1.171 +          val (d, (q', r')) = find_most_significant q r
   1.172 +        in
   1.173 +          if precision <= 0 then
   1.174 +            let
   1.175 +              val x1 = seq2bin [d]
   1.176 +            in
   1.177 +              if q' = 0 then
   1.178 +                (x1, x1)
   1.179 +              else
   1.180 +                (x1, seq2bin [d + 1])
   1.181 +            end
   1.182 +          else
   1.183 +            approx [d] d precision (q', r')
   1.184 +        end
   1.185 +  in
   1.186 +    if 0 <= q then
   1.187 +      approx_start n (q,r)
   1.188 +    else
   1.189 +      let
   1.190 +        val ((a1,b1), (a2, b2)) = approx_start n (~ q, r)
   1.191 +      in
   1.192 +        ((~ a2, b2), (~ a1, b1))
   1.193 +      end
   1.194 +  end
   1.195 +
   1.196 +fun approx_decstr_by_bin n decstr =
   1.197 +  let
   1.198 +    fun str2int s = the_default 0 (Int.fromString s)
   1.199 +    fun signint p x = if p then x else ~ x
   1.200 +
   1.201 +    val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr
   1.202 +    val s = size d2
   1.203 +
   1.204 +    val q = signint p (str2int d1 * exp10 s + str2int d2)
   1.205 +    val r = signint ep (str2int e) - s
   1.206 +  in
   1.207 +    approx_dec_by_bin n (q,r)
   1.208 +  end
   1.209 +
   1.210 +fun mk_float (a, b) = @{term "float"} $
   1.211 +  HOLogic.mk_prod (pairself (HOLogic.mk_number HOLogic.intT) (a, b));
   1.212 +
   1.213 +fun dest_float (Const ("Float.float", _) $ (Const ("Pair", _) $ a $ b)) =
   1.214 +      pairself (snd o HOLogic.dest_number) (a, b)
   1.215 +  | dest_float t = ((snd o HOLogic.dest_number) t, 0);
   1.216 +
   1.217 +fun approx_float prec f value =
   1.218 +  let
   1.219 +    val interval = approx_decstr_by_bin prec value
   1.220 +    val (flower, fupper) = f interval
   1.221 +  in
   1.222 +    (mk_float flower, mk_float fupper)
   1.223 +  end;
   1.224 +
   1.225 +end;