src/HOL/Matrix/ExactFloatingPoint.ML
author obua
Fri Sep 03 17:10:36 2004 +0200 (2004-09-03)
changeset 15178 5f621aa35c25
permissions -rw-r--r--
Matrix theory, linear programming
obua@15178
     1
structure ExactFloatingPoint :
obua@15178
     2
sig
obua@15178
     3
    exception Destruct_floatstr of string
obua@15178
     4
				   
obua@15178
     5
    val destruct_floatstr : (char -> bool) -> (char -> bool) -> string -> bool * string * string * bool * string
obua@15178
     6
									  
obua@15178
     7
    exception Floating_point of string
obua@15178
     8
				
obua@15178
     9
    type floatrep = IntInf.int * IntInf.int
obua@15178
    10
    val approx_dec_by_bin : IntInf.int -> floatrep -> floatrep * floatrep
obua@15178
    11
    val approx_decstr_by_bin : int -> string -> floatrep * floatrep
obua@15178
    12
end 
obua@15178
    13
=
obua@15178
    14
struct
obua@15178
    15
obua@15178
    16
fun fst (a,b) = a
obua@15178
    17
fun snd (a,b) = b
obua@15178
    18
obua@15178
    19
val filter = List.filter;
obua@15178
    20
obua@15178
    21
exception Destruct_floatstr of string;
obua@15178
    22
obua@15178
    23
fun destruct_floatstr isDigit isExp number = 
obua@15178
    24
    let
obua@15178
    25
	val numlist = filter (not o Char.isSpace) (String.explode number)
obua@15178
    26
	
obua@15178
    27
	fun countsigns ((#"+")::cs) = countsigns cs
obua@15178
    28
	  | countsigns ((#"-")::cs) = 
obua@15178
    29
	    let	
obua@15178
    30
		val (positive, rest) = countsigns cs 
obua@15178
    31
	    in
obua@15178
    32
		(not positive, rest)
obua@15178
    33
	    end
obua@15178
    34
	  | countsigns cs = (true, cs)
obua@15178
    35
obua@15178
    36
	fun readdigits [] = ([], [])
obua@15178
    37
	  | readdigits (q as c::cs) = 
obua@15178
    38
	    if (isDigit c) then 
obua@15178
    39
		let
obua@15178
    40
		    val (digits, rest) = readdigits cs
obua@15178
    41
		in
obua@15178
    42
		    (c::digits, rest)
obua@15178
    43
		end
obua@15178
    44
	    else
obua@15178
    45
		([], q)		
obua@15178
    46
obua@15178
    47
	fun readfromexp_helper cs =
obua@15178
    48
	    let
obua@15178
    49
		val (positive, rest) = countsigns cs
obua@15178
    50
		val (digits, rest') = readdigits rest
obua@15178
    51
	    in
obua@15178
    52
		case rest' of
obua@15178
    53
		    [] => (positive, digits)
obua@15178
    54
		  | _ => raise (Destruct_floatstr number)
obua@15178
    55
	    end	    
obua@15178
    56
obua@15178
    57
	fun readfromexp [] = (true, [])
obua@15178
    58
	  | readfromexp (c::cs) = 
obua@15178
    59
	    if isExp c then
obua@15178
    60
		readfromexp_helper cs
obua@15178
    61
	    else 
obua@15178
    62
		raise (Destruct_floatstr number)		
obua@15178
    63
obua@15178
    64
	fun readfromdot [] = ([], readfromexp [])
obua@15178
    65
	  | readfromdot ((#".")::cs) = 
obua@15178
    66
	    let		
obua@15178
    67
		val (digits, rest) = readdigits cs
obua@15178
    68
		val exp = readfromexp rest
obua@15178
    69
	    in
obua@15178
    70
		(digits, exp)
obua@15178
    71
	    end		
obua@15178
    72
	  | readfromdot cs = readfromdot ((#".")::cs)
obua@15178
    73
			    
obua@15178
    74
	val (positive, numlist) = countsigns numlist				 
obua@15178
    75
	val (digits1, numlist) = readdigits numlist				 
obua@15178
    76
 	val (digits2, exp) = readfromdot numlist
obua@15178
    77
    in
obua@15178
    78
	(positive, String.implode digits1, String.implode digits2, fst exp, String.implode (snd exp))
obua@15178
    79
    end
obua@15178
    80
obua@15178
    81
type floatrep = IntInf.int * IntInf.int
obua@15178
    82
obua@15178
    83
exception Floating_point of string;
obua@15178
    84
obua@15178
    85
val ln2_10 = (Math.ln 10.0)/(Math.ln 2.0)
obua@15178
    86
	
obua@15178
    87
fun intmul a b = IntInf.* (a,b)
obua@15178
    88
fun intsub a b = IntInf.- (a,b)	
obua@15178
    89
fun intadd a b = IntInf.+ (a,b) 		 
obua@15178
    90
fun intpow a b = IntInf.pow (a, IntInf.toInt b);
obua@15178
    91
fun intle a b = IntInf.<= (a, b);
obua@15178
    92
fun intless a b = IntInf.< (a, b);
obua@15178
    93
fun intneg a = IntInf.~ a;
obua@15178
    94
val zero = IntInf.fromInt 0;
obua@15178
    95
val one = IntInf.fromInt 1;
obua@15178
    96
val two = IntInf.fromInt 2;
obua@15178
    97
val ten = IntInf.fromInt 10;
obua@15178
    98
val five = IntInf.fromInt 5;
obua@15178
    99
obua@15178
   100
fun find_most_significant q r = 
obua@15178
   101
    let 
obua@15178
   102
	fun int2real i = 
obua@15178
   103
	    case Real.fromString (IntInf.toString i) of 
obua@15178
   104
		SOME r => r 
obua@15178
   105
	      | NONE => raise (Floating_point "int2real")	
obua@15178
   106
	fun subtract (q, r) (q', r') = 
obua@15178
   107
	    if intle r r' then
obua@15178
   108
		(intsub q (intmul q' (intpow ten (intsub r' r))), r)
obua@15178
   109
	    else
obua@15178
   110
		(intsub (intmul q (intpow ten (intsub r r'))) q', r')
obua@15178
   111
	fun bin2dec d =
obua@15178
   112
	    if intle zero d then 
obua@15178
   113
		(intpow two d, zero)
obua@15178
   114
	    else
obua@15178
   115
		(intpow five (intneg d), d)				
obua@15178
   116
		
obua@15178
   117
	val L = IntInf.fromInt (Real.floor (int2real (IntInf.fromInt (IntInf.log2 q)) + (int2real r) * ln2_10))	
obua@15178
   118
	val L1 = intadd L one
obua@15178
   119
obua@15178
   120
	val (q1, r1) = subtract (q, r) (bin2dec L1) 		
obua@15178
   121
    in
obua@15178
   122
	if intle zero q1 then 
obua@15178
   123
	    let
obua@15178
   124
		val (q2, r2) = subtract (q, r) (bin2dec (intadd L1 one))
obua@15178
   125
	    in
obua@15178
   126
		if intle zero q2 then 
obua@15178
   127
		    raise (Floating_point "find_most_significant")
obua@15178
   128
		else
obua@15178
   129
		    (L1, (q1, r1))
obua@15178
   130
	    end
obua@15178
   131
	else
obua@15178
   132
	    let
obua@15178
   133
		val (q0, r0) = subtract (q, r) (bin2dec L)
obua@15178
   134
	    in
obua@15178
   135
		if intle zero q0 then
obua@15178
   136
		    (L, (q0, r0))
obua@15178
   137
		else
obua@15178
   138
		    raise (Floating_point "find_most_significant")
obua@15178
   139
	    end		    
obua@15178
   140
    end
obua@15178
   141
obua@15178
   142
fun approx_dec_by_bin n (q,r) =
obua@15178
   143
    let	
obua@15178
   144
	fun addseq acc d' [] = acc
obua@15178
   145
	  | addseq acc d' (d::ds) = addseq (intadd acc (intpow two (intsub d d'))) d' ds
obua@15178
   146
obua@15178
   147
	fun seq2bin [] = (zero, zero)
obua@15178
   148
	  | seq2bin (d::ds) = (intadd (addseq zero d ds) one, d)
obua@15178
   149
obua@15178
   150
	fun approx d_seq d0 precision (q,r) = 
obua@15178
   151
	    if q = zero then 
obua@15178
   152
		let val x = seq2bin d_seq in
obua@15178
   153
		    (x, x)
obua@15178
   154
		end
obua@15178
   155
	    else    
obua@15178
   156
		let 
obua@15178
   157
		    val (d, (q', r')) = find_most_significant q r
obua@15178
   158
		in	
obua@15178
   159
		    if intless precision (intsub d0 d) then 
obua@15178
   160
			let 
obua@15178
   161
			    val d' = intsub d0 precision
obua@15178
   162
			    val x1 = seq2bin (d_seq)
obua@15178
   163
			    val x2 = (intadd (intmul (fst x1) (intpow two (intsub (snd x1) d'))) one,  d') (* = seq2bin (d'::d_seq) *) 
obua@15178
   164
			in
obua@15178
   165
			    (x1, x2)
obua@15178
   166
			end
obua@15178
   167
		    else
obua@15178
   168
			approx (d::d_seq) d0 precision (q', r') 						    		
obua@15178
   169
		end		
obua@15178
   170
	
obua@15178
   171
	fun approx_start precision (q, r) =
obua@15178
   172
	    if q = zero then 
obua@15178
   173
		((zero, zero), (zero, zero))
obua@15178
   174
	    else
obua@15178
   175
		let 
obua@15178
   176
		    val (d, (q', r')) = find_most_significant q r
obua@15178
   177
		in	
obua@15178
   178
		    if intle precision zero then 
obua@15178
   179
			let
obua@15178
   180
			    val x1 = seq2bin [d]
obua@15178
   181
			in
obua@15178
   182
			    if q' = zero then 
obua@15178
   183
				(x1, x1)
obua@15178
   184
			    else
obua@15178
   185
				(x1, seq2bin [intadd d one])
obua@15178
   186
			end
obua@15178
   187
		    else
obua@15178
   188
			approx [d] d precision (q', r')
obua@15178
   189
		end		
obua@15178
   190
    in
obua@15178
   191
	if intle zero q then 
obua@15178
   192
	    approx_start n (q,r)
obua@15178
   193
	else
obua@15178
   194
	    let 
obua@15178
   195
		val ((a1,b1), (a2, b2)) = approx_start n (intneg q, r) 
obua@15178
   196
	    in
obua@15178
   197
		((intneg a2, b2), (intneg a1, b1))
obua@15178
   198
	    end					
obua@15178
   199
    end
obua@15178
   200
obua@15178
   201
fun approx_decstr_by_bin n decstr =
obua@15178
   202
    let 
obua@15178
   203
	fun str2int s = case IntInf.fromString s of SOME x => x | NONE => zero 
obua@15178
   204
	fun signint p x = if p then x else intneg x
obua@15178
   205
obua@15178
   206
	val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr
obua@15178
   207
	val s = IntInf.fromInt (size d2)
obua@15178
   208
		
obua@15178
   209
	val q = signint p (intadd (intmul (str2int d1) (intpow ten s)) (str2int d2))
obua@15178
   210
	val r = intsub (signint ep (str2int e)) s
obua@15178
   211
    in
obua@15178
   212
	approx_dec_by_bin (IntInf.fromInt n) (q,r)
obua@15178
   213
    end
obua@15178
   214
obua@15178
   215
end;