Matrix theory, linear programming
authorobua
Fri Sep 03 17:10:36 2004 +0200 (2004-09-03)
changeset 151785f621aa35c25
parent 15177 e7616269fdca
child 15179 b8ef323fd41e
Matrix theory, linear programming
src/HOL/IsaMakefile
src/HOL/Matrix/Cplex.ML
src/HOL/Matrix/CplexMatrixConverter.ML
src/HOL/Matrix/ExactFloatingPoint.ML
src/HOL/Matrix/Float.thy
src/HOL/Matrix/FloatArith.ML
src/HOL/Matrix/LinProg.thy
src/HOL/Matrix/Matrix.thy
src/HOL/Matrix/MatrixGeneral.thy
src/HOL/Matrix/MatrixLP.ML
src/HOL/Matrix/MatrixLP.thy
src/HOL/Matrix/MatrixLP_gensimp.ML
src/HOL/Matrix/ROOT.ML
src/HOL/Matrix/SparseMatrix.thy
src/HOL/Matrix/codegen_prep.ML
src/HOL/Matrix/conv.ML
src/HOL/Matrix/eq_codegen.ML
src/HOL/Matrix/fspmlp.ML
src/HOL/OrderedGroup.thy
src/HOL/Ring_and_Field.thy
     1.1 --- a/src/HOL/IsaMakefile	Fri Sep 03 10:27:05 2004 +0200
     1.2 +++ b/src/HOL/IsaMakefile	Fri Sep 03 17:10:36 2004 +0200
     1.3 @@ -1,5 +1,3 @@
     1.4 -#
     1.5 -# $Id$
     1.6  #
     1.7  # IsaMakefile for HOL
     1.8  #
     1.9 @@ -80,8 +78,7 @@
    1.10    $(SRC)/Provers/quasi.ML \
    1.11    $(SRC)/Provers/simplifier.ML $(SRC)/Provers/splitter.ML \
    1.12    $(SRC)/Provers/trancl.ML \
    1.13 -  $(SRC)/TFL/isand.ML $(SRC)/TFL/casesplit.ML $(SRC)/TFL/dcterm.ML\
    1.14 -  $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \
    1.15 +  $(SRC)/TFL/dcterm.ML $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \
    1.16    $(SRC)/TFL/tfl.ML $(SRC)/TFL/thms.ML $(SRC)/TFL/thry.ML \
    1.17    $(SRC)/TFL/usyntax.ML $(SRC)/TFL/utils.ML \
    1.18    Datatype.thy Datatype_Universe.ML Datatype_Universe.thy \
    1.19 @@ -97,8 +94,7 @@
    1.20    Nat.thy NatArith.ML NatArith.thy \
    1.21    Power.thy PreList.thy Product_Type.ML Product_Type.thy \
    1.22    Refute.thy ROOT.ML \
    1.23 -  Recdef.thy Reconstruction.thy Record.thy\
    1.24 -  Relation.ML Relation.thy Relation_Power.ML \
    1.25 +  Recdef.thy Record.thy Relation.ML Relation.thy Relation_Power.ML \
    1.26    Relation_Power.thy LOrder.thy OrderedGroup.thy OrderedGroup.ML Ring_and_Field.thy\
    1.27    Set.ML Set.thy SetInterval.ML SetInterval.thy \
    1.28    Sum_Type.ML Sum_Type.thy Tools/datatype_abs_proofs.ML Tools/datatype_aux.ML \
    1.29 @@ -109,7 +105,7 @@
    1.30    Tools/primrec_package.ML \
    1.31    Tools/prop_logic.ML \
    1.32    Tools/recdef_package.ML Tools/recfun_codegen.ML \
    1.33 -  Tools/record_package.ML Tools/reconstruction.ML\
    1.34 +  Tools/record_package.ML \
    1.35    Tools/refute.ML Tools/refute_isar.ML \
    1.36    Tools/rewrite_hol_proof.ML \
    1.37    Tools/sat_solver.ML \
    1.38 @@ -635,14 +631,18 @@
    1.39  
    1.40  ## HOL-Matrix
    1.41  
    1.42 -HOL-Matrix: HOL $(LOG)/HOL-Matrix.gz
    1.43 +HOL-Matrix: HOL HOL-Complex $(LOG)/HOL-Matrix.gz
    1.44  
    1.45 -$(LOG)/HOL-Matrix.gz: $(OUT)/HOL Matrix/ROOT.ML \
    1.46 -  Matrix/Matrix.thy Matrix/LinProg.thy Matrix/MatrixGeneral.thy \
    1.47 -  Matrix/SparseMatrix.thy Matrix/document/root.tex
    1.48 -	@$(ISATOOL) usedir $(OUT)/HOL Matrix
    1.49 -
    1.50 -
    1.51 +$(LOG)/HOL-Matrix.gz: $(OUT)/HOL-Complex \
    1.52 +  Matrix/MatrixGeneral.thy Matrix/Matrix.thy Matrix/SparseMatrix.thy \
    1.53 +  Matrix/Float.thy Matrix/FloatArith.ML Matrix/ExactFloatingPoint.ML \
    1.54 +  Matrix/Cplex.ML Matrix/CplexMatrixConverter.ML \
    1.55 +  Matrix/FloatSparseMatrixBuilder.ML \
    1.56 +  Matrix/conv.ML Matrix/eq_codegen.ML Matrix/codegen_prep.ML \
    1.57 +  Matrix/fspmlp.ML \
    1.58 +  Matrix/MatrixLP_gensimp.ML Matrix/MatrixLP.ML Matrix/MatrixLP.thy \
    1.59 +  Matrix/document/root.tex Matrix/ROOT.ML
    1.60 +	@cd Matrix; $(ISATOOL) usedir -b $(OUT)/HOL-Complex HOL-Matrix
    1.61  
    1.62  ## TLA
    1.63  
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/src/HOL/Matrix/Cplex.ML	Fri Sep 03 17:10:36 2004 +0200
     2.3 @@ -0,0 +1,1027 @@
     2.4 +signature CPLEX =
     2.5 +sig
     2.6 +    
     2.7 +    datatype cplexTerm = cplexVar of string | cplexNum of string | cplexInf 
     2.8 +                       | cplexNeg of cplexTerm 
     2.9 +                       | cplexProd of cplexTerm * cplexTerm 
    2.10 +                       | cplexSum of (cplexTerm list) 
    2.11 +    
    2.12 +    datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq 
    2.13 +								  
    2.14 +    datatype cplexGoal = cplexMinimize of cplexTerm 
    2.15 +		       | cplexMaximize of cplexTerm
    2.16 +					  
    2.17 +    datatype cplexConstr = cplexConstr of cplexComp * 
    2.18 +					  (cplexTerm * cplexTerm)
    2.19 +					  
    2.20 +    datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm 
    2.21 +					  * cplexComp * cplexTerm
    2.22 +			 | cplexBound of cplexTerm * cplexComp * cplexTerm
    2.23 +					 
    2.24 +    datatype cplexProg = cplexProg of string 
    2.25 +				      * cplexGoal 
    2.26 +				      * ((string option * cplexConstr) 
    2.27 +					     list) 
    2.28 +				      * cplexBounds list
    2.29 +				      
    2.30 +    datatype cplexResult = Unbounded 
    2.31 +			 | Infeasible 
    2.32 +			 | Undefined
    2.33 +			 | Optimal of string * 
    2.34 +				      (((* name *) string * 
    2.35 +					(* value *) string) list)
    2.36 +						     
    2.37 +    exception Load_cplexFile of string
    2.38 +    exception Load_cplexResult of string
    2.39 +    exception Save_cplexFile of string
    2.40 +    exception Execute of string
    2.41 +
    2.42 +    val load_cplexFile : string -> cplexProg
    2.43 +
    2.44 +    val save_cplexFile : string -> cplexProg -> unit 
    2.45 +
    2.46 +    val elim_nonfree_bounds : cplexProg -> cplexProg
    2.47 +
    2.48 +    val relax_strict_ineqs : cplexProg -> cplexProg
    2.49 +
    2.50 +    val is_normed_cplexProg : cplexProg -> bool
    2.51 +
    2.52 +    val load_cplexResult : string -> cplexResult
    2.53 +
    2.54 +    val solve : cplexProg -> cplexResult
    2.55 +end;
    2.56 +
    2.57 +structure Cplex : CPLEX =
    2.58 +struct
    2.59 +
    2.60 +exception Load_cplexFile of string;
    2.61 +exception Load_cplexResult of string;
    2.62 +exception Save_cplexFile of string;
    2.63 +	  
    2.64 +datatype cplexTerm = cplexVar of string 
    2.65 +		   | cplexNum of string 
    2.66 +		   | cplexInf 
    2.67 +                   | cplexNeg of cplexTerm 
    2.68 +                   | cplexProd of cplexTerm * cplexTerm 
    2.69 +                   | cplexSum of (cplexTerm list) 
    2.70 +datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq 
    2.71 +datatype cplexGoal = cplexMinimize of cplexTerm | cplexMaximize of cplexTerm
    2.72 +datatype cplexConstr = cplexConstr of cplexComp * (cplexTerm * cplexTerm)
    2.73 +datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm 
    2.74 +				      * cplexComp * cplexTerm
    2.75 +                     | cplexBound of cplexTerm * cplexComp * cplexTerm 
    2.76 +datatype cplexProg = cplexProg of string 
    2.77 +				  * cplexGoal 
    2.78 +				  * ((string option * cplexConstr) list) 
    2.79 +				  * cplexBounds list
    2.80 +				  
    2.81 +fun rev_cmp cplexLe = cplexGe
    2.82 +  | rev_cmp cplexLeq = cplexGeq
    2.83 +  | rev_cmp cplexGe = cplexLe
    2.84 +  | rev_cmp cplexGeq = cplexLeq
    2.85 +  | rev_cmp cplexEq = cplexEq
    2.86 +
    2.87 +fun the None = raise (Load_cplexFile "Some expected")
    2.88 +  | the (Some x) = x; 
    2.89 +    
    2.90 +fun modulo_signed is_something (cplexNeg u) = is_something u
    2.91 +  | modulo_signed is_something u = is_something u
    2.92 +
    2.93 +fun is_Num (cplexNum _) = true
    2.94 +  | is_Num _ = false
    2.95 +
    2.96 +fun is_Inf cplexInf = true
    2.97 +  | is_Inf _ = false
    2.98 +
    2.99 +fun is_Var (cplexVar _) = true
   2.100 +  | is_Var _ = false
   2.101 +
   2.102 +fun is_Neg (cplexNeg x ) = true
   2.103 +  | is_Neg _ = false
   2.104 +
   2.105 +fun is_normed_Prod (cplexProd (t1, t2)) = 
   2.106 +    (is_Num t1) andalso (is_Var t2)
   2.107 +  | is_normed_Prod x = is_Var x
   2.108 +
   2.109 +fun is_normed_Sum (cplexSum ts) = 
   2.110 +    (ts <> []) andalso forall (modulo_signed is_normed_Prod) ts
   2.111 +  | is_normed_Sum x = modulo_signed is_normed_Prod x
   2.112 +
   2.113 +fun is_normed_Constr (cplexConstr (c, (t1, t2))) = 
   2.114 +    (is_normed_Sum t1) andalso (modulo_signed is_Num t2)
   2.115 +
   2.116 +fun is_Num_or_Inf x = is_Inf x orelse is_Num x
   2.117 +
   2.118 +fun is_normed_Bounds (cplexBounds (t1, c1, t2, c2, t3)) =
   2.119 +    (c1 = cplexLe orelse c1 = cplexLeq) andalso 
   2.120 +    (c2 = cplexLe orelse c2 = cplexLeq) andalso
   2.121 +    is_Var t2 andalso
   2.122 +    modulo_signed is_Num_or_Inf t1 andalso
   2.123 +    modulo_signed is_Num_or_Inf t3
   2.124 +  | is_normed_Bounds (cplexBound (t1, c, t2)) =
   2.125 +    (is_Var t1 andalso (modulo_signed is_Num_or_Inf t2)) 
   2.126 +    orelse
   2.127 +    (c <> cplexEq andalso 
   2.128 +     is_Var t2 andalso (modulo_signed is_Num_or_Inf t1))
   2.129 +	             		
   2.130 +fun term_of_goal (cplexMinimize x) = x
   2.131 +  | term_of_goal (cplexMaximize x) = x
   2.132 +
   2.133 +fun is_normed_cplexProg (cplexProg (name, goal, constraints, bounds)) = 
   2.134 +    is_normed_Sum (term_of_goal goal) andalso
   2.135 +    forall (fn (_,x) => is_normed_Constr x) constraints andalso
   2.136 +    forall is_normed_Bounds bounds
   2.137 +
   2.138 +fun is_NL s = s = "\n"
   2.139 +
   2.140 +fun is_blank s = forall (fn c => c <> #"\n" andalso Char.isSpace c) (String.explode s)
   2.141 +		 
   2.142 +fun is_num a = 
   2.143 +    let 
   2.144 +	val b = String.explode a 
   2.145 +	fun num4 cs = forall Char.isDigit cs		      
   2.146 +	fun num3 [] = true
   2.147 +	  | num3 (ds as (c::cs)) = 
   2.148 +	    if c = #"+" orelse c = #"-" then
   2.149 +		num4 cs
   2.150 +	    else 
   2.151 +		num4 ds		    
   2.152 +	fun num2 [] = true
   2.153 +	  | num2 (c::cs) =
   2.154 +	    if c = #"e" orelse c = #"E" then num3 cs
   2.155 +	    else (Char.isDigit c) andalso num2 cs			 
   2.156 +	fun num1 [] = true
   2.157 +	  | num1 (c::cs) = 
   2.158 +	    if c = #"." then num2 cs 
   2.159 +	    else if c = #"e" orelse c = #"E" then num3 cs 
   2.160 +	    else (Char.isDigit c) andalso num1 cs
   2.161 +	fun num [] = true
   2.162 +	  | num (c::cs) = 
   2.163 +	    if c = #"." then num2 cs 
   2.164 +	    else (Char.isDigit c) andalso num1 cs
   2.165 +    in
   2.166 +	num b
   2.167 +    end  
   2.168 +
   2.169 +fun is_delimiter s = s = "+" orelse s = "-" orelse s = ":"
   2.170 +		    
   2.171 +fun is_cmp s = s = "<" orelse s = ">" orelse s = "<=" 
   2.172 +		     orelse s = ">=" orelse s = "="
   2.173 +
   2.174 +fun is_symbol a = 
   2.175 +    let
   2.176 +	val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~"
   2.177 +	fun is_symbol_char c = Char.isAlphaNum c orelse 
   2.178 +			       exists (fn d => d=c) symbol_char
   2.179 +	fun is_symbol_start c = is_symbol_char c andalso 
   2.180 +				not (Char.isDigit c) andalso 
   2.181 +				not (c= #".")
   2.182 +	val b = String.explode a
   2.183 +    in
   2.184 +	b <> [] andalso is_symbol_start (hd b) andalso 
   2.185 +	forall is_symbol_char b
   2.186 +    end
   2.187 +
   2.188 +fun to_upper s = String.implode (map Char.toUpper (String.explode s))
   2.189 +
   2.190 +fun keyword x = 
   2.191 +    let 	
   2.192 +	val a = to_upper x 
   2.193 +    in
   2.194 +	if a = "BOUNDS" orelse a = "BOUND" then
   2.195 +	    Some "BOUNDS"
   2.196 +	else if a = "MINIMIZE" orelse a = "MINIMUM" orelse a = "MIN" then
   2.197 +	    Some "MINIMIZE"
   2.198 +	else if a = "MAXIMIZE" orelse a = "MAXIMUM" orelse a = "MAX" then
   2.199 +	    Some "MAXIMIZE"
   2.200 +	else if a = "ST" orelse a = "S.T." orelse a = "ST." then
   2.201 +	    Some "ST"
   2.202 +	else if a = "FREE" orelse a = "END" then
   2.203 +	    Some a
   2.204 +	else if a = "GENERAL" orelse a = "GENERALS" orelse a = "GEN" then
   2.205 +	    Some "GENERAL"
   2.206 +	else if a = "INTEGER" orelse a = "INTEGERS" orelse a = "INT" then
   2.207 +	    Some "INTEGER"
   2.208 +	else if a = "BINARY" orelse a = "BINARIES" orelse a = "BIN" then
   2.209 +	    Some "BINARY"
   2.210 +	else if a = "INF" orelse a = "INFINITY" then
   2.211 +	    Some "INF"
   2.212 +	else								   
   2.213 +	    None
   2.214 +    end
   2.215 +	        
   2.216 +val TOKEN_ERROR = ~1
   2.217 +val TOKEN_BLANK = 0
   2.218 +val TOKEN_NUM = 1
   2.219 +val TOKEN_DELIMITER = 2
   2.220 +val TOKEN_SYMBOL = 3
   2.221 +val TOKEN_LABEL = 4
   2.222 +val TOKEN_CMP = 5
   2.223 +val TOKEN_KEYWORD = 6
   2.224 +val TOKEN_NL = 7
   2.225 +		      
   2.226 +(* tokenize takes a list of chars as argument and returns a list of
   2.227 +   int * string pairs, each string representing a "cplex token", 
   2.228 +   and each int being one of TOKEN_NUM, TOKEN_DELIMITER, TOKEN_CMP
   2.229 +   or TOKEN_SYMBOL *)
   2.230 +fun tokenize s = 
   2.231 +    let
   2.232 +	val flist = [(is_NL, TOKEN_NL), 
   2.233 +		     (is_blank, TOKEN_BLANK), 
   2.234 +		     (is_num, TOKEN_NUM),  
   2.235 +                     (is_delimiter, TOKEN_DELIMITER), 
   2.236 +		     (is_cmp, TOKEN_CMP), 
   2.237 +		     (is_symbol, TOKEN_SYMBOL)]  
   2.238 +	fun match_helper [] s = (fn x => false, TOKEN_ERROR)
   2.239 +	  | match_helper (f::fs) s = 
   2.240 +	    if ((fst f) s) then f else match_helper fs s   
   2.241 +	fun match s = match_helper flist s
   2.242 +	fun tok s = 
   2.243 +	    if s = "" then [] else       
   2.244 +	    let 
   2.245 +		val h = String.substring (s,0,1) 
   2.246 +		val (f, j) = match h
   2.247 +		fun len i = 
   2.248 +		    if size s = i then i 
   2.249 +		    else if f (String.substring (s,0,i+1)) then 
   2.250 +			len (i+1) 
   2.251 +		    else i 
   2.252 +	    in
   2.253 +		if j < 0 then 
   2.254 +		    (if h = "\\" then [] 
   2.255 +		     else raise (Load_cplexFile ("token expected, found: "
   2.256 +						 ^s)))
   2.257 +		else 
   2.258 +		    let 
   2.259 +			val l = len 1 
   2.260 +			val u = String.substring (s,0,l)
   2.261 +			val v = String.extract (s,l,NONE)
   2.262 +		    in
   2.263 +			if j = 0 then tok v else (j, u) :: tok v
   2.264 +		    end
   2.265 +	    end                
   2.266 +    in
   2.267 +	tok s     
   2.268 +    end
   2.269 +
   2.270 +exception Tokenize of string;
   2.271 +
   2.272 +fun tokenize_general flist s = 
   2.273 +    let
   2.274 +	fun match_helper [] s = raise (Tokenize s)
   2.275 +	  | match_helper (f::fs) s = 
   2.276 +	    if ((fst f) s) then f else match_helper fs s   
   2.277 +	fun match s = match_helper flist s
   2.278 +	fun tok s = 
   2.279 +	    if s = "" then [] else       
   2.280 +	    let 
   2.281 +		val h = String.substring (s,0,1) 
   2.282 +		val (f, j) = match h
   2.283 +		fun len i = 
   2.284 +		    if size s = i then i 
   2.285 +		    else if f (String.substring (s,0,i+1)) then 
   2.286 +			len (i+1) 
   2.287 +		    else i 
   2.288 +		val l = len 1
   2.289 +	    in
   2.290 +		(j, String.substring (s,0,l)) :: tok (String.extract (s,l,NONE))
   2.291 +	    end
   2.292 +    in
   2.293 +	tok s     
   2.294 +    end
   2.295 +
   2.296 +fun load_cplexFile name = 
   2.297 +    let 
   2.298 +	val f = TextIO.openIn name  
   2.299 +        val ignore_NL = ref true  
   2.300 +	val rest = ref []
   2.301 +
   2.302 +	fun is_symbol s c = (fst c) = TOKEN_SYMBOL andalso (to_upper (snd c)) = s 
   2.303 +		   
   2.304 +	fun readToken_helper () =
   2.305 +	    if length (!rest) > 0 then  
   2.306 +		let val u = hd (!rest) in 
   2.307 +		    (
   2.308 +		     rest := tl (!rest); 
   2.309 +		     Some u
   2.310 +		    ) 
   2.311 +		end
   2.312 +	    else 
   2.313 +		let val s = TextIO.inputLine f in
   2.314 +		    if s = "" then None else
   2.315 +		    let val t = tokenize s in
   2.316 +			if (length t >= 2 andalso 
   2.317 +			    snd(hd (tl t)) = ":") 
   2.318 +			then 
   2.319 +			    rest := (TOKEN_LABEL, snd (hd t)) :: (tl (tl t))
   2.320 +			else if (length t >= 2) andalso is_symbol "SUBJECT" (hd (t)) 
   2.321 +				andalso is_symbol "TO" (hd (tl t)) 
   2.322 +			then
   2.323 +			    rest := (TOKEN_SYMBOL, "ST") :: (tl (tl t))
   2.324 +			else
   2.325 +			    rest := t;
   2.326 +			readToken_helper ()  
   2.327 +		    end
   2.328 +		end           
   2.329 +	
   2.330 +	fun readToken_helper2 () = 
   2.331 +		let val c = readToken_helper () in
   2.332 +		    if c = None then None
   2.333 +                    else if !ignore_NL andalso fst (the c) = TOKEN_NL then
   2.334 +			readToken_helper2 ()
   2.335 +		    else if fst (the c) = TOKEN_SYMBOL 
   2.336 +			    andalso keyword (snd (the c)) <> None
   2.337 +		    then Some (TOKEN_KEYWORD, the (keyword (snd (the c))))
   2.338 +		    else c
   2.339 +		end
   2.340 +		
   2.341 +	fun readToken () = readToken_helper2 ()
   2.342 +	    			   
   2.343 +	fun pushToken a = rest := (a::(!rest))
   2.344 +	    
   2.345 +	fun is_value token = 
   2.346 +	    fst token = TOKEN_NUM orelse (fst token = TOKEN_KEYWORD 
   2.347 +					  andalso snd token = "INF")	
   2.348 +
   2.349 +        fun get_value token = 
   2.350 +	    if fst token = TOKEN_NUM then
   2.351 +		cplexNum (snd token)
   2.352 +	    else if fst token = TOKEN_KEYWORD andalso snd token = "INF" 
   2.353 +	    then
   2.354 +		cplexInf
   2.355 +	    else
   2.356 +		raise (Load_cplexFile "num expected")
   2.357 +					    
   2.358 +	fun readTerm_Product only_num = 
   2.359 +	    let val c = readToken () in
   2.360 +		if c = None then None
   2.361 +		else if fst (the c) = TOKEN_SYMBOL 
   2.362 +		then (
   2.363 +		    if only_num then (pushToken (the c); None) 
   2.364 +		    else Some (cplexVar (snd (the c)))
   2.365 +		    ) 
   2.366 +		else if only_num andalso is_value (the c) then
   2.367 +		    Some (get_value (the c))
   2.368 +		else if is_value (the c) then  
   2.369 +		    let val t1 = get_value (the c) 
   2.370 +			val d = readToken () 
   2.371 +		    in
   2.372 +			if d = None then Some t1 
   2.373 +			else if fst (the d) = TOKEN_SYMBOL then 
   2.374 +			    Some (cplexProd (t1, cplexVar (snd (the d))))
   2.375 +			else
   2.376 +			    (pushToken (the d); Some t1)
   2.377 +		    end
   2.378 +		else (pushToken (the c); None)
   2.379 +	    end
   2.380 +	    
   2.381 +	fun readTerm_Signed only_signed only_num = 
   2.382 +	    let
   2.383 +		val c = readToken ()
   2.384 +	    in
   2.385 +		if c = None then None 
   2.386 +		else 
   2.387 +		    let val d = the c in 
   2.388 +			if d = (TOKEN_DELIMITER, "+") then 
   2.389 +			    readTerm_Product only_num
   2.390 +			 else if d = (TOKEN_DELIMITER, "-") then 
   2.391 +			     Some (cplexNeg (the (readTerm_Product 
   2.392 +						      only_num)))
   2.393 +			 else (pushToken d; 
   2.394 +			       if only_signed then None 
   2.395 +			       else readTerm_Product only_num)
   2.396 +		    end
   2.397 +	    end
   2.398 +    
   2.399 +	fun readTerm_Sum first_signed = 
   2.400 +	    let val c = readTerm_Signed first_signed false in
   2.401 +		if c = None then [] else (the c)::(readTerm_Sum true)
   2.402 +	    end
   2.403 +	    
   2.404 +	fun readTerm () = 
   2.405 +	    let val c = readTerm_Sum false in
   2.406 +		if c = [] then None 
   2.407 +		else if tl c = [] then Some (hd c)
   2.408 +		else Some (cplexSum c) 
   2.409 +	    end
   2.410 +	    
   2.411 +	fun readLabeledTerm () = 
   2.412 +	    let val c = readToken () in
   2.413 +		if c = None then (None, None) 
   2.414 +		else if fst (the c) = TOKEN_LABEL then 
   2.415 +		    let val t = readTerm () in
   2.416 +			if t = None then 
   2.417 +			    raise (Load_cplexFile ("term after label "^
   2.418 +						   (snd (the c))^
   2.419 +						   " expected"))
   2.420 +			else (Some (snd (the c)), t)
   2.421 +		    end
   2.422 +		else (pushToken (the c); (None, readTerm ()))
   2.423 +	    end
   2.424 +            
   2.425 +	fun readGoal () = 
   2.426 +	    let 
   2.427 +		val g = readToken ()
   2.428 +	    in
   2.429 +    		if g = Some (TOKEN_KEYWORD, "MAXIMIZE") then 
   2.430 +		    cplexMaximize (the (snd (readLabeledTerm ()))) 
   2.431 +		else if g = Some (TOKEN_KEYWORD, "MINIMIZE") then 
   2.432 +		    cplexMinimize (the (snd (readLabeledTerm ())))
   2.433 +		else raise (Load_cplexFile "MAXIMIZE or MINIMIZE expected")
   2.434 +	    end
   2.435 +	    
   2.436 +	fun str2cmp b = 
   2.437 +	    (case b of 
   2.438 +		 "<" => cplexLe 
   2.439 +	       | "<=" => cplexLeq 
   2.440 +	       | ">" => cplexGe 
   2.441 +	       | ">=" => cplexGeq 
   2.442 +               | "=" => cplexEq
   2.443 +	       | _ => raise (Load_cplexFile (b^" is no TOKEN_CMP")))
   2.444 +			    
   2.445 +	fun readConstraint () = 
   2.446 +            let 
   2.447 +		val t = readLabeledTerm () 
   2.448 +		fun make_constraint b t1 t2 =
   2.449 +                    cplexConstr 
   2.450 +			(str2cmp b,
   2.451 +			 (t1, t2))		
   2.452 +	    in
   2.453 +		if snd t = None then None
   2.454 +		else 
   2.455 +		    let val c = readToken () in
   2.456 +			if c = None orelse fst (the c) <> TOKEN_CMP 
   2.457 +			then raise (Load_cplexFile "TOKEN_CMP expected")
   2.458 +			else 
   2.459 +			    let val n = readTerm_Signed false true in
   2.460 +				if n = None then 
   2.461 +				    raise (Load_cplexFile "num expected")
   2.462 +				else
   2.463 +				    Some (fst t,
   2.464 +					  make_constraint (snd (the c)) 
   2.465 +							  (the (snd t)) 
   2.466 +							  (the n))
   2.467 +			    end
   2.468 +		    end            
   2.469 +	    end
   2.470 +        
   2.471 +        fun readST () = 
   2.472 +	    let 
   2.473 +		fun readbody () = 
   2.474 +		    let val t = readConstraint () in
   2.475 +			if t = None then []
   2.476 +			else if (is_normed_Constr (snd (the t))) then
   2.477 +			    (the t)::(readbody ())
   2.478 +			else if (fst (the t) <> None) then 
   2.479 +			    raise (Load_cplexFile 
   2.480 +				       ("constraint '"^(the (fst (the t)))^ 
   2.481 +					"'is not normed"))
   2.482 +			else
   2.483 +			    raise (Load_cplexFile 
   2.484 +				       "constraint is not normed")
   2.485 +		    end
   2.486 +	    in    		
   2.487 +		if readToken () = Some (TOKEN_KEYWORD, "ST") 
   2.488 +		then 
   2.489 +		    readbody ()
   2.490 +		else 
   2.491 +		    raise (Load_cplexFile "ST expected")
   2.492 +	    end
   2.493 +	    
   2.494 +	fun readCmp () = 
   2.495 +	    let val c = readToken () in
   2.496 +		if c = None then None
   2.497 +		else if fst (the c) = TOKEN_CMP then 
   2.498 +		    Some (str2cmp (snd (the c)))
   2.499 +		else (pushToken (the c); None)
   2.500 +	    end
   2.501 +	    
   2.502 +	fun skip_NL () =
   2.503 +	    let val c = readToken () in
   2.504 +		if c <> None andalso fst (the c) = TOKEN_NL then
   2.505 +		    skip_NL ()
   2.506 +		else
   2.507 +		    (pushToken (the c); ())
   2.508 +	    end
   2.509 +
   2.510 +	fun is_var (cplexVar _) = true
   2.511 +	  | is_var _ = false
   2.512 +	    
   2.513 +	fun make_bounds c t1 t2 = 
   2.514 +	    cplexBound (t1, c, t2)
   2.515 +
   2.516 +	fun readBound () = 	    
   2.517 +	    let 		
   2.518 +		val _ = skip_NL ()
   2.519 +		val t1 = readTerm ()
   2.520 +	    in
   2.521 +		if t1 = None then None 
   2.522 +		else 
   2.523 +		    let 
   2.524 +			val c1 = readCmp () 
   2.525 +		    in
   2.526 +			if c1 = None then
   2.527 +			    let 
   2.528 +				val c = readToken () 
   2.529 +			    in
   2.530 +				if c = Some (TOKEN_KEYWORD, "FREE") then
   2.531 +				    Some (
   2.532 +				    cplexBounds (cplexNeg cplexInf,
   2.533 +						 cplexLeq,
   2.534 +						 the t1,
   2.535 +						 cplexLeq,
   2.536 +						 cplexInf))
   2.537 +				else
   2.538 +				    raise (Load_cplexFile "FREE expected")
   2.539 +			    end
   2.540 +			else
   2.541 +			    let 
   2.542 +				val t2 = readTerm () 
   2.543 +			    in
   2.544 +				if t2 = None then
   2.545 +				    raise (Load_cplexFile "term expected")
   2.546 +				else
   2.547 +				    let val c2 = readCmp () in
   2.548 +					if c2 = None then
   2.549 +					    Some (make_bounds (the c1) 
   2.550 +							      (the t1)
   2.551 +							      (the t2))
   2.552 +					else
   2.553 +					    Some (
   2.554 +					    cplexBounds (the t1,
   2.555 +							 the c1,
   2.556 +							 the t2,
   2.557 +							 the c2,
   2.558 +							 the (readTerm())))
   2.559 +				    end
   2.560 +			    end
   2.561 +		    end
   2.562 +	    end
   2.563 +	    
   2.564 +	fun readBounds () =
   2.565 +	    let 
   2.566 +		fun makestring b = "?"
   2.567 +		fun readbody () = 
   2.568 +		    let 
   2.569 +			val b = readBound ()
   2.570 +		    in
   2.571 +			if b = None then []
   2.572 +			else if (is_normed_Bounds (the b)) then
   2.573 +			    (the b)::(readbody())
   2.574 +			else (
   2.575 +			    raise (Load_cplexFile 
   2.576 +				       ("bounds are not normed in: "^
   2.577 +					(makestring (the b)))))
   2.578 +		    end
   2.579 +	    in
   2.580 +		if readToken () = Some (TOKEN_KEYWORD, "BOUNDS") then 
   2.581 +		    readbody ()
   2.582 +		else raise (Load_cplexFile "BOUNDS expected")
   2.583 +	    end
   2.584 +
   2.585 +        fun readEnd () = 
   2.586 +	    if readToken () = Some (TOKEN_KEYWORD, "END") then ()
   2.587 +	    else raise (Load_cplexFile "END expected")
   2.588 +		    
   2.589 +	val result_Goal = readGoal ()  
   2.590 +	val result_ST = readST ()  
   2.591 +	val _ =	ignore_NL := false 
   2.592 +        val result_Bounds = readBounds ()
   2.593 +        val _ = ignore_NL := true
   2.594 +        val _ = readEnd ()
   2.595 +	val _ = TextIO.closeIn f
   2.596 +    in
   2.597 +	cplexProg (name, result_Goal, result_ST, result_Bounds)
   2.598 +    end
   2.599 +
   2.600 +fun save_cplexFile filename (cplexProg (name, goal, constraints, bounds)) =
   2.601 +    let
   2.602 +	val f = TextIO.openOut filename
   2.603 +	
   2.604 +	fun basic_write s = TextIO.output(f, s)
   2.605 +
   2.606 +	val linebuf = ref ""
   2.607 +	fun buf_flushline s = 
   2.608 +	    (basic_write (!linebuf); 
   2.609 +	     basic_write "\n";
   2.610 +	     linebuf := s)
   2.611 +	fun buf_add s = linebuf := (!linebuf) ^ s
   2.612 +
   2.613 +	fun write s = 
   2.614 +	    if (String.size s) + (String.size (!linebuf)) >= 250 then 
   2.615 +		buf_flushline ("    "^s)
   2.616 +	    else 
   2.617 +		buf_add s
   2.618 +
   2.619 +        fun writeln s = (buf_add s; buf_flushline "")
   2.620 +        
   2.621 +	fun write_term (cplexVar x) = write x
   2.622 +	  | write_term (cplexNum x) = write x
   2.623 +	  | write_term cplexInf = write "inf"
   2.624 +	  | write_term (cplexProd (cplexNum "1", b)) = write_term b
   2.625 +	  | write_term (cplexProd (a, b)) = 
   2.626 +	    (write_term a; write " "; write_term b)
   2.627 +          | write_term (cplexNeg x) = (write " - "; write_term x)
   2.628 +          | write_term (cplexSum ts) = write_terms ts
   2.629 +	and write_terms [] = ()
   2.630 +	  | write_terms (t::ts) = 
   2.631 +	    (if (not (is_Neg t)) then write " + " else (); 
   2.632 +	     write_term t; write_terms ts)
   2.633 +
   2.634 +	fun write_goal (cplexMaximize term) = 
   2.635 +	    (writeln "MAXIMIZE"; write_term term; writeln "")
   2.636 +	  | write_goal (cplexMinimize term) = 
   2.637 +	    (writeln "MINIMIZE"; write_term term; writeln "")
   2.638 +
   2.639 +	fun write_cmp cplexLe = write "<"
   2.640 +	  | write_cmp cplexLeq = write "<="
   2.641 +	  | write_cmp cplexEq = write "="
   2.642 +	  | write_cmp cplexGe = write ">"
   2.643 +	  | write_cmp cplexGeq = write ">="
   2.644 +        
   2.645 +	fun write_constr (cplexConstr (cmp, (a,b))) = 
   2.646 +	    (write_term a;
   2.647 +	     write " ";
   2.648 +	     write_cmp cmp;
   2.649 +	     write " ";
   2.650 +	     write_term b)
   2.651 +        
   2.652 +	fun write_constraints [] = ()
   2.653 +	  | write_constraints (c::cs) = 
   2.654 +	    (if (fst c <> None) 
   2.655 +	     then 
   2.656 +		 (write (the (fst c)); write ": ") 
   2.657 +	     else 
   2.658 +		 ();
   2.659 +	     write_constr (snd c);
   2.660 +	     writeln "";
   2.661 +	     write_constraints cs)
   2.662 +
   2.663 +	fun write_bounds [] = ()
   2.664 +	  | write_bounds ((cplexBounds (t1,c1,t2,c2,t3))::bs) =
   2.665 +	    ((if t1 = cplexNeg cplexInf andalso t3 = cplexInf
   2.666 +		 andalso (c1 = cplexLeq orelse c1 = cplexLe) 
   2.667 +		 andalso (c2 = cplexLeq orelse c2 = cplexLe) 
   2.668 +	      then
   2.669 +		  (write_term t2; write " free")
   2.670 +	      else		 		 
   2.671 +		  (write_term t1; write " "; write_cmp c1; write " ";
   2.672 +		   write_term t2; write " "; write_cmp c2; write " ";
   2.673 +		   write_term t3)
   2.674 +	     ); writeln ""; write_bounds bs)
   2.675 +	  | write_bounds ((cplexBound (t1, c, t2)) :: bs) = 
   2.676 +	    (write_term t1; write " "; 
   2.677 +	     write_cmp c; write " ";
   2.678 +	     write_term t2; writeln ""; write_bounds bs)
   2.679 +	     
   2.680 +	val _ = write_goal goal
   2.681 +        val _ = (writeln ""; writeln "ST")
   2.682 +	val _ = write_constraints constraints
   2.683 +        val _ = (writeln ""; writeln "BOUNDS")
   2.684 +	val _ = write_bounds bounds
   2.685 +        val _ = (writeln ""; writeln "END") 
   2.686 +        val _ = TextIO.closeOut f
   2.687 +    in
   2.688 +	()
   2.689 +    end
   2.690 +
   2.691 +fun norm_Constr (constr as cplexConstr (c, (t1, t2))) = 
   2.692 +    if not (modulo_signed is_Num t2) andalso  
   2.693 +       modulo_signed is_Num t1 
   2.694 +    then
   2.695 +	[cplexConstr (rev_cmp c, (t2, t1))]
   2.696 +    else if (c = cplexLe orelse c = cplexLeq) andalso
   2.697 +	    (t1 = (cplexNeg cplexInf) orelse t2 = cplexInf)	    
   2.698 +    then 
   2.699 +	[]
   2.700 +    else if (c = cplexGe orelse c = cplexGeq) andalso
   2.701 +	    (t1 = cplexInf orelse t2 = cplexNeg cplexInf)
   2.702 +    then
   2.703 +	[]
   2.704 +    else
   2.705 +	[constr]
   2.706 +
   2.707 +fun bound2constr (cplexBounds (t1,c1,t2,c2,t3)) =
   2.708 +    (norm_Constr(cplexConstr (c1, (t1, t2))))
   2.709 +    @ (norm_Constr(cplexConstr (c2, (t2, t3))))
   2.710 +  | bound2constr (cplexBound (t1, cplexEq, t2)) =
   2.711 +    (norm_Constr(cplexConstr (cplexLeq, (t1, t2))))
   2.712 +    @ (norm_Constr(cplexConstr (cplexLeq, (t2, t1))))
   2.713 +  | bound2constr (cplexBound (t1, c1, t2)) = 
   2.714 +    norm_Constr(cplexConstr (c1, (t1,t2)))
   2.715 +
   2.716 +val emptyset = Symtab.empty
   2.717 +
   2.718 +fun singleton v = Symtab.update ((v, ()), emptyset)
   2.719 +
   2.720 +fun merge a b = Symtab.merge (op =) (a, b)
   2.721 +
   2.722 +fun mergemap f ts = foldl
   2.723 +			(fn (table, x) => merge table (f x))
   2.724 +			(Symtab.empty, ts)
   2.725 +
   2.726 +fun diff a b = Symtab.foldl (fn (a, (k,v)) => 
   2.727 +			 (Symtab.delete k a) handle UNDEF => a) 
   2.728 +		     (a, b)
   2.729 +		    
   2.730 +fun collect_vars (cplexVar v) = singleton v
   2.731 +  | collect_vars (cplexNeg t) = collect_vars t
   2.732 +  | collect_vars (cplexProd (t1, t2)) = 
   2.733 +    merge (collect_vars t1) (collect_vars t2)
   2.734 +  | collect_vars (cplexSum ts) = mergemap collect_vars ts
   2.735 +  | collect_vars _ = emptyset
   2.736 +
   2.737 +(* Eliminates all nonfree bounds from the linear program and produces an
   2.738 +   equivalent program with only free bounds 
   2.739 +   IF for the input program P holds: is_normed_cplexProg P *)
   2.740 +fun elim_nonfree_bounds (cplexProg (name, goal, constraints, bounds)) =
   2.741 +    let
   2.742 +	fun collect_constr_vars (_, cplexConstr (c, (t1,_))) =
   2.743 +	    (collect_vars t1)   
   2.744 +	
   2.745 +	val cvars = merge (collect_vars (term_of_goal goal)) 
   2.746 +			  (mergemap collect_constr_vars constraints)         
   2.747 +		    
   2.748 +	fun collect_lower_bounded_vars 
   2.749 +	    (cplexBounds (t1, c1, cplexVar v, c2, t3)) =
   2.750 +	    singleton v
   2.751 +	  |  collect_lower_bounded_vars 
   2.752 +		 (cplexBound (_, cplexLe, cplexVar v)) =
   2.753 +	     singleton v
   2.754 +	  |  collect_lower_bounded_vars 
   2.755 +		 (cplexBound (_, cplexLeq, cplexVar v)) = 
   2.756 +	     singleton v
   2.757 +	  |  collect_lower_bounded_vars 
   2.758 +		 (cplexBound (cplexVar v, cplexGe,_)) = 
   2.759 +	     singleton v
   2.760 +	  |  collect_lower_bounded_vars 
   2.761 +		 (cplexBound (cplexVar v, cplexGeq, _)) = 
   2.762 +	     singleton v
   2.763 +	  | collect_lower_bounded_vars
   2.764 +	    (cplexBound (cplexVar v, cplexEq, _)) = 
   2.765 +	    singleton v
   2.766 +	  |  collect_lower_bounded_vars _ = emptyset
   2.767 +	
   2.768 +	val lvars = mergemap collect_lower_bounded_vars bounds        
   2.769 +	val positive_vars = diff cvars lvars			   
   2.770 +	val zero = cplexNum "0"
   2.771 +	
   2.772 +	fun make_pos_constr v = 
   2.773 +	    (None, cplexConstr (cplexGeq, ((cplexVar v), zero)))
   2.774 +	
   2.775 +	fun make_free_bound v = 
   2.776 +	    cplexBounds (cplexNeg cplexInf, cplexLeq, 
   2.777 +			 cplexVar v, cplexLeq,
   2.778 +			 cplexInf)
   2.779 +	
   2.780 +	val pos_constrs = rev(Symtab.foldl 
   2.781 +			      (fn (l, (k,v)) => (make_pos_constr k) :: l)
   2.782 +			      ([], positive_vars))
   2.783 +        val bound_constrs = map (fn x => (None, x)) 
   2.784 +				(flat (map bound2constr bounds))
   2.785 +	val constraints' = constraints @ pos_constrs @ bound_constrs	   
   2.786 +	val bounds' = rev(Symtab.foldl (fn (l, (v,_)) => 
   2.787 +					   (make_free_bound v)::l)
   2.788 +				       ([], cvars))			  
   2.789 +    in
   2.790 +	cplexProg (name, goal, constraints', bounds')
   2.791 +    end
   2.792 +
   2.793 +fun relax_strict_ineqs (cplexProg (name, goals, constrs, bounds)) = 
   2.794 +    let
   2.795 +	fun relax cplexLe = cplexLeq
   2.796 +	  | relax cplexGe = cplexGeq
   2.797 +	  | relax x = x
   2.798 +	
   2.799 +	fun relax_constr (n, cplexConstr(c, (t1, t2))) = 
   2.800 +	    (n, cplexConstr(relax c, (t1, t2)))
   2.801 +	    
   2.802 +	fun relax_bounds (cplexBounds (t1, c1, t2, c2, t3)) = 
   2.803 +	    cplexBounds (t1, relax c1, t2, relax c2, t3)
   2.804 +	  | relax_bounds (cplexBound (t1, c, t2)) = 
   2.805 +	    cplexBound (t1, relax c, t2)
   2.806 +    in
   2.807 +	cplexProg (name, 
   2.808 +		   goals, 
   2.809 +		   map relax_constr constrs, 
   2.810 +		   map relax_bounds bounds)
   2.811 +    end
   2.812 +
   2.813 +datatype cplexResult = Unbounded 
   2.814 +		     | Infeasible 
   2.815 +		     | Undefined
   2.816 +		     | Optimal of string * ((string * string) list)
   2.817 +
   2.818 +fun is_separator x = forall (fn c => c = #"-") (String.explode x)
   2.819 +
   2.820 +fun is_sign x = (x = "+" orelse x = "-")
   2.821 +
   2.822 +fun is_colon x = (x = ":")
   2.823 +
   2.824 +fun is_resultsymbol a = 
   2.825 +    let
   2.826 +	val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~-"
   2.827 +	fun is_symbol_char c = Char.isAlphaNum c orelse 
   2.828 +			       exists (fn d => d=c) symbol_char
   2.829 +	fun is_symbol_start c = is_symbol_char c andalso 
   2.830 +				not (Char.isDigit c) andalso 
   2.831 +				not (c= #".") andalso
   2.832 +				not (c= #"-")
   2.833 +	val b = String.explode a
   2.834 +    in
   2.835 +	b <> [] andalso is_symbol_start (hd b) andalso 
   2.836 +	forall is_symbol_char b
   2.837 +    end
   2.838 +
   2.839 +val TOKEN_SIGN = 100
   2.840 +val TOKEN_COLON = 101
   2.841 +val TOKEN_SEPARATOR = 102
   2.842 +
   2.843 +fun load_cplexResult name = 
   2.844 +    let
   2.845 +	val flist = [(is_NL, TOKEN_NL), 
   2.846 +		     (is_blank, TOKEN_BLANK),
   2.847 +		     (is_num, TOKEN_NUM),
   2.848 +		     (is_sign, TOKEN_SIGN),
   2.849 +                     (is_colon, TOKEN_COLON),
   2.850 +		     (is_cmp, TOKEN_CMP),
   2.851 +		     (is_resultsymbol, TOKEN_SYMBOL),
   2.852 +		     (is_separator, TOKEN_SEPARATOR)]
   2.853 +
   2.854 +	val tokenize = tokenize_general flist
   2.855 +
   2.856 +	val f = TextIO.openIn name  
   2.857 +
   2.858 +	val rest = ref []
   2.859 +		   
   2.860 +	fun readToken_helper () =
   2.861 +	    if length (!rest) > 0 then  
   2.862 +		let val u = hd (!rest) in 
   2.863 +		    (
   2.864 +		     rest := tl (!rest); 
   2.865 +		     Some u
   2.866 +		    ) 
   2.867 +		end
   2.868 +	    else 
   2.869 +		let val s = TextIO.inputLine f in
   2.870 +		    if s = "" then None else (rest := tokenize s; readToken_helper())					     
   2.871 +		end           
   2.872 +		
   2.873 +	fun is_tt tok ty = (tok <> None andalso (fst (the tok)) = ty)
   2.874 +
   2.875 +	fun pushToken a = if a = None then () else (rest := ((the a)::(!rest)))
   2.876 +	
   2.877 +	fun readToken () =
   2.878 +	    let val t = readToken_helper () in
   2.879 +		if is_tt t TOKEN_BLANK then 
   2.880 +		    readToken ()
   2.881 +		else if is_tt t TOKEN_NL then 
   2.882 +		    let val t2 = readToken_helper () in
   2.883 +			if is_tt t2 TOKEN_SIGN then 
   2.884 +			    (pushToken (Some (TOKEN_SEPARATOR, "-")); t)
   2.885 +			else 
   2.886 +			    (pushToken t2; t)
   2.887 +		    end
   2.888 +		else if is_tt t TOKEN_SIGN then
   2.889 +		    let val t2 = readToken_helper () in
   2.890 +			if is_tt t2 TOKEN_NUM then
   2.891 +			    (Some (TOKEN_NUM, (snd (the t))^(snd (the t2))))
   2.892 +			else
   2.893 +			    (pushToken t2; t)
   2.894 +		    end
   2.895 +		else
   2.896 +		    t
   2.897 +	    end	    		
   2.898 +
   2.899 +        fun readRestOfLine P = 
   2.900 +	    let 
   2.901 +		val t = readToken ()
   2.902 +	    in
   2.903 +		if is_tt t TOKEN_NL orelse t = None 
   2.904 +		then P
   2.905 +		else readRestOfLine P
   2.906 +	    end			
   2.907 +							 
   2.908 +	fun readHeader () = 
   2.909 +	    let
   2.910 +		fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ())))
   2.911 +		fun readObjective () = readRestOfLine ("OBJECTIVE", snd (the (readToken (); readToken (); readToken ())))
   2.912 +		val t1 = readToken ()
   2.913 +		val t2 = readToken ()
   2.914 +	    in
   2.915 +		if is_tt t1 TOKEN_SYMBOL andalso is_tt t2 TOKEN_COLON 
   2.916 +		then
   2.917 +		    case to_upper (snd (the t1)) of
   2.918 +			"STATUS" => (readStatus ())::(readHeader ())
   2.919 +		      | "OBJECTIVE" => (readObjective())::(readHeader ())
   2.920 +		      | _ => (readRestOfLine (); readHeader ())
   2.921 +		else
   2.922 +		    (pushToken t2; pushToken t1; [])
   2.923 +	    end
   2.924 +
   2.925 +	fun skip_until_sep () = 
   2.926 +	    let val x = readToken () in
   2.927 +		if is_tt x TOKEN_SEPARATOR then
   2.928 +		    readRestOfLine ()
   2.929 +		else 
   2.930 +		    skip_until_sep ()
   2.931 +	    end
   2.932 +
   2.933 +	fun load_value () =
   2.934 +	    let 
   2.935 +		val t1 = readToken ()
   2.936 +		val t2 = readToken ()
   2.937 +	    in
   2.938 +		if is_tt t1 TOKEN_NUM andalso is_tt t2 TOKEN_SYMBOL then
   2.939 +		    let 
   2.940 +			val t = readToken ()
   2.941 +			val state = if is_tt t TOKEN_NL then readToken () else t
   2.942 +			val _ = if is_tt state TOKEN_SYMBOL then () else raise (Load_cplexResult "state expected")
   2.943 +			val k = readToken ()
   2.944 +		    in
   2.945 +			if is_tt k TOKEN_NUM then
   2.946 +			    readRestOfLine (Some (snd (the t2), snd (the k)))
   2.947 +			else
   2.948 +			    raise (Load_cplexResult "number expected")
   2.949 +		    end				
   2.950 +		else
   2.951 +		    (pushToken t2; pushToken t1; None)
   2.952 +	    end
   2.953 +
   2.954 +	fun load_values () = 
   2.955 +	    let val v = load_value () in
   2.956 +		if v = None then [] else (the v)::(load_values ())
   2.957 +	    end
   2.958 +	
   2.959 +	val header = readHeader () 
   2.960 +
   2.961 +	val result = 
   2.962 +	    case assoc (header, "STATUS") of
   2.963 +		Some "INFEASIBLE" => Infeasible
   2.964 +	      | Some "UNBOUNDED" => Unbounded
   2.965 +	      | Some "OPTIMAL" => Optimal (the (assoc (header, "OBJECTIVE")), 
   2.966 +					   (skip_until_sep (); 
   2.967 +					    skip_until_sep ();
   2.968 +					    load_values ()))
   2.969 +	      | _ => Undefined
   2.970 +
   2.971 +	val _ = TextIO.closeIn f
   2.972 +    in
   2.973 +	result
   2.974 +    end
   2.975 +    handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s))
   2.976 +	 | OPTION => raise (Load_cplexResult "OPTION")
   2.977 +	 | x => raise x
   2.978 +
   2.979 +exception Execute of string;
   2.980 +
   2.981 +fun execute_cplex lpname resultname = 
   2.982 +let
   2.983 +    fun wrap s = "\""^s^"\""
   2.984 +    val cplex_path = getenv "CPLEX"
   2.985 +    val cplex = if cplex_path = "" then "glpsol" else cplex_path
   2.986 +    val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname)
   2.987 +in
   2.988 +    execute command
   2.989 +end
   2.990 +
   2.991 +fun tmp_file s = Path.pack (Path.expand (File.tmp_path (Path.make [s])))
   2.992 +
   2.993 +fun solve prog =
   2.994 +    let 
   2.995 +	val name = LargeInt.toString (Time.toMicroseconds (Time.now ()))
   2.996 +	val lpname = tmp_file (name^".lp")
   2.997 +	val resultname = tmp_file (name^".result")
   2.998 +	val _ = save_cplexFile lpname prog
   2.999 +	val answer = execute_cplex lpname resultname
  2.1000 +    in
  2.1001 +	let
  2.1002 +	    val result = load_cplexResult resultname
  2.1003 +	    val _ = OS.FileSys.remove lpname
  2.1004 +	    val _ = OS.FileSys.remove resultname
  2.1005 +	in
  2.1006 +	    result
  2.1007 +	end 
  2.1008 +	handle (Load_cplexResult s) => raise (Execute ("Load_cplexResult: "^s^"\nExecute: "^answer))
  2.1009 +	     | _ => raise (Execute answer)
  2.1010 +    end			     
  2.1011 +end;
  2.1012 +
  2.1013 +val demofile = "/home/obua/flyspeck/kepler/LP/cplexPent2.lp45"
  2.1014 +val demoout = "/home/obua/flyspeck/kepler/LP/test.out"
  2.1015 +val demoresult = "/home/obua/flyspeck/kepler/LP/try/test2.sol"
  2.1016 +	   
  2.1017 +fun loadcplex () = Cplex.relax_strict_ineqs 
  2.1018 +		   (Cplex.load_cplexFile demofile)
  2.1019 +
  2.1020 +fun writecplex lp = Cplex.save_cplexFile demoout lp
  2.1021 +
  2.1022 +fun test () = 
  2.1023 +    let
  2.1024 +	val lp = loadcplex ()
  2.1025 +	val lp2 = Cplex.elim_nonfree_bounds lp
  2.1026 +    in
  2.1027 +	writecplex lp2
  2.1028 +    end
  2.1029 +
  2.1030 +fun loadresult () = Cplex.load_cplexResult demoresult;
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/src/HOL/Matrix/CplexMatrixConverter.ML	Fri Sep 03 17:10:36 2004 +0200
     3.3 @@ -0,0 +1,132 @@
     3.4 +signature MATRIX_BUILDER =
     3.5 +sig
     3.6 +    type vector
     3.7 +    type matrix
     3.8 +    
     3.9 +    val empty_vector : vector
    3.10 +    val empty_matrix : matrix
    3.11 +		       
    3.12 +    exception Nat_expected of int
    3.13 +    val set_elem : vector -> int -> string -> vector			
    3.14 +    val set_vector : matrix -> int -> vector -> matrix
    3.15 +end;
    3.16 +
    3.17 +signature CPLEX_MATRIX_CONVERTER = 
    3.18 +sig
    3.19 +    structure cplex : CPLEX
    3.20 +    structure matrix_builder : MATRIX_BUILDER 
    3.21 +    type vector = matrix_builder.vector
    3.22 +    type matrix = matrix_builder.matrix
    3.23 +    type naming = int * (int -> string) * (string -> int)
    3.24 +    
    3.25 +    exception Converter of string
    3.26 +
    3.27 +    (* program must fulfill is_normed_cplexProg and must be an element of the image of elim_nonfree_bounds *)
    3.28 +    (* convert_prog maximize c A b naming *)
    3.29 +    val convert_prog : cplex.cplexProg -> bool * vector * matrix * vector * naming
    3.30 +
    3.31 +    (* results must be optimal, converts_results returns the optimal value as string and the solution as vector *)
    3.32 +    (* convert_results results name2index *)
    3.33 +    val convert_results : cplex.cplexResult -> (string -> int) -> string * vector
    3.34 +end;
    3.35 +
    3.36 +functor MAKE_CPLEX_MATRIX_CONVERTER (structure cplex: CPLEX and matrix_builder: MATRIX_BUILDER) : CPLEX_MATRIX_CONVERTER =
    3.37 +struct
    3.38 +
    3.39 +structure cplex = cplex
    3.40 +structure matrix_builder = matrix_builder
    3.41 +type matrix = matrix_builder.matrix
    3.42 +type vector = matrix_builder.vector		
    3.43 +type naming = int * (int -> string) * (string -> int)
    3.44 +		     
    3.45 +open matrix_builder 
    3.46 +open cplex
    3.47 +     
    3.48 +exception Converter of string;
    3.49 +
    3.50 +structure Inttab = TableFun(type key = int val ord = int_ord);
    3.51 +
    3.52 +fun neg_term (cplexNeg t) = t
    3.53 +  | neg_term (cplexSum ts) = cplexSum (map neg_term ts)
    3.54 +  | neg_term t = cplexNeg t 
    3.55 +	      
    3.56 +fun convert_prog (cplexProg (s, goal, constrs, bounds)) = 
    3.57 +    let	
    3.58 +	fun build_naming index i2s s2i [] = (index, i2s, s2i)
    3.59 +	  | build_naming index i2s s2i (cplexBounds (cplexNeg cplexInf, cplexLeq, cplexVar v, cplexLeq, cplexInf)::bounds)
    3.60 +	    = build_naming (index+1) (Inttab.update ((index, v), i2s)) (Symtab.update_new ((v, index), s2i)) bounds
    3.61 +	  | build_naming _ _ _ _ = raise (Converter "nonfree bound")
    3.62 +
    3.63 +	val (varcount, i2s_tab, s2i_tab) = build_naming 0 Inttab.empty Symtab.empty bounds
    3.64 +
    3.65 +	fun i2s i = case Inttab.lookup (i2s_tab, i) of None => raise (Converter "index not found")
    3.66 +						     | Some n => n
    3.67 +	fun s2i s = case Symtab.lookup (s2i_tab, s) of None => raise (Converter ("name not found: "^s))
    3.68 +						     | Some i => i
    3.69 +	fun num2str positive (cplexNeg t) = num2str (not positive) t
    3.70 +	  | num2str positive (cplexNum num) = if positive then num else "-"^num			
    3.71 +	  | num2str _ _ = raise (Converter "term is not a (possibly signed) number")
    3.72 +
    3.73 +	fun setprod vec positive (cplexNeg t) = setprod vec (not positive) t  
    3.74 +	  | setprod vec positive (cplexVar v) = set_elem vec (s2i v) (if positive then "1" else "-1")
    3.75 +	  | setprod vec positive (cplexProd (cplexNum num, cplexVar v)) = 
    3.76 +	    set_elem vec (s2i v) (if positive then num else "-"^num)
    3.77 +	  | setprod _ _ _ = raise (Converter "term is not a normed product")	
    3.78 +
    3.79 +	fun sum2vec (cplexSum ts) = foldl (fn (vec, t) => setprod vec true t) (empty_vector, ts)
    3.80 +	  | sum2vec t = setprod empty_vector true t						
    3.81 +
    3.82 +	fun constrs2Ab j A b [] = (A, b)
    3.83 +	  | constrs2Ab j A b ((_, cplexConstr (cplexLeq, (t1,t2)))::cs) = 
    3.84 +	    constrs2Ab (j+1) (set_vector A j (sum2vec t1)) (set_elem b j (num2str true t2)) cs
    3.85 +	  | constrs2Ab j A b ((_, cplexConstr (cplexGeq, (t1,t2)))::cs) = 
    3.86 +	    constrs2Ab (j+1) (set_vector A j (sum2vec (neg_term t1))) (set_elem b j (num2str true (neg_term t2))) cs
    3.87 +	  | constrs2Ab j A b ((_, cplexConstr (cplexEq, (t1,t2)))::cs) =
    3.88 +	    constrs2Ab j A b ((None, cplexConstr (cplexLeq, (t1,t2)))::
    3.89 +			      (None, cplexConstr (cplexGeq, (t1, t2)))::cs)
    3.90 +	  | constrs2Ab _ _ _ _ = raise (Converter "no strict constraints allowed")
    3.91 +
    3.92 +	val (A, b) = constrs2Ab 0 empty_matrix empty_vector constrs
    3.93 +								 
    3.94 +	val (goal_maximize, goal_term) = 
    3.95 +	    case goal of
    3.96 +		(cplexMaximize t) => (true, t)
    3.97 +	      | (cplexMinimize t) => (false, t)				     
    3.98 +    in	  
    3.99 +	(goal_maximize, sum2vec goal_term, A, b, (varcount, i2s, s2i))
   3.100 +    end
   3.101 +
   3.102 +fun convert_results (cplex.Optimal (opt, entries)) name2index =
   3.103 +    let
   3.104 +	fun setv (v, (name, value)) = (matrix_builder.set_elem v (name2index name) value) 
   3.105 +    in
   3.106 +	(opt, foldl setv (matrix_builder.empty_vector, entries))
   3.107 +    end
   3.108 +  | convert_results _ _ = raise (Converter "No optimal result")
   3.109 +    
   3.110 +
   3.111 +end;
   3.112 +
   3.113 +structure SimpleMatrixBuilder : MATRIX_BUILDER = 
   3.114 +struct
   3.115 +type vector = (int * string) list
   3.116 +type matrix = (int * vector) list
   3.117 +
   3.118 +val empty_matrix = []
   3.119 +val empty_vector = []
   3.120 +
   3.121 +exception Nat_expected of int;
   3.122 +
   3.123 +fun set_elem v i s = v @ [(i, s)] 
   3.124 +
   3.125 +fun set_vector m i v = m @ [(i, v)]
   3.126 +    
   3.127 +end;	      
   3.128 +
   3.129 +
   3.130 +
   3.131 +structure SimpleCplexMatrixConverter = MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = SimpleMatrixBuilder);
   3.132 +
   3.133 +
   3.134 +
   3.135 +
     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/src/HOL/Matrix/ExactFloatingPoint.ML	Fri Sep 03 17:10:36 2004 +0200
     4.3 @@ -0,0 +1,215 @@
     4.4 +structure ExactFloatingPoint :
     4.5 +sig
     4.6 +    exception Destruct_floatstr of string
     4.7 +				   
     4.8 +    val destruct_floatstr : (char -> bool) -> (char -> bool) -> string -> bool * string * string * bool * string
     4.9 +									  
    4.10 +    exception Floating_point of string
    4.11 +				
    4.12 +    type floatrep = IntInf.int * IntInf.int
    4.13 +    val approx_dec_by_bin : IntInf.int -> floatrep -> floatrep * floatrep
    4.14 +    val approx_decstr_by_bin : int -> string -> floatrep * floatrep
    4.15 +end 
    4.16 +=
    4.17 +struct
    4.18 +
    4.19 +fun fst (a,b) = a
    4.20 +fun snd (a,b) = b
    4.21 +
    4.22 +val filter = List.filter;
    4.23 +
    4.24 +exception Destruct_floatstr of string;
    4.25 +
    4.26 +fun destruct_floatstr isDigit isExp number = 
    4.27 +    let
    4.28 +	val numlist = filter (not o Char.isSpace) (String.explode number)
    4.29 +	
    4.30 +	fun countsigns ((#"+")::cs) = countsigns cs
    4.31 +	  | countsigns ((#"-")::cs) = 
    4.32 +	    let	
    4.33 +		val (positive, rest) = countsigns cs 
    4.34 +	    in
    4.35 +		(not positive, rest)
    4.36 +	    end
    4.37 +	  | countsigns cs = (true, cs)
    4.38 +
    4.39 +	fun readdigits [] = ([], [])
    4.40 +	  | readdigits (q as c::cs) = 
    4.41 +	    if (isDigit c) then 
    4.42 +		let
    4.43 +		    val (digits, rest) = readdigits cs
    4.44 +		in
    4.45 +		    (c::digits, rest)
    4.46 +		end
    4.47 +	    else
    4.48 +		([], q)		
    4.49 +
    4.50 +	fun readfromexp_helper cs =
    4.51 +	    let
    4.52 +		val (positive, rest) = countsigns cs
    4.53 +		val (digits, rest') = readdigits rest
    4.54 +	    in
    4.55 +		case rest' of
    4.56 +		    [] => (positive, digits)
    4.57 +		  | _ => raise (Destruct_floatstr number)
    4.58 +	    end	    
    4.59 +
    4.60 +	fun readfromexp [] = (true, [])
    4.61 +	  | readfromexp (c::cs) = 
    4.62 +	    if isExp c then
    4.63 +		readfromexp_helper cs
    4.64 +	    else 
    4.65 +		raise (Destruct_floatstr number)		
    4.66 +
    4.67 +	fun readfromdot [] = ([], readfromexp [])
    4.68 +	  | readfromdot ((#".")::cs) = 
    4.69 +	    let		
    4.70 +		val (digits, rest) = readdigits cs
    4.71 +		val exp = readfromexp rest
    4.72 +	    in
    4.73 +		(digits, exp)
    4.74 +	    end		
    4.75 +	  | readfromdot cs = readfromdot ((#".")::cs)
    4.76 +			    
    4.77 +	val (positive, numlist) = countsigns numlist				 
    4.78 +	val (digits1, numlist) = readdigits numlist				 
    4.79 + 	val (digits2, exp) = readfromdot numlist
    4.80 +    in
    4.81 +	(positive, String.implode digits1, String.implode digits2, fst exp, String.implode (snd exp))
    4.82 +    end
    4.83 +
    4.84 +type floatrep = IntInf.int * IntInf.int
    4.85 +
    4.86 +exception Floating_point of string;
    4.87 +
    4.88 +val ln2_10 = (Math.ln 10.0)/(Math.ln 2.0)
    4.89 +	
    4.90 +fun intmul a b = IntInf.* (a,b)
    4.91 +fun intsub a b = IntInf.- (a,b)	
    4.92 +fun intadd a b = IntInf.+ (a,b) 		 
    4.93 +fun intpow a b = IntInf.pow (a, IntInf.toInt b);
    4.94 +fun intle a b = IntInf.<= (a, b);
    4.95 +fun intless a b = IntInf.< (a, b);
    4.96 +fun intneg a = IntInf.~ a;
    4.97 +val zero = IntInf.fromInt 0;
    4.98 +val one = IntInf.fromInt 1;
    4.99 +val two = IntInf.fromInt 2;
   4.100 +val ten = IntInf.fromInt 10;
   4.101 +val five = IntInf.fromInt 5;
   4.102 +
   4.103 +fun find_most_significant q r = 
   4.104 +    let 
   4.105 +	fun int2real i = 
   4.106 +	    case Real.fromString (IntInf.toString i) of 
   4.107 +		SOME r => r 
   4.108 +	      | NONE => raise (Floating_point "int2real")	
   4.109 +	fun subtract (q, r) (q', r') = 
   4.110 +	    if intle r r' then
   4.111 +		(intsub q (intmul q' (intpow ten (intsub r' r))), r)
   4.112 +	    else
   4.113 +		(intsub (intmul q (intpow ten (intsub r r'))) q', r')
   4.114 +	fun bin2dec d =
   4.115 +	    if intle zero d then 
   4.116 +		(intpow two d, zero)
   4.117 +	    else
   4.118 +		(intpow five (intneg d), d)				
   4.119 +		
   4.120 +	val L = IntInf.fromInt (Real.floor (int2real (IntInf.fromInt (IntInf.log2 q)) + (int2real r) * ln2_10))	
   4.121 +	val L1 = intadd L one
   4.122 +
   4.123 +	val (q1, r1) = subtract (q, r) (bin2dec L1) 		
   4.124 +    in
   4.125 +	if intle zero q1 then 
   4.126 +	    let
   4.127 +		val (q2, r2) = subtract (q, r) (bin2dec (intadd L1 one))
   4.128 +	    in
   4.129 +		if intle zero q2 then 
   4.130 +		    raise (Floating_point "find_most_significant")
   4.131 +		else
   4.132 +		    (L1, (q1, r1))
   4.133 +	    end
   4.134 +	else
   4.135 +	    let
   4.136 +		val (q0, r0) = subtract (q, r) (bin2dec L)
   4.137 +	    in
   4.138 +		if intle zero q0 then
   4.139 +		    (L, (q0, r0))
   4.140 +		else
   4.141 +		    raise (Floating_point "find_most_significant")
   4.142 +	    end		    
   4.143 +    end
   4.144 +
   4.145 +fun approx_dec_by_bin n (q,r) =
   4.146 +    let	
   4.147 +	fun addseq acc d' [] = acc
   4.148 +	  | addseq acc d' (d::ds) = addseq (intadd acc (intpow two (intsub d d'))) d' ds
   4.149 +
   4.150 +	fun seq2bin [] = (zero, zero)
   4.151 +	  | seq2bin (d::ds) = (intadd (addseq zero d ds) one, d)
   4.152 +
   4.153 +	fun approx d_seq d0 precision (q,r) = 
   4.154 +	    if q = zero then 
   4.155 +		let val x = seq2bin d_seq in
   4.156 +		    (x, x)
   4.157 +		end
   4.158 +	    else    
   4.159 +		let 
   4.160 +		    val (d, (q', r')) = find_most_significant q r
   4.161 +		in	
   4.162 +		    if intless precision (intsub d0 d) then 
   4.163 +			let 
   4.164 +			    val d' = intsub d0 precision
   4.165 +			    val x1 = seq2bin (d_seq)
   4.166 +			    val x2 = (intadd (intmul (fst x1) (intpow two (intsub (snd x1) d'))) one,  d') (* = seq2bin (d'::d_seq) *) 
   4.167 +			in
   4.168 +			    (x1, x2)
   4.169 +			end
   4.170 +		    else
   4.171 +			approx (d::d_seq) d0 precision (q', r') 						    		
   4.172 +		end		
   4.173 +	
   4.174 +	fun approx_start precision (q, r) =
   4.175 +	    if q = zero then 
   4.176 +		((zero, zero), (zero, zero))
   4.177 +	    else
   4.178 +		let 
   4.179 +		    val (d, (q', r')) = find_most_significant q r
   4.180 +		in	
   4.181 +		    if intle precision zero then 
   4.182 +			let
   4.183 +			    val x1 = seq2bin [d]
   4.184 +			in
   4.185 +			    if q' = zero then 
   4.186 +				(x1, x1)
   4.187 +			    else
   4.188 +				(x1, seq2bin [intadd d one])
   4.189 +			end
   4.190 +		    else
   4.191 +			approx [d] d precision (q', r')
   4.192 +		end		
   4.193 +    in
   4.194 +	if intle zero q then 
   4.195 +	    approx_start n (q,r)
   4.196 +	else
   4.197 +	    let 
   4.198 +		val ((a1,b1), (a2, b2)) = approx_start n (intneg q, r) 
   4.199 +	    in
   4.200 +		((intneg a2, b2), (intneg a1, b1))
   4.201 +	    end					
   4.202 +    end
   4.203 +
   4.204 +fun approx_decstr_by_bin n decstr =
   4.205 +    let 
   4.206 +	fun str2int s = case IntInf.fromString s of SOME x => x | NONE => zero 
   4.207 +	fun signint p x = if p then x else intneg x
   4.208 +
   4.209 +	val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr
   4.210 +	val s = IntInf.fromInt (size d2)
   4.211 +		
   4.212 +	val q = signint p (intadd (intmul (str2int d1) (intpow ten s)) (str2int d2))
   4.213 +	val r = intsub (signint ep (str2int e)) s
   4.214 +    in
   4.215 +	approx_dec_by_bin (IntInf.fromInt n) (q,r)
   4.216 +    end
   4.217 +
   4.218 +end;
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/src/HOL/Matrix/Float.thy	Fri Sep 03 17:10:36 2004 +0200
     5.3 @@ -0,0 +1,508 @@
     5.4 +theory Float = Hyperreal:
     5.5 +
     5.6 +constdefs  
     5.7 +  pow2 :: "int \<Rightarrow> real"
     5.8 +  "pow2 a == if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a))))" 
     5.9 +  float :: "int * int \<Rightarrow> real"
    5.10 +  "float x == (real (fst x)) * (pow2 (snd x))"
    5.11 +
    5.12 +lemma pow2_0[simp]: "pow2 0 = 1"
    5.13 +by (simp add: pow2_def)
    5.14 +
    5.15 +lemma pow2_1[simp]: "pow2 1 = 2"
    5.16 +by (simp add: pow2_def)
    5.17 +
    5.18 +lemma pow2_neg: "pow2 x = inverse (pow2 (-x))"
    5.19 +by (simp add: pow2_def)
    5.20 +
    5.21 +lemma pow2_add1: "pow2 (1 + a) = 2 * (pow2 a)" 
    5.22 +proof -
    5.23 +  have h: "! n. nat (2 + int n) - Suc 0 = nat (1 + int n)" by arith
    5.24 +  have g: "! a b. a - -1 = a + (1::int)" by arith
    5.25 +  have pos: "! n. pow2 (int n + 1) = 2 * pow2 (int n)"
    5.26 +    apply (auto, induct_tac n)
    5.27 +    apply (simp_all add: pow2_def)
    5.28 +    apply (rule_tac m1="2" and n1="nat (2 + int na)" in ssubst[OF realpow_num_eq_if])
    5.29 +    apply (auto simp add: h)
    5.30 +    apply arith
    5.31 +    done  
    5.32 +  show ?thesis
    5.33 +  proof (induct a)
    5.34 +    case (1 n)
    5.35 +    from pos show ?case by (simp add: ring_eq_simps)
    5.36 +  next
    5.37 +    case (2 n)
    5.38 +    show ?case
    5.39 +      apply (auto)
    5.40 +      apply (subst pow2_neg[of "- int n"])
    5.41 +      apply (subst pow2_neg[of "-1 - int n"])
    5.42 +      apply (auto simp add: g pos)
    5.43 +      done
    5.44 +  qed  
    5.45 +qed
    5.46 +  
    5.47 +lemma pow2_add: "pow2 (a+b) = (pow2 a) * (pow2 b)"
    5.48 +proof (induct b)
    5.49 +  case (1 n) 
    5.50 +  show ?case
    5.51 +  proof (induct n)
    5.52 +    case 0
    5.53 +    show ?case by simp
    5.54 +  next
    5.55 +    case (Suc m)
    5.56 +    show ?case by (auto simp add: ring_eq_simps pow2_add1 prems)
    5.57 +  qed
    5.58 +next
    5.59 +  case (2 n)
    5.60 +  show ?case 
    5.61 +  proof (induct n)
    5.62 +    case 0
    5.63 +    show ?case 
    5.64 +      apply (auto)
    5.65 +      apply (subst pow2_neg[of "a + -1"])
    5.66 +      apply (subst pow2_neg[of "-1"])
    5.67 +      apply (simp)
    5.68 +      apply (insert pow2_add1[of "-a"])
    5.69 +      apply (simp add: ring_eq_simps)
    5.70 +      apply (subst pow2_neg[of "-a"])
    5.71 +      apply (simp)
    5.72 +      done
    5.73 +    case (Suc m)
    5.74 +    have a: "int m - (a + -2) =  1 + (int m - a + 1)" by arith	
    5.75 +    have b: "int m - -2 = 1 + (int m + 1)" by arith
    5.76 +    show ?case
    5.77 +      apply (auto)
    5.78 +      apply (subst pow2_neg[of "a + (-2 - int m)"])
    5.79 +      apply (subst pow2_neg[of "-2 - int m"])
    5.80 +      apply (auto simp add: ring_eq_simps)
    5.81 +      apply (subst a)
    5.82 +      apply (subst b)
    5.83 +      apply (simp only: pow2_add1)
    5.84 +      apply (subst pow2_neg[of "int m - a + 1"])
    5.85 +      apply (subst pow2_neg[of "int m + 1"])
    5.86 +      apply auto
    5.87 +      apply (insert prems)
    5.88 +      apply (auto simp add: ring_eq_simps)
    5.89 +      done
    5.90 +  qed
    5.91 +qed
    5.92 +
    5.93 +lemma "float (a, e) + float (b, e) = float (a + b, e)"  
    5.94 +by (simp add: float_def ring_eq_simps)
    5.95 +
    5.96 +constdefs 
    5.97 +  int_of_real :: "real \<Rightarrow> int"
    5.98 +  "int_of_real x == SOME y. real y = x"  
    5.99 +  real_is_int :: "real \<Rightarrow> bool"
   5.100 +  "real_is_int x == ? (u::int). x = real u" 
   5.101 +
   5.102 +lemma real_is_int_def2: "real_is_int x = (x = real (int_of_real x))"
   5.103 +by (auto simp add: real_is_int_def int_of_real_def)
   5.104 +
   5.105 +lemma float_transfer: "real_is_int ((real a)*(pow2 c)) \<Longrightarrow> float (a, b) = float (int_of_real ((real a)*(pow2 c)), b - c)"
   5.106 +by (simp add: float_def real_is_int_def2 pow2_add[symmetric])
   5.107 +
   5.108 +lemma pow2_int: "pow2 (int c) = (2::real)^c"
   5.109 +by (simp add: pow2_def)
   5.110 +
   5.111 +lemma float_transfer_nat: "float (a, b) = float (a * 2^c, b - int c)" 
   5.112 +by (simp add: float_def pow2_int[symmetric] pow2_add[symmetric])
   5.113 +
   5.114 +lemma real_is_int_real[simp]: "real_is_int (real (x::int))"
   5.115 +by (auto simp add: real_is_int_def int_of_real_def)
   5.116 +
   5.117 +lemma int_of_real_real[simp]: "int_of_real (real x) = x"
   5.118 +by (simp add: int_of_real_def)
   5.119 +
   5.120 +lemma real_int_of_real[simp]: "real_is_int x \<Longrightarrow> real (int_of_real x) = x"
   5.121 +by (auto simp add: int_of_real_def real_is_int_def)
   5.122 +
   5.123 +lemma real_is_int_add_int_of_real: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> (int_of_real (a+b)) = (int_of_real a) + (int_of_real b)"
   5.124 +by (auto simp add: int_of_real_def real_is_int_def)
   5.125 +
   5.126 +lemma real_is_int_add[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a+b)"
   5.127 +apply (subst real_is_int_def2)
   5.128 +apply (simp add: real_is_int_add_int_of_real real_int_of_real)
   5.129 +done
   5.130 +
   5.131 +lemma int_of_real_sub: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> (int_of_real (a-b)) = (int_of_real a) - (int_of_real b)"
   5.132 +by (auto simp add: int_of_real_def real_is_int_def)
   5.133 +
   5.134 +lemma real_is_int_sub[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a-b)"
   5.135 +apply (subst real_is_int_def2)
   5.136 +apply (simp add: int_of_real_sub real_int_of_real)
   5.137 +done
   5.138 +
   5.139 +lemma real_is_int_rep: "real_is_int x \<Longrightarrow> ?! (a::int). real a = x"
   5.140 +by (auto simp add: real_is_int_def)
   5.141 +
   5.142 +lemma int_of_real_mult: 
   5.143 +  assumes "real_is_int a" "real_is_int b"
   5.144 +  shows "(int_of_real (a*b)) = (int_of_real a) * (int_of_real b)"
   5.145 +proof -
   5.146 +  from prems have a: "?! (a'::int). real a' = a" by (rule_tac real_is_int_rep, auto)
   5.147 +  from prems have b: "?! (b'::int). real b' = b" by (rule_tac real_is_int_rep, auto)
   5.148 +  from a obtain a'::int where a':"a = real a'" by auto
   5.149 +  from b obtain b'::int where b':"b = real b'" by auto
   5.150 +  have r: "real a' * real b' = real (a' * b')" by auto
   5.151 +  show ?thesis
   5.152 +    apply (simp add: a' b')
   5.153 +    apply (subst r)
   5.154 +    apply (simp only: int_of_real_real)
   5.155 +    done
   5.156 +qed
   5.157 +
   5.158 +lemma real_is_int_mult[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a*b)"
   5.159 +apply (subst real_is_int_def2)
   5.160 +apply (simp add: int_of_real_mult)
   5.161 +done
   5.162 +
   5.163 +lemma real_is_int_0[simp]: "real_is_int (0::real)"
   5.164 +by (simp add: real_is_int_def int_of_real_def)
   5.165 +
   5.166 +lemma real_is_int_1[simp]: "real_is_int (1::real)"
   5.167 +proof -
   5.168 +  have "real_is_int (1::real) = real_is_int(real (1::int))" by auto
   5.169 +  also have "\<dots> = True" by (simp only: real_is_int_real)
   5.170 +  ultimately show ?thesis by auto
   5.171 +qed
   5.172 +
   5.173 +lemma real_is_int_n1: "real_is_int (-1::real)"
   5.174 +proof -
   5.175 +  have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto
   5.176 +  also have "\<dots> = True" by (simp only: real_is_int_real)
   5.177 +  ultimately show ?thesis by auto
   5.178 +qed
   5.179 +
   5.180 +lemma real_is_int_number_of[simp]: "real_is_int ((number_of::bin\<Rightarrow>real) x)"
   5.181 +proof -
   5.182 +  have neg1: "real_is_int (-1::real)"
   5.183 +  proof -
   5.184 +    have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto
   5.185 +    also have "\<dots> = True" by (simp only: real_is_int_real)
   5.186 +    ultimately show ?thesis by auto
   5.187 +  qed
   5.188 +  
   5.189 +  { 
   5.190 +    fix x::int
   5.191 +    have "!! y. real_is_int ((number_of::bin\<Rightarrow>real) (Abs_Bin x))"
   5.192 +      apply (simp add: number_of_eq)
   5.193 +      apply (subst Abs_Bin_inverse)
   5.194 +      apply (simp add: Bin_def)
   5.195 +      apply (induct x)
   5.196 +      apply (induct_tac n)
   5.197 +      apply (simp)
   5.198 +      apply (simp)
   5.199 +      apply (induct_tac n)
   5.200 +      apply (simp add: neg1)
   5.201 +    proof -
   5.202 +      fix n :: nat
   5.203 +      assume rn: "(real_is_int (of_int (- (int (Suc n)))))"
   5.204 +      have s: "-(int (Suc (Suc n))) = -1 + - (int (Suc n))" by simp
   5.205 +      show "real_is_int (of_int (- (int (Suc (Suc n)))))"
   5.206 +	apply (simp only: s of_int_add)
   5.207 +	apply (rule real_is_int_add)
   5.208 +	apply (simp add: neg1)
   5.209 +	apply (simp only: rn)
   5.210 +	done
   5.211 +    qed
   5.212 +  }
   5.213 +  note Abs_Bin = this
   5.214 +  {
   5.215 +    fix x :: bin
   5.216 +    have "? u. x = Abs_Bin u"
   5.217 +      apply (rule exI[where x = "Rep_Bin x"])
   5.218 +      apply (simp add: Rep_Bin_inverse)
   5.219 +      done
   5.220 +  }
   5.221 +  then obtain u::int where "x = Abs_Bin u" by auto
   5.222 +  with Abs_Bin show ?thesis by auto
   5.223 +qed
   5.224 +
   5.225 +lemma int_of_real_0[simp]: "int_of_real (0::real) = (0::int)"
   5.226 +by (simp add: int_of_real_def)
   5.227 +
   5.228 +lemma int_of_real_1[simp]: "int_of_real (1::real) = (1::int)"
   5.229 +proof - 
   5.230 +  have 1: "(1::real) = real (1::int)" by auto
   5.231 +  show ?thesis by (simp only: 1 int_of_real_real)
   5.232 +qed
   5.233 +
   5.234 +lemma int_of_real_number_of[simp]: "int_of_real (number_of b) = number_of b"
   5.235 +proof -
   5.236 +  have "real_is_int (number_of b)" by simp
   5.237 +  then have uu: "?! u::int. number_of b = real u" by (auto simp add: real_is_int_rep)
   5.238 +  then obtain u::int where u:"number_of b = real u" by auto
   5.239 +  have "number_of b = real ((number_of b)::int)" 
   5.240 +    by (simp add: number_of_eq real_of_int_def)
   5.241 +  have ub: "number_of b = real ((number_of b)::int)" 
   5.242 +    by (simp add: number_of_eq real_of_int_def)
   5.243 +  from uu u ub have unb: "u = number_of b"
   5.244 +    by blast
   5.245 +  have "int_of_real (number_of b) = u" by (simp add: u)
   5.246 +  with unb show ?thesis by simp
   5.247 +qed
   5.248 +
   5.249 +lemma float_transfer_even: "even a \<Longrightarrow> float (a, b) = float (a div 2, b+1)"
   5.250 +  apply (subst float_transfer[where a="a" and b="b" and c="-1", simplified])
   5.251 +  apply (simp_all add: pow2_def even_def real_is_int_def ring_eq_simps)
   5.252 +  apply (auto)
   5.253 +proof -
   5.254 +  fix q::int
   5.255 +  have a:"b - (-1\<Colon>int) = (1\<Colon>int) + b" by arith
   5.256 +  show "(float (q, (b - (-1\<Colon>int)))) = (float (q, ((1\<Colon>int) + b)))" 
   5.257 +    by (simp add: a)
   5.258 +qed
   5.259 +    
   5.260 +consts
   5.261 +  norm_float :: "int*int \<Rightarrow> int*int"
   5.262 +
   5.263 +lemma int_div_zdiv: "int (a div b) = (int a) div (int b)"
   5.264 +apply (subst split_div, auto)
   5.265 +apply (subst split_zdiv, auto)
   5.266 +apply (rule_tac a="int (b * i) + int j" and b="int b" and r="int j" and r'=ja in IntDiv.unique_quotient)
   5.267 +apply (auto simp add: IntDiv.quorem_def int_eq_of_nat)
   5.268 +done
   5.269 +
   5.270 +lemma int_mod_zmod: "int (a mod b) = (int a) mod (int b)"
   5.271 +apply (subst split_mod, auto)
   5.272 +apply (subst split_zmod, auto)
   5.273 +apply (rule_tac a="int (b * i) + int j" and b="int b" and q="int i" and q'=ia in IntDiv.unique_remainder)
   5.274 +apply (auto simp add: IntDiv.quorem_def int_eq_of_nat)
   5.275 +done
   5.276 +
   5.277 +lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> abs((a::int) div 2) < abs a"
   5.278 +by arith
   5.279 +
   5.280 +lemma terminating_norm_float: "\<forall>a. (a::int) \<noteq> 0 \<and> even a \<longrightarrow> a \<noteq> 0 \<and> \<bar>a div 2\<bar> < \<bar>a\<bar>"
   5.281 +apply (auto)
   5.282 +apply (rule abs_div_2_less)
   5.283 +apply (auto)
   5.284 +done
   5.285 +
   5.286 +ML {* simp_depth_limit := 2 *} 
   5.287 +recdef norm_float "measure (% (a,b). nat (abs a))"
   5.288 +  "norm_float (a,b) = (if (a \<noteq> 0) & (even a) then norm_float (a div 2, b+1) else (if a=0 then (0,0) else (a,b)))"
   5.289 +(hints simp: terminating_norm_float)
   5.290 +ML {* simp_depth_limit := 1000 *}
   5.291 +
   5.292 +
   5.293 +lemma norm_float: "float x = float (norm_float x)"
   5.294 +proof -
   5.295 +  {
   5.296 +    fix a b :: int 
   5.297 +    have norm_float_pair: "float (a,b) = float (norm_float (a,b))" 
   5.298 +    proof (induct a b rule: norm_float.induct)
   5.299 +      case (1 u v)
   5.300 +      show ?case 
   5.301 +      proof cases
   5.302 +	assume u: "u \<noteq> 0 \<and> even u"
   5.303 +	with prems have ind: "float (u div 2, v + 1) = float (norm_float (u div 2, v + 1))" by auto
   5.304 +	with u have "float (u,v) = float (u div 2, v+1)" by (simp add: float_transfer_even) 
   5.305 +	then show ?thesis
   5.306 +	  apply (subst norm_float.simps)
   5.307 +	  apply (simp add: ind)
   5.308 +	  done
   5.309 +      next
   5.310 +	assume "~(u \<noteq> 0 \<and> even u)"
   5.311 +	then show ?thesis
   5.312 +	  by (simp add: prems float_def)
   5.313 +      qed
   5.314 +    qed
   5.315 +  }
   5.316 +  note helper = this
   5.317 +  have "? a b. x = (a,b)" by auto
   5.318 +  then obtain a b where "x = (a, b)" by blast
   5.319 +  then show ?thesis by (simp only: helper)
   5.320 +qed
   5.321 +
   5.322 +lemma pow2_int: "pow2 (int n) = 2^n"
   5.323 +  by (simp add: pow2_def)
   5.324 +
   5.325 +lemma float_add: 
   5.326 +  "float (a1, e1) + float (a2, e2) = 
   5.327 +  (if e1<=e2 then float (a1+a2*2^(nat(e2-e1)), e1) 
   5.328 +  else float (a1*2^(nat (e1-e2))+a2, e2))"
   5.329 +  apply (simp add: float_def ring_eq_simps)
   5.330 +  apply (auto simp add: pow2_int[symmetric] pow2_add[symmetric])
   5.331 +  done
   5.332 +
   5.333 +lemma float_mult:
   5.334 +  "float (a1, e1) * float (a2, e2) = 
   5.335 +  (float (a1 * a2, e1 + e2))"
   5.336 +  by (simp add: float_def pow2_add)
   5.337 +
   5.338 +lemma float_minus:
   5.339 +  "- (float (a,b)) = float (-a, b)"
   5.340 +  by (simp add: float_def)
   5.341 +
   5.342 +lemma zero_less_pow2:
   5.343 +  "0 < pow2 x"
   5.344 +proof -
   5.345 +  {
   5.346 +    fix y
   5.347 +    have "0 <= y \<Longrightarrow> 0 < pow2 y"    
   5.348 +      by (induct y, induct_tac n, simp_all add: pow2_add)
   5.349 +  }
   5.350 +  note helper=this
   5.351 +  show ?thesis
   5.352 +    apply (case_tac "0 <= x")
   5.353 +    apply (simp add: helper)
   5.354 +    apply (subst pow2_neg)
   5.355 +    apply (simp add: helper)
   5.356 +    done
   5.357 +qed
   5.358 +
   5.359 +lemma zero_le_float:
   5.360 +  "(0 <= float (a,b)) = (0 <= a)"
   5.361 +  apply (auto simp add: float_def)
   5.362 +  apply (auto simp add: zero_le_mult_iff zero_less_pow2) 
   5.363 +  apply (insert zero_less_pow2[of b])
   5.364 +  apply (simp_all)
   5.365 +  done
   5.366 +
   5.367 +lemma float_abs:
   5.368 +  "abs (float (a,b)) = (if 0 <= a then (float (a,b)) else (float (-a,b)))"
   5.369 +  apply (auto simp add: abs_if)
   5.370 +  apply (simp_all add: zero_le_float[symmetric, of a b] float_minus)
   5.371 +  done
   5.372 +
   5.373 +lemma norm_0_1: "(0::_::number_ring) = Numeral0 & (1::_::number_ring) = Numeral1"
   5.374 +  by auto
   5.375 +  
   5.376 +lemma add_left_zero: "0 + a = (a::'a::comm_monoid_add)"
   5.377 +  by simp
   5.378 +
   5.379 +lemma add_right_zero: "a + 0 = (a::'a::comm_monoid_add)"
   5.380 +  by simp
   5.381 +
   5.382 +lemma mult_left_one: "1 * a = (a::'a::semiring_1)"
   5.383 +  by simp
   5.384 +
   5.385 +lemma mult_right_one: "a * 1 = (a::'a::semiring_1)"
   5.386 +  by simp
   5.387 +
   5.388 +lemma int_pow_0: "(a::int)^(Numeral0) = 1"
   5.389 +  by simp
   5.390 +
   5.391 +lemma int_pow_1: "(a::int)^(Numeral1) = a"
   5.392 +  by simp
   5.393 +
   5.394 +lemma zero_eq_Numeral0_nring: "(0::'a::number_ring) = Numeral0"
   5.395 +  by simp
   5.396 +
   5.397 +lemma one_eq_Numeral1_nring: "(1::'a::number_ring) = Numeral1"
   5.398 +  by simp
   5.399 +
   5.400 +lemma zero_eq_Numeral0_nat: "(0::nat) = Numeral0"
   5.401 +  by simp
   5.402 +
   5.403 +lemma one_eq_Numeral1_nat: "(1::nat) = Numeral1"
   5.404 +  by simp
   5.405 +
   5.406 +lemma zpower_Pls: "(z::int)^Numeral0 = Numeral1"
   5.407 +  by simp
   5.408 +
   5.409 +lemma zpower_Min: "(z::int)^((-1)::nat) = Numeral1"
   5.410 +proof -
   5.411 +  have 1:"((-1)::nat) = 0"
   5.412 +    by simp
   5.413 +  show ?thesis by (simp add: 1)
   5.414 +qed
   5.415 +
   5.416 +lemma fst_cong: "a=a' \<Longrightarrow> fst (a,b) = fst (a',b)"
   5.417 +  by simp
   5.418 +
   5.419 +lemma snd_cong: "b=b' \<Longrightarrow> snd (a,b) = snd (a,b')"
   5.420 +  by simp
   5.421 +
   5.422 +lemma lift_bool: "x \<Longrightarrow> x=True"
   5.423 +  by simp
   5.424 +
   5.425 +lemma nlift_bool: "~x \<Longrightarrow> x=False"
   5.426 +  by simp
   5.427 +
   5.428 +lemma not_false_eq_true: "(~ False) = True" by simp
   5.429 +
   5.430 +lemma not_true_eq_false: "(~ True) = False" by simp
   5.431 +
   5.432 +
   5.433 +lemmas binarith = 
   5.434 +  Pls_0_eq Min_1_eq
   5.435 +  bin_pred_Pls bin_pred_Min bin_pred_1 bin_pred_0     
   5.436 +  bin_succ_Pls bin_succ_Min bin_succ_1 bin_succ_0
   5.437 +  bin_add_Pls bin_add_Min bin_add_BIT_0 bin_add_BIT_10
   5.438 +  bin_add_BIT_11 bin_minus_Pls bin_minus_Min bin_minus_1 
   5.439 +  bin_minus_0 bin_mult_Pls bin_mult_Min bin_mult_1 bin_mult_0 
   5.440 +  bin_add_Pls_right bin_add_Min_right
   5.441 +
   5.442 +thm eq_number_of_eq
   5.443 +
   5.444 +lemma int_eq_number_of_eq: "(((number_of v)::int)=(number_of w)) = iszero ((number_of (bin_add v (bin_minus w)))::int)"
   5.445 +  by simp
   5.446 +
   5.447 +lemma int_iszero_number_of_Pls: "iszero (Numeral0::int)" 
   5.448 +  by (simp only: iszero_number_of_Pls)
   5.449 +
   5.450 +lemma int_nonzero_number_of_Min: "~(iszero ((-1)::int))"
   5.451 +  by simp
   5.452 +thm iszero_number_of_1
   5.453 +
   5.454 +lemma int_iszero_number_of_0: "iszero ((number_of (w BIT False))::int) = iszero ((number_of w)::int)"
   5.455 +  by simp
   5.456 +
   5.457 +lemma int_iszero_number_of_1: "\<not> iszero ((number_of (w BIT True))::int)" 
   5.458 +  by simp
   5.459 +
   5.460 +lemma int_less_number_of_eq_neg: "(((number_of x)::int) < number_of y) = neg ((number_of (bin_add x (bin_minus y)))::int)"
   5.461 +  by simp
   5.462 +
   5.463 +lemma int_not_neg_number_of_Pls: "\<not> (neg (Numeral0::int))" 
   5.464 +  by simp
   5.465 +
   5.466 +lemma int_neg_number_of_Min: "neg (-1::int)"
   5.467 +  by simp
   5.468 +
   5.469 +lemma int_neg_number_of_BIT: "neg ((number_of (w BIT x))::int) = neg ((number_of w)::int)"
   5.470 +  by simp
   5.471 +
   5.472 +lemma int_le_number_of_eq: "(((number_of x)::int) \<le> number_of y) = (\<not> neg ((number_of (bin_add y (bin_minus x)))::int))"
   5.473 +  by simp
   5.474 +
   5.475 +lemmas intarithrel = 
   5.476 +  (*abs_zero abs_one*)
   5.477 +  int_eq_number_of_eq 
   5.478 +  lift_bool[OF int_iszero_number_of_Pls] nlift_bool[OF int_nonzero_number_of_Min] int_iszero_number_of_0 
   5.479 +  lift_bool[OF int_iszero_number_of_1] int_less_number_of_eq_neg nlift_bool[OF int_not_neg_number_of_Pls] lift_bool[OF int_neg_number_of_Min]
   5.480 +  int_neg_number_of_BIT int_le_number_of_eq
   5.481 +
   5.482 +lemma int_number_of_add_sym: "((number_of v)::int) + number_of w = number_of (bin_add v w)"
   5.483 +  by simp
   5.484 +
   5.485 +lemma int_number_of_diff_sym: "((number_of v)::int) - number_of w = number_of (bin_add v (bin_minus w))"
   5.486 +  by simp
   5.487 +
   5.488 +lemma int_number_of_mult_sym: "((number_of v)::int) * number_of w = number_of (bin_mult v w)"
   5.489 +  by simp
   5.490 +
   5.491 +lemma int_number_of_minus_sym: "- ((number_of v)::int) = number_of (bin_minus v)"
   5.492 +  by simp
   5.493 +
   5.494 +lemmas intarith = int_number_of_add_sym int_number_of_minus_sym int_number_of_diff_sym int_number_of_mult_sym
   5.495 +
   5.496 +lemmas natarith = (*zero_eq_Numeral0_nat one_eq_Numeral1_nat*) add_nat_number_of 
   5.497 +  diff_nat_number_of mult_nat_number_of eq_nat_number_of less_nat_number_of
   5.498 +
   5.499 +lemmas powerarith = nat_number_of zpower_number_of_even 
   5.500 +  zpower_number_of_odd[simplified zero_eq_Numeral0_nring one_eq_Numeral1_nring]   
   5.501 +  zpower_Pls zpower_Min
   5.502 +
   5.503 +lemmas floatarith[simplified norm_0_1] = float_add float_mult float_minus float_abs zero_le_float
   5.504 +
   5.505 +lemmas arith = binarith intarith intarithrel natarith powerarith floatarith not_false_eq_true not_true_eq_false
   5.506 +
   5.507 +(* needed for the verifying code generator *)
   5.508 +lemmas arith_no_let = arith[simplified Let_def]
   5.509 +
   5.510 +end
   5.511 + 
     6.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     6.2 +++ b/src/HOL/Matrix/FloatArith.ML	Fri Sep 03 17:10:36 2004 +0200
     6.3 @@ -0,0 +1,56 @@
     6.4 +structure FloatArith = 
     6.5 +struct
     6.6 +
     6.7 +type float = IntInf.int * IntInf.int 
     6.8 +
     6.9 +val izero = IntInf.fromInt 0
    6.10 +fun imul a b = IntInf.* (a,b)
    6.11 +fun isub a b = IntInf.- (a,b)
    6.12 +fun iadd a b = IntInf.+ (a,b)
    6.13 +
    6.14 +val floatzero = (izero, izero)
    6.15 +
    6.16 +fun positive_part (a,b) = 
    6.17 +    (if IntInf.< (a,izero) then izero else a, b)
    6.18 +
    6.19 +fun negative_part (a,b) = 
    6.20 +    (if IntInf.< (a,izero) then a else izero, b)
    6.21 +
    6.22 +fun is_negative (a,b) = 
    6.23 +    if IntInf.< (a, izero) then true else false
    6.24 +
    6.25 +fun is_positive (a,b) = 
    6.26 +    if IntInf.< (izero, a) then true else false
    6.27 +
    6.28 +fun is_zero (a,b) = 
    6.29 +    if a = izero then true else false
    6.30 +
    6.31 +fun ipow2 a = IntInf.pow ((IntInf.fromInt 2), IntInf.toInt a)
    6.32 +
    6.33 +fun add (a1, b1) (a2, b2) = 
    6.34 +    if b1 < b2 then
    6.35 +	(iadd a1 (imul a2 (ipow2 (isub b2 b1))), b1)
    6.36 +    else
    6.37 +	(iadd (imul a1 (ipow2 (isub b1 b2))) a2, b2)
    6.38 +
    6.39 +fun sub (a1, b1) (a2, b2) = 
    6.40 +    if b1 < b2 then
    6.41 +	(isub a1 (imul a2 (ipow2 (isub b2 b1))), b1)
    6.42 +    else
    6.43 +	(isub (imul a1 (ipow2 (isub b1 b2))) a2, b2)
    6.44 +
    6.45 +fun neg (a, b) = (IntInf.~ a, b)
    6.46 +
    6.47 +fun is_equal a b = is_zero (sub a b)
    6.48 +
    6.49 +fun is_less a b = is_negative (sub a b)
    6.50 +
    6.51 +fun max a b = if is_less a b then b else a
    6.52 +
    6.53 +fun min a b = if is_less a b then a else b
    6.54 +
    6.55 +fun abs a = if is_negative a then neg a else a
    6.56 +
    6.57 +fun mul (a1, b1) (a2, b2) = (imul a1 a2, iadd b1 b2)
    6.58 +
    6.59 +end;
     7.1 --- a/src/HOL/Matrix/LinProg.thy	Fri Sep 03 10:27:05 2004 +0200
     7.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     7.3 @@ -1,79 +0,0 @@
     7.4 -(*  Title:      HOL/Matrix/LinProg.thy
     7.5 -    ID:         $Id$
     7.6 -    Author:     Steven Obua
     7.7 -*)
     7.8 -
     7.9 -theory LinProg = Matrix:
    7.10 -
    7.11 -lemma linprog_by_duality_approx_general:
    7.12 -  assumes
    7.13 -  fmuladdprops:
    7.14 -  "! a b c d. a <= b & c <= d \<longrightarrow> fadd a c <= fadd b d"
    7.15 -  "! c a b. 0 <= c & a <= b \<longrightarrow> fmul c a <= fmul c b"
    7.16 -  "! a. fmul 0 a = 0"
    7.17 -  "! a. fmul a 0 = 0"
    7.18 -  "fadd 0 0 = 0"
    7.19 -  "associative fadd"
    7.20 -  "commutative fadd"
    7.21 -  "associative fmul"
    7.22 -  "distributive fmul fadd"
    7.23 -  and specificprops:
    7.24 -  "mult_matrix fmul fadd (combine_matrix fadd A dA) x <= (b::('a::{order, zero}) matrix)"
    7.25 -  "mult_matrix fmul fadd y A = c"
    7.26 -  "0 <= y"
    7.27 -  shows
    7.28 -  "combine_matrix fadd (mult_matrix fmul fadd c x) (mult_matrix fmul fadd (mult_matrix fmul fadd y dA) x)
    7.29 -  <= mult_matrix fmul fadd y b"
    7.30 -proof -
    7.31 -  let ?mul = "mult_matrix fmul fadd"
    7.32 -  let ?add = "combine_matrix fadd"
    7.33 -  let ?t1 = "?mul y (?mul (?add A dA) x)"
    7.34 -  have a: "?t1 <= ?mul y b" by (rule le_left_mult, simp_all!)
    7.35 -  have b: "?t1 = ?mul (?mul y (?add A dA)) x" by (simp! add: mult_matrix_assoc_simple[THEN sym])
    7.36 -  have assoc: "associative ?add" by (simp! add: combine_matrix_assoc)
    7.37 -  have r_distr: "r_distributive ?mul ?add"
    7.38 -    apply (rule r_distributive_matrix)
    7.39 -    by (simp! add: distributive_def)+
    7.40 -  have l_distr: "l_distributive ?mul ?add"
    7.41 -    apply (rule l_distributive_matrix)
    7.42 -    by (simp! add: distributive_def)+
    7.43 -  have c:"?mul y (?add A dA) = ?add (?mul y A) (?mul y dA)"
    7.44 -    by (insert r_distr, simp add: r_distributive_def)
    7.45 -  have d:"?mul (?add (?mul y A) (?mul y dA)) x = ?add (?mul (?mul y A) x) (?mul (?mul y dA) x)"
    7.46 -    by (insert l_distr, simp add: l_distributive_def)
    7.47 -  have e:"\<dots> = ?add (?mul c x) (?mul (?mul y dA) x)" by (simp!)
    7.48 -  from a b c d e show "?add (?mul c x) (?mul (?mul y dA) x) <= ?mul y b" by simp
    7.49 -qed
    7.50 -
    7.51 -lemma linprog_by_duality_approx:
    7.52 -  assumes
    7.53 -  "(A + dA) * x <= (b::('a::lordered_ring) matrix)"
    7.54 -  "y * A = c"
    7.55 -  "0 <= y"
    7.56 -  shows
    7.57 -  "c * x  + (y * dA) * x <= y * b"
    7.58 -apply (simp add: times_matrix_def plus_matrix_def)
    7.59 -apply (rule linprog_by_duality_approx_general)
    7.60 -apply (simp_all)
    7.61 -apply (simp_all add: associative_def add_assoc mult_assoc)
    7.62 -apply (simp_all add: commutative_def add_commute)
    7.63 -apply (simp_all add: distributive_def l_distributive_def r_distributive_def left_distrib right_distrib)
    7.64 -apply (simp_all! add: plus_matrix_def times_matrix_def)
    7.65 -apply (simp add: add_mono)
    7.66 -apply (simp add: mult_left_mono)
    7.67 -done
    7.68 -
    7.69 -lemma linprog_by_duality:
    7.70 -  assumes
    7.71 -  "A * x <= (b::('a::lordered_ring) matrix)"
    7.72 -  "y * A = c"
    7.73 -  "0 <= y"
    7.74 -  shows
    7.75 -  "c * x  <= y * b"
    7.76 -proof -
    7.77 -  have a:"(A + 0) * x <= b" by (simp!)
    7.78 -  have b:"c * x + (y*0)*x <= y * b" by (rule linprog_by_duality_approx, simp_all!)
    7.79 -  show "c * x <= y*b" by (insert b, simp)
    7.80 -qed
    7.81 -
    7.82 -end
     8.1 --- a/src/HOL/Matrix/Matrix.thy	Fri Sep 03 10:27:05 2004 +0200
     8.2 +++ b/src/HOL/Matrix/Matrix.thy	Fri Sep 03 17:10:36 2004 +0200
     8.3 @@ -21,7 +21,10 @@
     8.4    times_matrix_def: "A * B == mult_matrix (op *) (op +) A B"
     8.5  
     8.6  lemma is_meet_combine_matrix_meet: "is_meet (combine_matrix meet)"
     8.7 -by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le)
     8.8 +  by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le)
     8.9 +
    8.10 +lemma is_join_combine_matrix_join: "is_join (combine_matrix join)"
    8.11 +  by (simp_all add: is_join_def le_matrix_def join_left_le join_right_le join_imp_le)
    8.12  
    8.13  instance matrix :: (lordered_ab_group) lordered_ab_group_meet
    8.14  proof 
    8.15 @@ -284,7 +287,20 @@
    8.16    apply (auto)
    8.17    done
    8.18  
    8.19 +lemma Rep_minus[simp]: "Rep_matrix (-(A::_::lordered_ab_group)) x y = - (Rep_matrix A x y)"
    8.20 +  by (simp add: minus_matrix_def)
    8.21  
    8.22 +lemma join_matrix: "join (A::('a::lordered_ring) matrix) B = combine_matrix join A B"  
    8.23 +  apply (insert join_unique[of "(combine_matrix join)::('a matrix \<Rightarrow> 'a matrix \<Rightarrow> 'a matrix)", simplified is_join_combine_matrix_join])
    8.24 +  apply (simp)
    8.25 +  done
    8.26  
    8.27 +lemma meet_matrix: "meet (A::('a::lordered_ring) matrix) B = combine_matrix meet A B"
    8.28 +  apply (insert meet_unique[of "(combine_matrix meet)::('a matrix \<Rightarrow> 'a matrix \<Rightarrow> 'a matrix)", simplified is_meet_combine_matrix_meet])
    8.29 +  apply (simp)
    8.30 +  done
    8.31 +
    8.32 +lemma Rep_abs[simp]: "Rep_matrix (abs (A::_::lordered_ring)) x y = abs (Rep_matrix A x y)"
    8.33 +  by (simp add: abs_lattice join_matrix)
    8.34  
    8.35  end
     9.1 --- a/src/HOL/Matrix/MatrixGeneral.thy	Fri Sep 03 10:27:05 2004 +0200
     9.2 +++ b/src/HOL/Matrix/MatrixGeneral.thy	Fri Sep 03 17:10:36 2004 +0200
     9.3 @@ -1392,7 +1392,44 @@
     9.4    apply (rule le_foldseq)
     9.5    by (auto)
     9.6  
     9.7 +lemma spec2: "! j i. P j i \<Longrightarrow> P j i" by blast
     9.8 +lemma neg_imp: "(\<not> Q \<longrightarrow> \<not> P) \<Longrightarrow> P \<longrightarrow> Q" by blast
     9.9 +
    9.10  lemma singleton_matrix_le[simp]: "(singleton_matrix j i a <= singleton_matrix j i b) = (a <= (b::_::order))"
    9.11    by (auto simp add: le_matrix_def)
    9.12  
    9.13 +lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))"
    9.14 +  apply (auto)
    9.15 +  apply (simp add: le_matrix_def)
    9.16 +  apply (drule_tac j=j and i=i in spec2)
    9.17 +  apply (simp)
    9.18 +  apply (simp add: le_matrix_def)
    9.19 +  done
    9.20 +
    9.21 +lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)"
    9.22 +  apply (auto)
    9.23 +  apply (simp add: le_matrix_def)
    9.24 +  apply (drule_tac j=j and i=i in spec2)
    9.25 +  apply (simp)
    9.26 +  apply (simp add: le_matrix_def)
    9.27 +  done
    9.28 +
    9.29 +lemma move_matrix_le_zero[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))"
    9.30 +  apply (auto simp add: le_matrix_def neg_def)
    9.31 +  apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
    9.32 +  apply (auto)
    9.33 +  done
    9.34 +
    9.35 +lemma move_matrix_zero_le[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)"
    9.36 +  apply (auto simp add: le_matrix_def neg_def)
    9.37 +  apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
    9.38 +  apply (auto)
    9.39 +  done
    9.40 +
    9.41 +lemma move_matrix_le_move_matrix_iff[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= move_matrix B j i) = (A <= (B::('a::{order,zero}) matrix))"
    9.42 +  apply (auto simp add: le_matrix_def neg_def)
    9.43 +  apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
    9.44 +  apply (auto)
    9.45 +  done  
    9.46 +
    9.47  end
    10.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    10.2 +++ b/src/HOL/Matrix/MatrixLP.ML	Fri Sep 03 17:10:36 2004 +0200
    10.3 @@ -0,0 +1,83 @@
    10.4 +use "MatrixLP_gensimp.ML";
    10.5 +
    10.6 +signature MATRIX_LP = 
    10.7 +sig
    10.8 +  val lp_dual_estimate : string -> int -> thm
    10.9 +  val simplify : thm -> thm
   10.10 +end
   10.11 +
   10.12 +structure MatrixLP : MATRIX_LP =
   10.13 +struct
   10.14 +
   10.15 +val _ = SimprocsCodegen.use_simprocs_code sg geninput
   10.16 +
   10.17 +val simp_le_spmat = simp_SparseMatrix_le_spmat_'_mult__'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_bool'
   10.18 +val simp_add_spmat = simp_SparseMatrix_add_spmat_'_mult__'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''
   10.19 +val simp_diff_spmat = simp_SparseMatrix_diff_spmat_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''' 
   10.20 +val simp_abs_spmat = simp_SparseMatrix_abs_spmat_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''' 
   10.21 +val simp_mult_spmat = simp_SparseMatrix_mult_spmat_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''
   10.22 +val simp_sorted_sparse_matrix = simp_SparseMatrix_sorted_sparse_matrix
   10.23 +val simp_sorted_spvec = simp_SparseMatrix_sorted_spvec
   10.24 +	   
   10.25 +fun single_arg simp ct = snd (simp (snd (Thm.dest_comb ct)))
   10.26 +
   10.27 +fun double_arg simp ct =  
   10.28 +	let
   10.29 +	    val (a, y) = Thm.dest_comb ct
   10.30 +	    val (_, x) = Thm.dest_comb a
   10.31 +	in
   10.32 +	    snd (simp x y)
   10.33 +	end      
   10.34 +
   10.35 +fun spmc ct =
   10.36 +    (case term_of ct of	 
   10.37 +	 ((Const ("SparseMatrix.le_spmat", _)) $ _) =>
   10.38 +	 single_arg simp_le_spmat ct
   10.39 +       | ((Const ("SparseMatrix.add_spmat", _)) $ _) =>
   10.40 +	 single_arg simp_add_spmat ct
   10.41 +       | ((Const ("SparseMatrix.diff_spmat", _)) $ _ $ _) =>
   10.42 +	 double_arg simp_diff_spmat ct
   10.43 +       | ((Const ("SparseMatrix.abs_spmat", _)) $ _) =>
   10.44 +	 single_arg simp_abs_spmat ct
   10.45 +       | ((Const ("SparseMatrix.mult_spmat", _)) $ _ $ _) =>
   10.46 +	 double_arg simp_mult_spmat ct
   10.47 +       | ((Const ("SparseMatrix.sorted_sparse_matrix", _)) $ _) =>
   10.48 +	 single_arg simp_sorted_sparse_matrix ct
   10.49 +       | ((Const ("SparseMatrix.sorted_spvec", ty)) $ _) =>
   10.50 +	 single_arg simp_sorted_spvec ct 
   10.51 +      | _ => raise Fail_conv)
   10.52 +
   10.53 +fun lp_dual_estimate lptfile prec =
   10.54 +    let
   10.55 +	val th = inst_real (thm "SparseMatrix.spm_linprog_dual_estimate_1")
   10.56 +	val sg = sign_of (theory "MatrixLP")
   10.57 +	val (y, (A1, A2), (c1, c2), b, r) = 
   10.58 +	    let
   10.59 +		open fspmlp
   10.60 +		val l = load lptfile prec false
   10.61 +	    in
   10.62 +		(y l, A l, c l, b l, r l)
   10.63 +	    end		
   10.64 +	fun var s x = (cterm_of sg (Var ((s,0), FloatSparseMatrixBuilder.real_spmatT)), x)
   10.65 +	val th = Thm.instantiate ([], [var "A1" A1, var "A2" A2, var "y" y, var "c1" c1, var "c2" c2, var "r" r, var "b" b]) th
   10.66 +    in
   10.67 +	th
   10.68 +    end
   10.69 +
   10.70 +fun simplify th = 
   10.71 +    let
   10.72 +	val simp_th = botc spmc (cprop_of th)
   10.73 +	val th = strip_shyps (equal_elim simp_th th)
   10.74 +	fun removeTrue th = removeTrue (implies_elim th TrueI) handle _ => th
   10.75 +    in
   10.76 +	removeTrue th
   10.77 +    end
   10.78 +
   10.79 +fun float2real (x,y) = (Real.fromInt x) * (Math.pow (2.0, Real.fromInt y))
   10.80 +
   10.81 +end
   10.82 +
   10.83 +val result = ref ([]:thm list)
   10.84 +
   10.85 +fun run c = timeit (fn () => (result := [MatrixLP.simplify c]; ()))
   10.86 +
    11.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    11.2 +++ b/src/HOL/Matrix/MatrixLP.thy	Fri Sep 03 17:10:36 2004 +0200
    11.3 @@ -0,0 +1,4 @@
    11.4 +theory MatrixLP = Float + SparseMatrix:
    11.5 +
    11.6 +end
    11.7 +
    12.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    12.2 +++ b/src/HOL/Matrix/MatrixLP_gensimp.ML	Fri Sep 03 17:10:36 2004 +0200
    12.3 @@ -0,0 +1,35 @@
    12.4 +open BasisLibrary;
    12.5 +
    12.6 +use "Cplex.ML";
    12.7 +use "CplexMatrixConverter.ML";
    12.8 +use "ExactFloatingPoint.ML";
    12.9 +use "FloatArith.ML";
   12.10 +use "FloatSparseMatrixBuilder.ML";
   12.11 +use "fspmlp.ML";
   12.12 +use "codegen_prep.ML";
   12.13 +use "eq_codegen.ML";
   12.14 +use "conv.ML";
   12.15 +
   12.16 +
   12.17 +val th = theory "MatrixLP"
   12.18 +val sg = sign_of th
   12.19 +	 
   12.20 +fun prep x = standard (mk_meta_eq x)
   12.21 +	     
   12.22 +fun inst_real thm = standard (Thm.instantiate ([(fst (hd (term_tvars (prop_of thm))),
   12.23 +						 ctyp_of sg HOLogic.realT)], []) thm);
   12.24 +    
   12.25 +val spm_lp_dual_estimate_simp = 
   12.26 +    (thms "Float.arith_no_let") @ 
   12.27 +    (map inst_real (thms "SparseMatrix.sparse_row_matrix_arith_simps")) @ 
   12.28 +    (thms "SparseMatrix.sorted_sp_simps") @ 
   12.29 +    [thm "Product_Type.fst_conv", thm "Product_Type.snd_conv"] @ 
   12.30 +    (thms "SparseMatrix.boolarith") 
   12.31 +    
   12.32 +val simp_with = Simplifier.simplify (HOL_basic_ss addsimps [thm "SparseMatrix.if_case_eq", thm "Float.norm_0_1"])
   12.33 +		
   12.34 +val spm_lp_dual_estimate_simp = map simp_with spm_lp_dual_estimate_simp
   12.35 +				
   12.36 +val geninput = codegen_prep.prepare_thms spm_lp_dual_estimate_simp
   12.37 +
   12.38 +val _ = SimprocsCodegen.use_simprocs_code sg geninput
    13.1 --- a/src/HOL/Matrix/ROOT.ML	Fri Sep 03 10:27:05 2004 +0200
    13.2 +++ b/src/HOL/Matrix/ROOT.ML	Fri Sep 03 17:10:36 2004 +0200
    13.3 @@ -1,10 +1,1 @@
    13.4 -(*  Title:      HOL/Matrix/ROOT.ML
    13.5 -    ID:         $Id$
    13.6 -    Author:     Steven Obua
    13.7 -    License:    2004 Technische Universität München
    13.8 -
    13.9 -Theory of matrices with an application of matrix theory to linear
   13.10 -programming.
   13.11 -*)
   13.12 -
   13.13 -use_thy "LinProg";
   13.14 +use_thy "MatrixLP"
    14.1 --- a/src/HOL/Matrix/SparseMatrix.thy	Fri Sep 03 10:27:05 2004 +0200
    14.2 +++ b/src/HOL/Matrix/SparseMatrix.thy	Fri Sep 03 17:10:36 2004 +0200
    14.3 @@ -94,9 +94,62 @@
    14.4    done
    14.5  
    14.6  consts
    14.7 +  abs_spvec :: "('a::lordered_ring) spvec \<Rightarrow> 'a spvec"
    14.8 +  minus_spvec ::  "('a::lordered_ring) spvec \<Rightarrow> 'a spvec"
    14.9    smult_spvec :: "('a::lordered_ring) \<Rightarrow> 'a spvec \<Rightarrow> 'a spvec" 
   14.10    addmult_spvec :: "('a::lordered_ring) * 'a spvec * 'a spvec \<Rightarrow> 'a spvec"
   14.11  
   14.12 +primrec 
   14.13 +  "minus_spvec [] = []"
   14.14 +  "minus_spvec (a#as) = (fst a, -(snd a))#(minus_spvec as)"
   14.15 +
   14.16 +primrec 
   14.17 +  "abs_spvec [] = []"
   14.18 +  "abs_spvec (a#as) = (fst a, abs (snd a))#(abs_spvec as)"
   14.19 +
   14.20 +lemma sparse_row_vector_minus: 
   14.21 +  "sparse_row_vector (minus_spvec v) = - (sparse_row_vector v)"
   14.22 +  apply (induct v)
   14.23 +  apply (simp_all add: sparse_row_vector_cons)
   14.24 +  apply (simp add: Rep_matrix_inject[symmetric])
   14.25 +  apply (rule ext)+
   14.26 +  apply simp
   14.27 +  done
   14.28 +
   14.29 +lemma sparse_row_vector_abs:
   14.30 +  "sorted_spvec v \<Longrightarrow> sparse_row_vector (abs_spvec v) = abs (sparse_row_vector v)"
   14.31 +  apply (induct v)
   14.32 +  apply (simp_all add: sparse_row_vector_cons)
   14.33 +  apply (frule_tac sorted_spvec_cons1, simp)
   14.34 +  apply (simp only: Rep_matrix_inject[symmetric])
   14.35 +  apply (rule ext)+
   14.36 +  apply auto
   14.37 +  apply (subgoal_tac "Rep_matrix (sparse_row_vector list) 0 a = 0")
   14.38 +  apply (simp)
   14.39 +  apply (rule sorted_sparse_row_vector_zero)
   14.40 +  apply auto
   14.41 +  done
   14.42 +
   14.43 +lemma sorted_spvec_minus_spvec:
   14.44 +  "sorted_spvec v \<Longrightarrow> sorted_spvec (minus_spvec v)"
   14.45 +  apply (induct v)
   14.46 +  apply (simp)
   14.47 +  apply (frule sorted_spvec_cons1, simp)
   14.48 +  apply (simp add: sorted_spvec.simps)
   14.49 +  apply (case_tac list)
   14.50 +  apply (simp_all)
   14.51 +  done
   14.52 +
   14.53 +lemma sorted_spvec_abs_spvec:
   14.54 +  "sorted_spvec v \<Longrightarrow> sorted_spvec (abs_spvec v)"
   14.55 +  apply (induct v)
   14.56 +  apply (simp)
   14.57 +  apply (frule sorted_spvec_cons1, simp)
   14.58 +  apply (simp add: sorted_spvec.simps)
   14.59 +  apply (case_tac list)
   14.60 +  apply (simp_all)
   14.61 +  done
   14.62 +  
   14.63  defs
   14.64    smult_spvec_def: "smult_spvec y arr == map (% a. (fst a, y * snd a)) arr"  
   14.65  
   14.66 @@ -542,9 +595,6 @@
   14.67    else
   14.68      (le_spvec(snd a, snd b) & le_spmat (as, bs))))"
   14.69  
   14.70 -lemma spec2: "! j i. P j i \<Longrightarrow> P j i" by blast
   14.71 -lemma neg_imp: "(\<not> Q \<longrightarrow> \<not> P) \<Longrightarrow> P \<longrightarrow> Q" by blast
   14.72 -
   14.73  constdefs
   14.74    disj_matrices :: "('a::zero) matrix \<Rightarrow> 'a matrix \<Rightarrow> bool"
   14.75    "disj_matrices A B == (! j i. (Rep_matrix A j i \<noteq> 0) \<longrightarrow> (Rep_matrix B j i = 0)) & (! j i. (Rep_matrix B j i \<noteq> 0) \<longrightarrow> (Rep_matrix A j i = 0))"  
   14.76 @@ -597,22 +647,6 @@
   14.77    (B + A <= C) = (A <= C &  (B::('a::lordered_ab_group) matrix) <= 0)"
   14.78  by (auto simp add: disj_matrices_add[of B A 0 C,simplified] disj_matrices_commute)
   14.79  
   14.80 -lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))"
   14.81 -  apply (auto)
   14.82 -  apply (simp add: le_matrix_def)
   14.83 -  apply (drule_tac j=j and i=i in spec2)
   14.84 -  apply (simp)
   14.85 -  apply (simp add: le_matrix_def)
   14.86 -  done
   14.87 -
   14.88 -lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)"
   14.89 -  apply (auto)
   14.90 -  apply (simp add: le_matrix_def)
   14.91 -  apply (drule_tac j=j and i=i in spec2)
   14.92 -  apply (simp)
   14.93 -  apply (simp add: le_matrix_def)
   14.94 -  done
   14.95 -
   14.96  lemma disj_sparse_row_singleton: "i <= j \<Longrightarrow> sorted_spvec((j,y)#v) \<Longrightarrow> disj_matrices (sparse_row_vector v) (singleton_matrix 0 i x)"
   14.97    apply (simp add: disj_matrices_def)
   14.98    apply (rule conjI)
   14.99 @@ -641,27 +675,6 @@
  14.100  lemma disj_singleton_matrices[simp]: "disj_matrices (singleton_matrix j i x) (singleton_matrix u v y) = (j \<noteq> u | i \<noteq> v | x = 0 | y = 0)" 
  14.101    by (auto simp add: disj_matrices_def)
  14.102  
  14.103 -lemma le_spvec_iff_sparse_row_le[rule_format]: "(sorted_spvec a) \<longrightarrow> (sorted_spvec b) \<longrightarrow> (le_spvec (a,b)) = (sparse_row_vector a <= sparse_row_vector b)"
  14.104 -  apply (rule le_spvec.induct)
  14.105 -  apply (simp_all add: sorted_spvec_cons1 disj_matrices_add_le_zero disj_matrices_add_zero_le disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
  14.106 -  apply (rule conjI, intro strip)
  14.107 -  apply (simp add: sorted_spvec_cons1)
  14.108 -  apply (subst disj_matrices_add_x_le)
  14.109 -  apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute)
  14.110 -  apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
  14.111 -  apply (simp, blast)
  14.112 -  apply (intro strip, rule conjI, intro strip)
  14.113 -  apply (simp add: sorted_spvec_cons1)
  14.114 -  apply (subst disj_matrices_add_le_x)
  14.115 -  apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_sparse_row_singleton[OF less_imp_le] disj_matrices_commute disj_matrices_x_add)
  14.116 -  apply (blast)
  14.117 -  apply (intro strip)
  14.118 -  apply (simp add: sorted_spvec_cons1)
  14.119 -  apply (case_tac "a=aa", simp_all)
  14.120 -  apply (subst disj_matrices_add)
  14.121 -  apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
  14.122 -  done
  14.123 -
  14.124  lemma disj_move_sparse_vec_mat[simplified disj_matrices_commute]: 
  14.125    "j <= a \<Longrightarrow> sorted_spvec((a,c)#as) \<Longrightarrow> disj_matrices (move_matrix (sparse_row_vector b) (int j) i) (sparse_row_matrix as)"
  14.126    apply (auto simp add: neg_def disj_matrices_def)
  14.127 @@ -685,24 +698,28 @@
  14.128    apply (rule nrows, rule order_trans[of _ 1], simp, drule nrows_notzero, drule less_le_trans[OF _ nrows_spvec], arith)+
  14.129    done
  14.130  
  14.131 -lemma move_matrix_le_zero[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))"
  14.132 -  apply (auto simp add: le_matrix_def neg_def)
  14.133 -  apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
  14.134 -  apply (auto)
  14.135 +lemma le_spvec_iff_sparse_row_le[rule_format]: "(sorted_spvec a) \<longrightarrow> (sorted_spvec b) \<longrightarrow> (le_spvec (a,b)) = (sparse_row_vector a <= sparse_row_vector b)"
  14.136 +  apply (rule le_spvec.induct)
  14.137 +  apply (simp_all add: sorted_spvec_cons1 disj_matrices_add_le_zero disj_matrices_add_zero_le 
  14.138 +    disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
  14.139 +  apply (rule conjI, intro strip)
  14.140 +  apply (simp add: sorted_spvec_cons1)
  14.141 +  apply (subst disj_matrices_add_x_le)
  14.142 +  apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute)
  14.143 +  apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
  14.144 +  apply (simp, blast)
  14.145 +  apply (intro strip, rule conjI, intro strip)
  14.146 +  apply (simp add: sorted_spvec_cons1)
  14.147 +  apply (subst disj_matrices_add_le_x)
  14.148 +  apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_sparse_row_singleton[OF less_imp_le] disj_matrices_commute disj_matrices_x_add)
  14.149 +  apply (blast)
  14.150 +  apply (intro strip)
  14.151 +  apply (simp add: sorted_spvec_cons1)
  14.152 +  apply (case_tac "a=aa", simp_all)
  14.153 +  apply (subst disj_matrices_add)
  14.154 +  apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
  14.155    done
  14.156  
  14.157 -lemma move_matrix_zero_le[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)"
  14.158 -  apply (auto simp add: le_matrix_def neg_def)
  14.159 -  apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
  14.160 -  apply (auto)
  14.161 -  done
  14.162 -
  14.163 -lemma move_matrix_le_move_matrix_iff[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= move_matrix B j i) = (A <= (B::('a::{order,zero}) matrix))"
  14.164 -  apply (auto simp add: le_matrix_def neg_def)
  14.165 -  apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
  14.166 -  apply (auto)
  14.167 -  done  
  14.168 -
  14.169  lemma le_spvec_empty2_sparse_row[rule_format]: "(sorted_spvec b) \<longrightarrow> (le_spvec (b,[]) = (sparse_row_vector b <= 0))"
  14.170    apply (induct b)
  14.171    apply (simp_all add: sorted_spvec_cons1)
  14.172 @@ -752,39 +769,169 @@
  14.173    apply (simp add: sorted_spvec_cons1 le_spvec_iff_sparse_row_le)
  14.174    done
  14.175  
  14.176 -term smult_spvec
  14.177 -term addmult_spvec
  14.178 -term add_spvec
  14.179 -term mult_spvec_spmat
  14.180 -term mult_spmat
  14.181 -term add_spmat
  14.182 -term le_spvec
  14.183 -term le_spmat
  14.184 -term sorted_spvec
  14.185 -term sorted_spmat
  14.186 +ML {* simp_depth_limit := 999 *}
  14.187 +
  14.188 +consts 
  14.189 +   abs_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat"
  14.190 +   minus_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat"
  14.191 +
  14.192 +primrec
  14.193 +  "abs_spmat [] = []"
  14.194 +  "abs_spmat (a#as) = (fst a, abs_spvec (snd a))#(abs_spmat as)"
  14.195 +
  14.196 +primrec
  14.197 +  "minus_spmat [] = []"
  14.198 +  "minus_spmat (a#as) = (fst a, minus_spvec (snd a))#(minus_spmat as)"
  14.199 +
  14.200 +lemma sparse_row_matrix_minus:
  14.201 +  "sparse_row_matrix (minus_spmat A) = - (sparse_row_matrix A)"
  14.202 +  apply (induct A)
  14.203 +  apply (simp_all add: sparse_row_vector_minus sparse_row_matrix_cons)
  14.204 +  apply (subst Rep_matrix_inject[symmetric])
  14.205 +  apply (rule ext)+
  14.206 +  apply simp
  14.207 +  done
  14.208  
  14.209 -thm sparse_row_mult_spmat
  14.210 -thm sparse_row_add_spmat
  14.211 -thm le_spmat_iff_sparse_row_le
  14.212 +lemma Rep_sparse_row_vector_zero: "x \<noteq> 0 \<Longrightarrow> Rep_matrix (sparse_row_vector v) x y = 0"
  14.213 +proof -
  14.214 +  assume x:"x \<noteq> 0"
  14.215 +  have r:"nrows (sparse_row_vector v) <= Suc 0" by (rule nrows_spvec)
  14.216 +  show ?thesis
  14.217 +    apply (rule nrows)
  14.218 +    apply (subgoal_tac "Suc 0 <= x")
  14.219 +    apply (insert r)
  14.220 +    apply (simp only:)
  14.221 +    apply (insert x)
  14.222 +    apply arith
  14.223 +    done
  14.224 +qed
  14.225 +    
  14.226 +lemma sparse_row_matrix_abs:
  14.227 +  "sorted_spvec A \<Longrightarrow> sorted_spmat A \<Longrightarrow> sparse_row_matrix (abs_spmat A) = abs (sparse_row_matrix A)"
  14.228 +  apply (induct A)
  14.229 +  apply (simp_all add: sparse_row_vector_abs sparse_row_matrix_cons)
  14.230 +  apply (frule_tac sorted_spvec_cons1, simp)
  14.231 +  apply (subst Rep_matrix_inject[symmetric])
  14.232 +  apply (rule ext)+
  14.233 +  apply auto
  14.234 +  apply (case_tac "x=a")
  14.235 +  apply (simp)
  14.236 +  apply (subst sorted_sparse_row_matrix_zero)
  14.237 +  apply auto
  14.238 +  apply (subst Rep_sparse_row_vector_zero)
  14.239 +  apply (simp_all add: neg_def)
  14.240 +  done
  14.241 +
  14.242 +lemma sorted_spvec_minus_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec (minus_spmat A)"
  14.243 +  apply (induct A)
  14.244 +  apply (simp)
  14.245 +  apply (frule sorted_spvec_cons1, simp)
  14.246 +  apply (simp add: sorted_spvec.simps)
  14.247 +  apply (case_tac list)
  14.248 +  apply (simp_all)
  14.249 +  done 
  14.250 +
  14.251 +lemma sorted_spvec_abs_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec (abs_spmat A)" 
  14.252 +  apply (induct A)
  14.253 +  apply (simp)
  14.254 +  apply (frule sorted_spvec_cons1, simp)
  14.255 +  apply (simp add: sorted_spvec.simps)
  14.256 +  apply (case_tac list)
  14.257 +  apply (simp_all)
  14.258 +  done
  14.259 +
  14.260 +lemma sorted_spmat_minus_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat (minus_spmat A)"
  14.261 +  apply (induct A)
  14.262 +  apply (simp_all add: sorted_spvec_minus_spvec)
  14.263 +  done
  14.264 +
  14.265 +lemma sorted_spmat_abs_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat (abs_spmat A)"
  14.266 +  apply (induct A)
  14.267 +  apply (simp_all add: sorted_spvec_abs_spvec)
  14.268 +  done
  14.269  
  14.270 -thm sorted_spvec_mult_spmat
  14.271 -thm sorted_spmat_mult_spmat
  14.272 -thm sorted_spvec_add_spmat
  14.273 -thm sorted_spmat_add_spmat
  14.274 +constdefs
  14.275 +  diff_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat \<Rightarrow> 'a spmat"
  14.276 +  "diff_spmat A B == add_spmat (A, minus_spmat B)"
  14.277 +
  14.278 +lemma sorted_spmat_diff_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat B \<Longrightarrow> sorted_spmat (diff_spmat A B)"
  14.279 +  by (simp add: diff_spmat_def sorted_spmat_minus_spmat sorted_spmat_add_spmat)
  14.280 +
  14.281 +lemma sorted_spvec_diff_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec B \<Longrightarrow> sorted_spvec (diff_spmat A B)"
  14.282 +  by (simp add: diff_spmat_def sorted_spvec_minus_spmat sorted_spvec_add_spmat)
  14.283 +
  14.284 +lemma sparse_row_diff_spmat: "sparse_row_matrix (diff_spmat A B ) = (sparse_row_matrix A) - (sparse_row_matrix B)"
  14.285 +  by (simp add: diff_spmat_def sparse_row_add_spmat sparse_row_matrix_minus)
  14.286 +
  14.287 +constdefs
  14.288 +  sorted_sparse_matrix :: "'a spmat \<Rightarrow> bool"
  14.289 +  "sorted_sparse_matrix A == (sorted_spvec A) & (sorted_spmat A)"
  14.290 +
  14.291 +lemma sorted_sparse_matrix_imp_spvec: "sorted_sparse_matrix A \<Longrightarrow> sorted_spvec A"
  14.292 +  by (simp add: sorted_sparse_matrix_def)
  14.293 +
  14.294 +lemma sorted_sparse_matrix_imp_spmat: "sorted_sparse_matrix A \<Longrightarrow> sorted_spmat A"
  14.295 +  by (simp add: sorted_sparse_matrix_def)
  14.296 +
  14.297 +lemmas sparse_row_matrix_op_simps =
  14.298 +  sorted_sparse_matrix_imp_spmat sorted_sparse_matrix_imp_spvec
  14.299 +  sparse_row_add_spmat sorted_spvec_add_spmat sorted_spmat_add_spmat
  14.300 +  sparse_row_diff_spmat sorted_spvec_diff_spmat sorted_spmat_diff_spmat
  14.301 +  sparse_row_matrix_minus sorted_spvec_minus_spmat sorted_spmat_minus_spmat
  14.302 +  sparse_row_mult_spmat sorted_spvec_mult_spmat sorted_spmat_mult_spmat
  14.303 +  sparse_row_matrix_abs sorted_spvec_abs_spmat sorted_spmat_abs_spmat
  14.304 +  le_spmat_iff_sparse_row_le
  14.305 +
  14.306 +lemma zero_eq_Numeral0: "(0::_::number_ring) = Numeral0" by simp
  14.307  
  14.308 -thm smult_spvec_empty
  14.309 -thm smult_spvec_cons
  14.310 -thm addmult_spvec.simps
  14.311 -thm add_spvec.simps
  14.312 -thm add_spmat.simps
  14.313 -thm mult_spvec_spmat.simps
  14.314 -thm mult_spmat.simps
  14.315 -thm le_spvec.simps
  14.316 -thm le_spmat.simps
  14.317 -thm sorted_spvec.simps
  14.318 -thm sorted_spmat.simps
  14.319 +lemmas sparse_row_matrix_arith_simps[simplified zero_eq_Numeral0] = 
  14.320 +  mult_spmat.simps mult_spvec_spmat.simps 
  14.321 +  addmult_spvec.simps 
  14.322 +  smult_spvec_empty smult_spvec_cons
  14.323 +  add_spmat.simps add_spvec.simps
  14.324 +  minus_spmat.simps minus_spvec.simps
  14.325 +  abs_spmat.simps abs_spvec.simps
  14.326 +  diff_spmat_def
  14.327 +  le_spmat.simps le_spvec.simps
  14.328 +
  14.329 +lemmas sorted_sp_simps = 
  14.330 +  sorted_spvec.simps
  14.331 +  sorted_spmat.simps
  14.332 +  sorted_sparse_matrix_def
  14.333 +
  14.334 +lemma bool1: "(\<not> True) = False"  by blast
  14.335 +lemma bool2: "(\<not> False) = True"  by blast
  14.336 +lemma bool3: "((P\<Colon>bool) \<and> True) = P" by blast
  14.337 +lemma bool4: "(True \<and> (P\<Colon>bool)) = P" by blast
  14.338 +lemma bool5: "((P\<Colon>bool) \<and> False) = False" by blast
  14.339 +lemma bool6: "(False \<and> (P\<Colon>bool)) = False" by blast
  14.340 +lemma bool7: "((P\<Colon>bool) \<or> True) = True" by blast
  14.341 +lemma bool8: "(True \<or> (P\<Colon>bool)) = True" by blast
  14.342 +lemma bool9: "((P\<Colon>bool) \<or> False) = P" by blast
  14.343 +lemma bool10: "(False \<or> (P\<Colon>bool)) = P" by blast
  14.344 +lemmas boolarith = bool1 bool2 bool3 bool4 bool5 bool6 bool7 bool8 bool9 bool10
  14.345 +
  14.346 +lemma if_case_eq: "(if b then x else y) = (case b of True => x | False => y)" by simp
  14.347 +
  14.348 +lemma spm_linprog_dual_estimate_1:
  14.349 +  assumes  
  14.350 +  "sorted_sparse_matrix A1"
  14.351 +  "sorted_sparse_matrix A2"
  14.352 +  "sorted_sparse_matrix c1"
  14.353 +  "sorted_sparse_matrix c2"
  14.354 +  "sorted_sparse_matrix y"
  14.355 +  "sorted_spvec b"
  14.356 +  "sorted_spvec r"
  14.357 +  "le_spmat ([], y)"
  14.358 +  "A * x \<le> sparse_row_matrix (b::('a::lordered_ring) spmat)"
  14.359 +  "sparse_row_matrix A1 <= A"
  14.360 +  "A <= sparse_row_matrix A2"
  14.361 +  "sparse_row_matrix c1 <= c"
  14.362 +  "c <= sparse_row_matrix c2"
  14.363 +  "abs x \<le> sparse_row_matrix r"
  14.364 +  shows
  14.365 +  "c * x \<le> sparse_row_matrix (add_spmat (mult_spmat y b, mult_spmat (add_spmat (add_spmat (mult_spmat y (diff_spmat A2 A1), 
  14.366 +  abs_spmat (diff_spmat (mult_spmat y A1) c1)), diff_spmat c2 c1)) r))"
  14.367 +  by (insert prems, simp add: sparse_row_matrix_op_simps linprog_dual_estimate_1[where A=A])
  14.368  
  14.369  end
  14.370 -
  14.371 -
  14.372 -
    15.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    15.2 +++ b/src/HOL/Matrix/codegen_prep.ML	Fri Sep 03 17:10:36 2004 +0200
    15.3 @@ -0,0 +1,110 @@
    15.4 +structure codegen_prep = 
    15.5 +struct
    15.6 +
    15.7 +exception Prepare_thms of string;
    15.8 +
    15.9 +fun is_meta_eq th =
   15.10 +    let
   15.11 +	fun check ((Const ("==", _)) $ _ $ _) = true
   15.12 +	  | check _ = false
   15.13 +    in
   15.14 +	check (concl_of th)
   15.15 +    end
   15.16 +
   15.17 +fun printty ty = 
   15.18 +    let
   15.19 +	fun immerse s [] = []
   15.20 +	  | immerse s [x] = [x]
   15.21 +	  | immerse s (x::xs) = x::s::(immerse s xs)
   15.22 +
   15.23 +	fun py t = 
   15.24 +	    let
   15.25 +		val (name, args) =  dest_Type t
   15.26 +		val args' = map printty args
   15.27 +	    in
   15.28 +		concat (immerse "_" (name::args'))
   15.29 +	    end
   15.30 +
   15.31 +	val (args, res) = strip_type ty
   15.32 +	val tystrs = map py (args@[res])
   15.33 +			  
   15.34 +	val s = "'"^(concat (immerse "_" tystrs))^"'"
   15.35 +    in
   15.36 +	s
   15.37 +    end
   15.38 +
   15.39 +fun head_name th = 
   15.40 +    let val s = 
   15.41 +	    let
   15.42 +		val h = fst (strip_comb (hd (snd (strip_comb (concl_of th)))))
   15.43 +		val n = fst (dest_Const h)
   15.44 +		val ty = snd (dest_Const h)
   15.45 +		val hn = fst (dest_Const h)
   15.46 +	    in
   15.47 +		hn^"_"^(printty ty) handle _ => (writeln ("warning: polymorphic "^hn); hn)	
   15.48 +	    end
   15.49 +    in
   15.50 +	s
   15.51 +    end
   15.52 +	    
   15.53 +val mangle_id  = 
   15.54 +    let
   15.55 +	fun tr #"=" = "_eq_"
   15.56 +	  | tr #"." = "_"
   15.57 +	  | tr #" " = "_"
   15.58 +	  | tr #"<" = "_le_"
   15.59 +	  | tr #">" = "_ge_"
   15.60 +	  | tr #"(" = "_bro_"
   15.61 +	  | tr #")" = "_brc_"
   15.62 +	  | tr #"+" = "_plus_"
   15.63 +	  | tr #"*" = "_mult_"
   15.64 +	  | tr #"-" = "_minus_"
   15.65 +	  | tr #"|" = "_or_"
   15.66 +	  | tr #"&" = "_and_"
   15.67 +	  | tr x = Char.toString x
   15.68 +	fun m s = "simp_"^(String.translate tr s)
   15.69 +    in
   15.70 +	m
   15.71 +    end
   15.72 +
   15.73 +fun group f [] = []
   15.74 +  | group f (x::xs) = 
   15.75 +    let
   15.76 +	val g = group f xs
   15.77 +	val key = f x
   15.78 +    in
   15.79 +	case assoc (g, key) of
   15.80 +	    None => (key, [x])::g
   15.81 +	  | Some l => overwrite (g, (key, x::l))
   15.82 +    end
   15.83 +    
   15.84 +fun prepare_thms ths = 
   15.85 +    let
   15.86 +	val ths = (filter is_meta_eq ths) @ (map (standard o mk_meta_eq) (filter (not o is_meta_eq) ths))		      
   15.87 +	val _ = if forall Thm.no_prems ths then () else raise (Prepare_thms "premisse found")
   15.88 +	val thmgroups = group head_name ths
   15.89 +	val test_clashgroups = group (fn (a,b) => mangle_id a) thmgroups
   15.90 +	val _ = if (length thmgroups <> length test_clashgroups) then raise (Prepare_thms "clash after name mangling") else ()	
   15.91 +		
   15.92 +	fun prep (name, ths) = 
   15.93 +	    let
   15.94 +		val m = mangle_id name
   15.95 +			
   15.96 +	    in
   15.97 +		(m, true, ths)
   15.98 +	    end
   15.99 +	    
  15.100 +	val thmgroups = map prep thmgroups
  15.101 +    in
  15.102 +	(thmgroups)
  15.103 +    end
  15.104 +    
  15.105 +fun writestr filename s = 
  15.106 +    let 
  15.107 +	val f = TextIO.openOut filename
  15.108 +	val _ = TextIO.output(f, s)
  15.109 +	val _ = TextIO.closeOut f
  15.110 +    in 
  15.111 +	() 
  15.112 +    end
  15.113 +end
    16.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    16.2 +++ b/src/HOL/Matrix/conv.ML	Fri Sep 03 17:10:36 2004 +0200
    16.3 @@ -0,0 +1,26 @@
    16.4 +exception Fail_conv;
    16.5 +
    16.6 +fun orelsec conv1 conv2 ct = conv1 ct handle Fail_conv => conv2 ct
    16.7 +
    16.8 +val allc = Thm.reflexive
    16.9 +
   16.10 +fun thenc conv1 conv2 ct = 
   16.11 +    let 
   16.12 +	fun rhs_of t = snd (Thm.dest_comb (strip_imp_concl (cprop_of t)))
   16.13 +	    
   16.14 +	val eq = conv1 ct
   16.15 +    in
   16.16 +	Thm.transitive eq (conv2 (rhs_of eq))
   16.17 +    end
   16.18 +
   16.19 +fun subc conv ct = 
   16.20 +    case term_of ct of
   16.21 +	_ $ _ => 
   16.22 +	let 
   16.23 +	    val (ct1, ct2) = Thm.dest_comb ct
   16.24 +	in
   16.25 +	    Thm.combination (conv ct1) (conv ct2)
   16.26 +	end
   16.27 +      | _ => allc ct
   16.28 +
   16.29 +fun botc conv ct = thenc (subc (botc conv)) (orelsec conv allc) ct
   16.30 \ No newline at end of file
    17.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    17.2 +++ b/src/HOL/Matrix/eq_codegen.ML	Fri Sep 03 17:10:36 2004 +0200
    17.3 @@ -0,0 +1,493 @@
    17.4 +fun inst_cterm inst ct = fst (Drule.dest_equals
    17.5 +  (Thm.cprop_of (Thm.instantiate inst (reflexive ct))));
    17.6 +fun tyinst_cterm tyinst = inst_cterm (tyinst, []);
    17.7 +
    17.8 +val bla = ref ([] : term list);
    17.9 +
   17.10 +(******************************************************)
   17.11 +(*        Code generator for equational proofs        *)
   17.12 +(******************************************************)
   17.13 +fun my_mk_meta_eq thm =
   17.14 +  let
   17.15 +    val (_, eq) = Thm.dest_comb (cprop_of thm);
   17.16 +    val (ct, rhs) = Thm.dest_comb eq;
   17.17 +    val (_, lhs) = Thm.dest_comb ct
   17.18 +  in Thm.implies_elim (Drule.instantiate' [Some (ctyp_of_term lhs)]
   17.19 +    [Some lhs, Some rhs] eq_reflection) thm
   17.20 +  end; 
   17.21 +
   17.22 +structure SimprocsCodegen =
   17.23 +struct
   17.24 +
   17.25 +val simp_thms = ref ([] : thm list);
   17.26 +
   17.27 +fun parens b = if b then Pretty.enclose "(" ")" else Pretty.block;
   17.28 +
   17.29 +fun gen_mk_val f xs ps = Pretty.block ([Pretty.str "val ",
   17.30 +  f (length xs > 1) (flat
   17.31 +    (separate [Pretty.str ",", Pretty.brk 1] (map (single o Pretty.str) xs))),
   17.32 +  Pretty.str " =", Pretty.brk 1] @ ps @ [Pretty.str ";"]);
   17.33 +
   17.34 +val mk_val = gen_mk_val parens;
   17.35 +val mk_vall = gen_mk_val (K (Pretty.enclose "[" "]"));
   17.36 +
   17.37 +fun rename s = if s mem ThmDatabase.ml_reserved then s ^ "'" else s;
   17.38 +
   17.39 +fun mk_decomp_name (Var ((s, i), _)) = rename (if i=0 then s else s ^ string_of_int i)
   17.40 +  | mk_decomp_name (Const (s, _)) = rename (Codegen.mk_id (Sign.base_name s))
   17.41 +  | mk_decomp_name _ = "ct";
   17.42 +
   17.43 +fun decomp_term_code cn ((vs, bs, ps), (v, t)) =
   17.44 +  if exists (equal t o fst) bs then (vs, bs, ps)
   17.45 +  else (case t of
   17.46 +      Var _ => (vs, bs @ [(t, v)], ps)
   17.47 +    | Const _ => (vs, if cn then bs @ [(t, v)] else bs, ps)
   17.48 +    | Bound _ => (vs, bs, ps)
   17.49 +    | Abs (s, T, t) =>
   17.50 +      let
   17.51 +        val v1 = variant vs s;
   17.52 +        val v2 = variant (v1 :: vs) (mk_decomp_name t)
   17.53 +      in
   17.54 +        decomp_term_code cn ((v1 :: v2 :: vs,
   17.55 +          bs @ [(Free (s, T), v1)],
   17.56 +          ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_abs", Pretty.brk 1,
   17.57 +            Pretty.str "None", Pretty.brk 1, Pretty.str v]]), (v2, t))
   17.58 +      end
   17.59 +    | t $ u =>
   17.60 +      let
   17.61 +        val v1 = variant vs (mk_decomp_name t);
   17.62 +        val v2 = variant (v1 :: vs) (mk_decomp_name u);
   17.63 +        val (vs', bs', ps') = decomp_term_code cn ((v1 :: v2 :: vs, bs,
   17.64 +          ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_comb", Pretty.brk 1,
   17.65 +            Pretty.str v]]), (v1, t));
   17.66 +        val (vs'', bs'', ps'') = decomp_term_code cn ((vs', bs', ps'), (v2, u))
   17.67 +      in
   17.68 +        if bs'' = bs then (vs, bs, ps) else (vs'', bs'', ps'')
   17.69 +      end);
   17.70 +
   17.71 +val strip_tv = implode o tl o explode;
   17.72 +
   17.73 +fun mk_decomp_tname (TVar ((s, i), _)) =
   17.74 +      strip_tv ((if i=0 then s else s ^ string_of_int i) ^ "T")
   17.75 +  | mk_decomp_tname (Type (s, _)) = Codegen.mk_id (Sign.base_name s) ^ "T"
   17.76 +  | mk_decomp_tname _ = "cT";
   17.77 +
   17.78 +fun decomp_type_code ((vs, bs, ps), (v, TVar (ixn, _))) =
   17.79 +      if exists (equal ixn o fst) bs then (vs, bs, ps)
   17.80 +      else (vs, bs @ [(ixn, v)], ps)
   17.81 +  | decomp_type_code ((vs, bs, ps), (v, Type (_, Ts))) =
   17.82 +      let
   17.83 +        val vs' = variantlist (map mk_decomp_tname Ts, vs);
   17.84 +        val (vs'', bs', ps') =
   17.85 +          foldl decomp_type_code ((vs @ vs', bs, ps @
   17.86 +            [mk_vall vs' [Pretty.str "Thm.dest_ctyp", Pretty.brk 1,
   17.87 +              Pretty.str v]]), vs' ~~ Ts)
   17.88 +      in
   17.89 +        if bs' = bs then (vs, bs, ps) else (vs'', bs', ps')
   17.90 +      end;
   17.91 +
   17.92 +fun gen_mk_bindings s dest decomp ((vs, bs, ps), (v, x)) =
   17.93 +  let
   17.94 +    val s' = variant vs s;
   17.95 +    val (vs', bs', ps') = decomp ((s' :: vs, bs, ps @
   17.96 +      [mk_val [s'] (dest v)]), (s', x))
   17.97 +  in
   17.98 +    if bs' = bs then (vs, bs, ps) else (vs', bs', ps')
   17.99 +  end;
  17.100 +
  17.101 +val mk_term_bindings = gen_mk_bindings "ct"
  17.102 +  (fn s => [Pretty.str "cprop_of", Pretty.brk 1, Pretty.str s])
  17.103 +  (decomp_term_code true);
  17.104 +
  17.105 +val mk_type_bindings = gen_mk_bindings "cT"
  17.106 +  (fn s => [Pretty.str "Thm.ctyp_of_term", Pretty.brk 1, Pretty.str s])
  17.107 +  decomp_type_code;
  17.108 +
  17.109 +fun pretty_pattern b (Const (s, _)) = Pretty.block [Pretty.str "Const",
  17.110 +      Pretty.brk 1, Pretty.str ("(\"" ^ s ^ "\", _)")]
  17.111 +  | pretty_pattern b (t as _ $ _) = parens b
  17.112 +      (flat (separate [Pretty.str " $", Pretty.brk 1]
  17.113 +        (map (single o pretty_pattern true) (op :: (strip_comb t)))))
  17.114 +  | pretty_pattern b _ = Pretty.str "_";
  17.115 +
  17.116 +fun term_consts' t = foldl_aterms
  17.117 +  (fn (cs, c as Const _) => c ins cs | (cs, _) => cs) ([], t);
  17.118 +
  17.119 +fun mk_apps s b p [] = p
  17.120 +  | mk_apps s b p (q :: qs) = 
  17.121 +      mk_apps s b (parens (b orelse not (null qs))
  17.122 +        [Pretty.str s, Pretty.brk 1, p, Pretty.brk 1, q]) qs;
  17.123 +
  17.124 +fun mk_refleq eq ct = mk_val [eq] [Pretty.str ("Thm.reflexive " ^ ct)];
  17.125 +
  17.126 +fun mk_tyinst ((s, i), s') =
  17.127 +  Pretty.block [Pretty.str ("((" ^ quote s ^ ","), Pretty.brk 1,
  17.128 +    Pretty.str (string_of_int i ^ "),"), Pretty.brk 1,
  17.129 +    Pretty.str (s' ^ ")")];
  17.130 +
  17.131 +fun inst_ty b ty_bs t s = (case term_tvars t of
  17.132 +    [] => Pretty.str s
  17.133 +  | Ts => parens b [Pretty.str "tyinst_cterm", Pretty.brk 1,
  17.134 +      Pretty.list "[" "]" (map (fn (ixn, _) => mk_tyinst
  17.135 +        (ixn, the (assoc (ty_bs, ixn)))) Ts),
  17.136 +      Pretty.brk 1, Pretty.str s]);
  17.137 +
  17.138 +fun mk_cterm_code b ty_bs ts xs (vals, t $ u) =
  17.139 +      let
  17.140 +        val (vals', p1) = mk_cterm_code true ty_bs ts xs (vals, t);
  17.141 +        val (vals'', p2) = mk_cterm_code true ty_bs ts xs (vals', u)
  17.142 +      in
  17.143 +        (vals'', parens b [Pretty.str "Thm.capply", Pretty.brk 1,
  17.144 +          p1, Pretty.brk 1, p2])
  17.145 +      end
  17.146 +  | mk_cterm_code b ty_bs ts xs (vals, Abs (s, T, t)) =
  17.147 +      let
  17.148 +        val u = Free (s, T);
  17.149 +        val Some s' = assoc (ts, u);
  17.150 +        val p = Pretty.str s';
  17.151 +        val (vals', p') = mk_cterm_code true ty_bs ts (p :: xs)
  17.152 +          (if null (typ_tvars T) then vals
  17.153 +           else vals @ [(u, (("", s'), [mk_val [s'] [inst_ty true ty_bs u s']]))], t)
  17.154 +      in (vals',
  17.155 +        parens b [Pretty.str "Thm.cabs", Pretty.brk 1, p, Pretty.brk 1, p'])
  17.156 +      end
  17.157 +  | mk_cterm_code b ty_bs ts xs (vals, Bound i) = (vals, nth_elem (i, xs))
  17.158 +  | mk_cterm_code b ty_bs ts xs (vals, t) = (case assoc (vals, t) of
  17.159 +        None =>
  17.160 +          let val Some s = assoc (ts, t)
  17.161 +          in (if is_Const t andalso not (null (term_tvars t)) then
  17.162 +              vals @ [(t, (("", s), [mk_val [s] [inst_ty true ty_bs t s]]))]
  17.163 +            else vals, Pretty.str s)
  17.164 +          end
  17.165 +      | Some ((_, s), _) => (vals, Pretty.str s));
  17.166 +
  17.167 +fun get_cases sg =
  17.168 +  Symtab.foldl (fn (tab, (k, {case_rewrites, ...})) => Symtab.update_new
  17.169 +    ((fst (dest_Const (head_of (fst (HOLogic.dest_eq (HOLogic.dest_Trueprop
  17.170 +      (prop_of (hd case_rewrites))))))), map my_mk_meta_eq case_rewrites), tab))
  17.171 +        (Symtab.empty, DatatypePackage.get_datatypes_sg sg);
  17.172 +
  17.173 +fun decomp_case th =
  17.174 +  let
  17.175 +    val (lhs, _) = Logic.dest_equals (prop_of th);
  17.176 +    val (f, ts) = strip_comb lhs;
  17.177 +    val (us, u) = split_last ts;
  17.178 +    val (Const (s, _), vs) = strip_comb u
  17.179 +  in (us, s, vs, u) end;
  17.180 +
  17.181 +fun rename vs t =
  17.182 +  let
  17.183 +    fun mk_subst ((vs, subs), Var ((s, i), T)) =
  17.184 +      let val s' = variant vs s
  17.185 +      in if s = s' then (vs, subs)
  17.186 +        else (s' :: vs, ((s, i), Var ((s', i), T)) :: subs)
  17.187 +      end;
  17.188 +    val (vs', subs) = foldl mk_subst ((vs, []), term_vars t)
  17.189 +  in (vs', subst_Vars subs t) end;
  17.190 +
  17.191 +fun is_instance sg t u = t = subst_TVars_Vartab
  17.192 +  (Type.typ_match (Sign.tsig_of sg) (Vartab.empty,
  17.193 +    (fastype_of u, fastype_of t))) u handle Type.TYPE_MATCH => false;
  17.194 +
  17.195 +(*
  17.196 +fun lookup sg fs t = apsome snd (Library.find_first
  17.197 +  (is_instance sg t o fst) fs);
  17.198 +*)
  17.199 +
  17.200 +fun lookup sg fs t = (case Library.find_first (is_instance sg t o fst) fs of
  17.201 +    None => (bla := (t ins !bla); None)
  17.202 +  | Some (_, x) => Some x);
  17.203 +
  17.204 +fun unint sg fs t = forall (is_none o lookup sg fs) (term_consts' t);
  17.205 +
  17.206 +fun mk_let s i xs ys =
  17.207 +  Pretty.blk (0, [Pretty.blk (i, separate Pretty.fbrk (Pretty.str s :: xs)),
  17.208 +    Pretty.fbrk,
  17.209 +    Pretty.blk (i, ([Pretty.str "in", Pretty.fbrk] @ ys)),
  17.210 +    Pretty.fbrk, Pretty.str "end"]);
  17.211 +
  17.212 +(*****************************************************************************)
  17.213 +(* Generate bindings for simplifying term t                                  *)
  17.214 +(* mkeq: whether to generate reflexivity theorem for uninterpreted terms     *)
  17.215 +(* fs:   interpreted functions                                               *)
  17.216 +(* ts:   atomic terms                                                        *)
  17.217 +(* vs:   used identifiers                                                    *)
  17.218 +(* vals: list of bindings of the form ((eq, ct), ps) where                   *)
  17.219 +(*       eq: name of equational theorem                                      *)
  17.220 +(*       ct: name of simplified cterm                                        *)
  17.221 +(*       ps: ML code for creating the above two items                        *)
  17.222 +(*****************************************************************************)
  17.223 +
  17.224 +fun mk_simpl_code sg case_tab mkeq fs ts ty_bs thm_bs ((vs, vals), t) =
  17.225 +  (case assoc (vals, t) of
  17.226 +    Some ((eq, ct), ps) =>  (* binding already generated *) 
  17.227 +      if mkeq andalso eq="" then
  17.228 +        let val eq' = variant vs "eq"
  17.229 +        in ((eq' :: vs, overwrite (vals,
  17.230 +          (t, ((eq', ct), ps @ [mk_refleq eq' ct])))), (eq', ct))
  17.231 +        end
  17.232 +      else ((vs, vals), (eq, ct))
  17.233 +  | None => (case assoc (ts, t) of
  17.234 +      Some v =>  (* atomic term *)
  17.235 +        let val xs = if not (null (term_tvars t)) andalso is_Const t then
  17.236 +          [mk_val [v] [inst_ty false ty_bs t v]] else []
  17.237 +        in
  17.238 +          if mkeq then
  17.239 +            let val eq = variant vs "eq"
  17.240 +            in ((eq :: vs, vals @
  17.241 +              [(t, ((eq, v), xs @ [mk_refleq eq v]))]), (eq, v))
  17.242 +            end
  17.243 +          else ((vs, if null xs then vals else vals @
  17.244 +            [(t, (("", v), xs))]), ("", v))
  17.245 +        end
  17.246 +    | None =>  (* complex term *)
  17.247 +        let val (f as Const (cname, _), us) = strip_comb t
  17.248 +        in case Symtab.lookup (case_tab, cname) of
  17.249 +            Some cases =>  (* case expression *)
  17.250 +              let
  17.251 +                val (us', u) = split_last us;
  17.252 +                val b = unint sg fs u;
  17.253 +                val ((vs1, vals1), (eq, ct)) =
  17.254 +                  mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs ((vs, vals), u);
  17.255 +                val xs = variantlist (replicate (length us') "f", vs1);
  17.256 +                val (vals2, ps) = foldl_map
  17.257 +                  (mk_cterm_code false ty_bs ts []) (vals1, us');
  17.258 +                val fvals = map (fn (x, p) => mk_val [x] [p]) (xs ~~ ps);
  17.259 +                val uT = fastype_of u;
  17.260 +                val (us'', _, _, u') = decomp_case (hd cases);
  17.261 +                val (vs2, ty_bs', ty_vals) = mk_type_bindings
  17.262 +                  (mk_type_bindings ((vs1 @ xs, [], []),
  17.263 +                    (hd xs, fastype_of (hd us''))), (ct, fastype_of u'));
  17.264 +                val insts1 = map mk_tyinst ty_bs';
  17.265 +                val i = length vals2;
  17.266 +   
  17.267 +                fun mk_case_code ((vs, vals), (f, (name, eqn))) =
  17.268 +                  let
  17.269 +                    val (fvs, cname, cvs, _) = decomp_case eqn;
  17.270 +                    val Ts = binder_types (fastype_of f);
  17.271 +                    val ys = variantlist (map (fst o fst o dest_Var) cvs, vs);
  17.272 +                    val cvs' = map Var (map (rpair 0) ys ~~ Ts);
  17.273 +                    val rs = cvs' ~~ cvs;
  17.274 +                    val lhs = list_comb (Const (cname, Ts ---> uT), cvs');
  17.275 +                    val rhs = foldl betapply (f, cvs');
  17.276 +                    val (vs', tm_bs, tm_vals) = decomp_term_code false
  17.277 +                      ((vs @ ys, [], []), (ct, lhs));
  17.278 +                    val ((vs'', all_vals), (eq', ct')) = mk_simpl_code sg case_tab
  17.279 +                      false fs (tm_bs @ ts) ty_bs thm_bs ((vs', vals), rhs);
  17.280 +                    val (old_vals, eq_vals) = splitAt (i, all_vals);
  17.281 +                    val vs''' = vs @ filter (fn v => exists
  17.282 +                      (fn (_, ((v', _), _)) => v = v') old_vals) (vs'' \\ vs');
  17.283 +                    val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(",
  17.284 +                      inst_ty false ty_bs' t (the (assoc (thm_bs, t))), Pretty.str ",",
  17.285 +                      Pretty.brk 1, Pretty.str (s ^ ")")]) ((fvs ~~ xs) @
  17.286 +                        (map (fn (v, s) => (the (assoc (rs, v)), s)) tm_bs));
  17.287 +                    val eq'' = if null insts1 andalso null insts2 then Pretty.str name
  17.288 +                      else parens (eq' <> "") [Pretty.str
  17.289 +                          (if null cvs then "Thm.instantiate" else "Drule.instantiate"),
  17.290 +                        Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1,
  17.291 +                        Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2,
  17.292 +                        Pretty.str ")", Pretty.brk 1, Pretty.str name];
  17.293 +                    val eq''' = if eq' = "" then eq'' else
  17.294 +                      Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1,
  17.295 +                        eq'', Pretty.brk 1, Pretty.str eq']
  17.296 +                  in
  17.297 +                    ((vs''', old_vals), Pretty.block [pretty_pattern false lhs,
  17.298 +                      Pretty.str " =>",
  17.299 +                      Pretty.brk 1, mk_let "let" 2 (tm_vals @ flat (map (snd o snd) eq_vals))
  17.300 +                        [Pretty.str ("(" ^ ct' ^ ","), Pretty.brk 1, eq''', Pretty.str ")"]])
  17.301 +                  end;
  17.302 +
  17.303 +                val case_names = map (fn i => Sign.base_name cname ^ "_" ^
  17.304 +                  string_of_int i) (1 upto length cases);
  17.305 +                val ((vs3, vals3), case_ps) = foldl_map mk_case_code
  17.306 +                  ((vs2, vals2), us' ~~ (case_names ~~ cases));
  17.307 +                val eq' = variant vs3 "eq";
  17.308 +                val ct' = variant (eq' :: vs3) "ct";
  17.309 +                val eq'' = variant (eq' :: ct' :: vs3) "eq";
  17.310 +                val case_vals =
  17.311 +                  fvals @ ty_vals @
  17.312 +                  [mk_val [ct', eq'] ([Pretty.str "(case", Pretty.brk 1,
  17.313 +                    Pretty.str ("term_of " ^ ct ^ " of"), Pretty.brk 1] @
  17.314 +                    flat (separate [Pretty.brk 1, Pretty.str "| "]
  17.315 +                      (map single case_ps)) @ [Pretty.str ")"])]
  17.316 +              in
  17.317 +                if b then
  17.318 +                  ((eq' :: ct' :: vs3, vals3 @
  17.319 +                     [(t, ((eq', ct'), case_vals))]), (eq', ct'))
  17.320 +                else
  17.321 +                  let val ((vs4, vals4), (_, ctcase)) = mk_simpl_code sg case_tab false
  17.322 +                    fs ts ty_bs thm_bs ((eq' :: eq'' :: ct' :: vs3, vals3), f)
  17.323 +                  in
  17.324 +                    ((vs4, vals4 @ [(t, ((eq'', ct'), case_vals @
  17.325 +                       [mk_val [eq''] [Pretty.str "Thm.transitive", Pretty.brk 1,
  17.326 +                          Pretty.str "(Thm.combination", Pretty.brk 1,
  17.327 +                          Pretty.str "(Thm.reflexive", Pretty.brk 1,
  17.328 +                          mk_apps "Thm.capply" true (Pretty.str ctcase)
  17.329 +                            (map Pretty.str xs),
  17.330 +                          Pretty.str ")", Pretty.brk 1, Pretty.str (eq ^ ")"),
  17.331 +                          Pretty.brk 1, Pretty.str eq']]))]), (eq'', ct'))
  17.332 +                  end
  17.333 +              end
  17.334 +          
  17.335 +          | None =>
  17.336 +            let
  17.337 +              val b = forall (unint sg fs) us;
  17.338 +              val (q, eqs) = foldl_map
  17.339 +                (mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs) ((vs, vals), us);
  17.340 +              val ((vs', vals'), (eqf, ctf)) = if is_some (lookup sg fs f) andalso b
  17.341 +                then (q, ("", ""))
  17.342 +                else mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs (q, f);
  17.343 +              val ct = variant vs' "ct";
  17.344 +              val eq = variant (ct :: vs') "eq";
  17.345 +              val ctv = mk_val [ct] [mk_apps "Thm.capply" false
  17.346 +                (Pretty.str ctf) (map (Pretty.str o snd) eqs)];
  17.347 +              fun combp b = mk_apps "Thm.combination" b
  17.348 +                (Pretty.str eqf) (map (Pretty.str o fst) eqs)
  17.349 +            in
  17.350 +              case (lookup sg fs f, b) of
  17.351 +                (None, true) =>  (* completely uninterpreted *)
  17.352 +                  if mkeq then ((ct :: eq :: vs', vals' @
  17.353 +                    [(t, ((eq, ct), [ctv, mk_refleq eq ct]))]), (eq, ct))
  17.354 +                  else ((ct :: vs', vals' @ [(t, (("", ct), [ctv]))]), ("", ct))
  17.355 +              | (None, false) =>  (* function uninterpreted *)
  17.356 +                  ((eq :: ct :: vs', vals' @
  17.357 +                     [(t, ((eq, ct), [ctv, mk_val [eq] [combp false]]))]), (eq, ct))
  17.358 +              | (Some (s, _, _), true) =>  (* arguments uninterpreted *)
  17.359 +                  ((eq :: ct :: vs', vals' @
  17.360 +                     [(t, ((eq, ct), [mk_val [ct, eq] (separate (Pretty.brk 1)
  17.361 +                       (Pretty.str s :: map (Pretty.str o snd) eqs))]))]), (eq, ct))
  17.362 +              | (Some (s, _, _), false) =>  (* function and arguments interpreted *)
  17.363 +                  let val eq' = variant (eq :: ct :: vs') "eq"
  17.364 +                  in ((eq' :: eq :: ct :: vs', vals' @ [(t, ((eq', ct),
  17.365 +                    [mk_val [ct, eq] (separate (Pretty.brk 1)
  17.366 +                       (Pretty.str s :: map (Pretty.str o snd) eqs)),
  17.367 +                     mk_val [eq'] [Pretty.str "Thm.transitive", Pretty.brk 1,
  17.368 +                       combp true, Pretty.brk 1, Pretty.str eq]]))]), (eq', ct))
  17.369 +                  end
  17.370 +            end
  17.371 +        end));
  17.372 +
  17.373 +fun lhs_of thm = fst (Logic.dest_equals (prop_of thm));
  17.374 +fun rhs_of thm = snd (Logic.dest_equals (prop_of thm));
  17.375 +
  17.376 +fun mk_funs_code sg case_tab fs fs' =
  17.377 +  let
  17.378 +    val case_thms = mapfilter (fn s => (case Symtab.lookup (case_tab, s) of
  17.379 +        None => None
  17.380 +      | Some thms => Some (unsuffix "_case" (Sign.base_name s) ^ ".cases",
  17.381 +          map (fn i => Sign.base_name s ^ "_" ^ string_of_int i)
  17.382 +            (1 upto length thms) ~~ thms)))
  17.383 +      (foldr add_term_consts (map (prop_of o snd)
  17.384 +        (flat (map (#3 o snd) fs')), []));
  17.385 +    val case_vals = map (fn (s, cs) => mk_vall (map fst cs)
  17.386 +      [Pretty.str "map my_mk_meta_eq", Pretty.brk 1,
  17.387 +       Pretty.str ("(thms \"" ^ s ^ "\")")]) case_thms;
  17.388 +    val (vs, thm_bs, thm_vals) = foldl mk_term_bindings (([], [], []),
  17.389 +      flat (map (map (apsnd prop_of) o #3 o snd) fs') @
  17.390 +      map (apsnd prop_of) (flat (map snd case_thms)));
  17.391 +
  17.392 +    fun mk_fun_code (prfx, (fname, d, eqns)) =
  17.393 +      let
  17.394 +        val (f, ts) = strip_comb (lhs_of (snd (hd eqns)));
  17.395 +        val args = variantlist (replicate (length ts) "ct", vs);
  17.396 +        val (vs', ty_bs, ty_vals) = foldl mk_type_bindings
  17.397 +          ((vs @ args, [], []), args ~~ map fastype_of ts);
  17.398 +        val insts1 = map mk_tyinst ty_bs;
  17.399 +
  17.400 +        fun mk_eqn_code (name, eqn) =
  17.401 +          let
  17.402 +            val (_, argts) = strip_comb (lhs_of eqn);
  17.403 +            val (vs'', tm_bs, tm_vals) = foldl (decomp_term_code false)
  17.404 +              ((vs', [], []), args ~~ argts);
  17.405 +            val ((vs''', eq_vals), (eq, ct)) = mk_simpl_code sg case_tab false fs
  17.406 +              (tm_bs @ filter_out (is_Var o fst) thm_bs) ty_bs thm_bs
  17.407 +              ((vs'', []), rhs_of eqn);
  17.408 +            val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(",
  17.409 +              inst_ty false ty_bs t (the (assoc (thm_bs, t))), Pretty.str ",", Pretty.brk 1,
  17.410 +              Pretty.str (s ^ ")")]) tm_bs
  17.411 +            val eq' = if null insts1 andalso null insts2 then Pretty.str name
  17.412 +              else parens (eq <> "") [Pretty.str "Thm.instantiate",
  17.413 +                Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1,
  17.414 +                Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2,
  17.415 +                Pretty.str ")", Pretty.brk 1, Pretty.str name];
  17.416 +            val eq'' = if eq = "" then eq' else
  17.417 +              Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1,
  17.418 +                eq', Pretty.brk 1, Pretty.str eq]
  17.419 +          in
  17.420 +            Pretty.block [parens (length argts > 1)
  17.421 +                (Pretty.commas (map (pretty_pattern false) argts)),
  17.422 +              Pretty.str " =>",
  17.423 +              Pretty.brk 1, mk_let "let" 2 (ty_vals @ tm_vals @ flat (map (snd o snd) eq_vals))
  17.424 +                [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1, eq'', Pretty.str ")"]]
  17.425 +          end;
  17.426 +
  17.427 +        val default = if d then
  17.428 +            let
  17.429 +              val Some s = assoc (thm_bs, f);
  17.430 +              val ct = variant vs' "ct"
  17.431 +            in [Pretty.brk 1, Pretty.str "handle", Pretty.brk 1,
  17.432 +              Pretty.str "Match =>", Pretty.brk 1, mk_let "let" 2
  17.433 +                (ty_vals @ (if null (term_tvars f) then [] else
  17.434 +                   [mk_val [s] [inst_ty false ty_bs f s]]) @
  17.435 +                 [mk_val [ct] [mk_apps "Thm.capply" false (Pretty.str s)
  17.436 +                    (map Pretty.str args)]])
  17.437 +                [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1,
  17.438 +                 Pretty.str "Thm.reflexive", Pretty.brk 1, Pretty.str (ct ^ ")")]]
  17.439 +            end
  17.440 +          else []
  17.441 +      in
  17.442 +        ("and ", Pretty.block (separate (Pretty.brk 1)
  17.443 +            (Pretty.str (prfx ^ fname) :: map Pretty.str args) @
  17.444 +          [Pretty.str " =", Pretty.brk 1, Pretty.str "(case", Pretty.brk 1,
  17.445 +           Pretty.list "(" ")" (map (fn s => Pretty.str ("term_of " ^ s)) args),
  17.446 +           Pretty.str " of", Pretty.brk 1] @
  17.447 +          flat (separate [Pretty.brk 1, Pretty.str "| "]
  17.448 +            (map (single o mk_eqn_code) eqns)) @ [Pretty.str ")"] @ default))
  17.449 +      end;
  17.450 +
  17.451 +    val (_, decls) = foldl_map mk_fun_code ("fun ", map snd fs')
  17.452 +  in
  17.453 +    mk_let "local" 2 (case_vals @ thm_vals) (separate Pretty.fbrk decls)
  17.454 +  end;
  17.455 +
  17.456 +fun mk_simprocs_code sg eqns =
  17.457 +  let
  17.458 +    val case_tab = get_cases sg;
  17.459 +    fun get_head th = head_of (fst (Logic.dest_equals (prop_of th)));
  17.460 +    fun attach_term (x as (_, _, (_, th) :: _)) = (get_head th, x);
  17.461 +    val eqns' = map attach_term eqns;
  17.462 +    fun mk_node (s, _, (_, th) :: _) = (s, get_head th);
  17.463 +    fun mk_edges (s, _, ths) = map (pair s) (distinct
  17.464 +      (mapfilter (fn t => apsome #1 (lookup sg eqns' t))
  17.465 +        (flat (map (term_consts' o prop_of o snd) ths))));
  17.466 +    val gr = foldr (uncurry Graph.add_edge)
  17.467 +      (map (pair "" o #1) eqns @ flat (map mk_edges eqns),
  17.468 +       foldr (uncurry Graph.new_node)
  17.469 +         (("", Bound 0) :: map mk_node eqns, Graph.empty));
  17.470 +    val keys = rev (Graph.all_succs gr [""] \ "");
  17.471 +    fun gr_ord (x :: _, y :: _) =
  17.472 +      int_ord (find_index (equal x) keys, find_index (equal y) keys);
  17.473 +    val scc = map (fn xs => filter (fn (_, (s, _, _)) => s mem xs) eqns')
  17.474 +      (sort gr_ord (Graph.strong_conn gr \ [""]));
  17.475 +  in
  17.476 +    flat (separate [Pretty.str ";", Pretty.fbrk, Pretty.str " ", Pretty.fbrk]
  17.477 +      (map (fn eqns'' => [mk_funs_code sg case_tab eqns' eqns'']) scc)) @
  17.478 +    [Pretty.str ";", Pretty.fbrk]
  17.479 +  end;
  17.480 +
  17.481 +fun use_simprocs_code sg eqns =
  17.482 +  let
  17.483 +    fun attach_name (i, x) = (i+1, ("simp_thm_" ^ string_of_int i, x));
  17.484 +    fun attach_names (i, (s, b, eqs)) =
  17.485 +      let val (i', eqs') = foldl_map attach_name (i, eqs)
  17.486 +      in (i', (s, b, eqs')) end;
  17.487 +    val (_, eqns') = foldl_map attach_names (1, eqns);
  17.488 +    val (names, thms) = split_list (flat (map #3 eqns'));
  17.489 +    val s = setmp print_mode [] Pretty.string_of
  17.490 +      (mk_let "local" 2 [mk_vall names [Pretty.str "!SimprocsCodegen.simp_thms"]]
  17.491 +        (mk_simprocs_code sg eqns'))
  17.492 +  in
  17.493 +    (simp_thms := thms; use_text Context.ml_output false s)
  17.494 +  end;
  17.495 +
  17.496 +end;
    18.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    18.2 +++ b/src/HOL/Matrix/fspmlp.ML	Fri Sep 03 17:10:36 2004 +0200
    18.3 @@ -0,0 +1,303 @@
    18.4 +signature FSPMLP = 
    18.5 +sig
    18.6 +    type linprog
    18.7 +
    18.8 +    val y : linprog -> cterm
    18.9 +    val A : linprog -> cterm * cterm
   18.10 +    val b : linprog -> cterm
   18.11 +    val c : linprog -> cterm * cterm
   18.12 +    val r : linprog -> cterm
   18.13 +
   18.14 +    exception Load of string
   18.15 +		       
   18.16 +    val load : string -> int -> bool -> linprog
   18.17 +end
   18.18 +
   18.19 +structure fspmlp : FSPMLP = 
   18.20 +struct
   18.21 +
   18.22 +type linprog = cterm * (cterm * cterm) * cterm * (cterm * cterm) * cterm 
   18.23 +
   18.24 +fun y (c1, c2, c3, c4, c5) = c1
   18.25 +fun A (c1, c2, c3, c4, c5) = c2
   18.26 +fun b (c1, c2, c3, c4, c5) = c3
   18.27 +fun c (c1, c2, c3, c4, c5) = c4
   18.28 +fun r (c1, c2, c3, c4, c5) = c5
   18.29 +
   18.30 +structure CplexFloatSparseMatrixConverter = 
   18.31 +MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder);
   18.32 +
   18.33 +datatype bound_type = LOWER | UPPER
   18.34 +
   18.35 +fun intbound_ord ((i1, b1),(i2,b2)) = 
   18.36 +    if i1 < i2 then LESS
   18.37 +    else if i1 = i2 then 
   18.38 +	(if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER)
   18.39 +    else GREATER
   18.40 +
   18.41 +structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord));
   18.42 +
   18.43 +structure VarGraph = TableFun(type key = int*bound_type val ord = intbound_ord);
   18.44 +(* key -> (float option) * (int -> (float * (((float * float) * key) list)))) *)
   18.45 +(* dest_key -> (sure_bound * (row_index -> (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *)
   18.46 +
   18.47 +exception Internal of string;
   18.48 +
   18.49 +fun add_row_bound g dest_key row_index row_bound = 
   18.50 +    let 
   18.51 +	val x = 
   18.52 +	    case VarGraph.lookup (g, dest_key) of
   18.53 +		None => (None, Inttab.update ((row_index, (row_bound, [])), Inttab.empty))
   18.54 +	      | Some (sure_bound, f) =>
   18.55 +		(sure_bound,
   18.56 +		 case Inttab.lookup (f, row_index) of
   18.57 +		     None => Inttab.update ((row_index, (row_bound, [])), f)
   18.58 +		   | Some _ => raise (Internal "add_row_bound"))				     
   18.59 +    in
   18.60 +	VarGraph.update ((dest_key, x), g)
   18.61 +    end    
   18.62 +
   18.63 +fun update_sure_bound g (key as (_, btype)) bound = 
   18.64 +    let
   18.65 +	val x = 
   18.66 +	    case VarGraph.lookup (g, key) of
   18.67 +		None => (Some bound, Inttab.empty)
   18.68 +	      | Some (None, f) => (Some bound, f)
   18.69 +	      | Some (Some old_bound, f) => 
   18.70 +		(Some ((case btype of 
   18.71 +			    UPPER => FloatArith.min 
   18.72 +			  | LOWER => FloatArith.max) 
   18.73 +			   old_bound bound), f)
   18.74 +    in
   18.75 +	VarGraph.update ((key, x), g)
   18.76 +    end
   18.77 +
   18.78 +fun get_sure_bound g key = 
   18.79 +    case VarGraph.lookup (g, key) of 
   18.80 +	None => None
   18.81 +      | Some (sure_bound, _) => sure_bound
   18.82 +
   18.83 +(*fun get_row_bound g key row_index = 
   18.84 +    case VarGraph.lookup (g, key) of
   18.85 +	None => None
   18.86 +      | Some (sure_bound, f) =>
   18.87 +	(case Inttab.lookup (f, row_index) of 
   18.88 +	     None => None
   18.89 +	   | Some (row_bound, _) => (sure_bound, row_bound))*)
   18.90 +    
   18.91 +fun add_edge g src_key dest_key row_index coeff = 
   18.92 +    case VarGraph.lookup (g, dest_key) of
   18.93 +	None => raise (Internal "add_edge: dest_key not found")
   18.94 +      | Some (sure_bound, f) =>
   18.95 +	(case Inttab.lookup (f, row_index) of
   18.96 +	     None => raise (Internal "add_edge: row_index not found")
   18.97 +	   | Some (row_bound, sources) => 
   18.98 +	     VarGraph.update ((dest_key, (sure_bound, Inttab.update ((row_index, (row_bound, (coeff, src_key) :: sources)), f))), g))
   18.99 +
  18.100 +fun split_graph g = 
  18.101 +    let
  18.102 +	fun split ((r1, r2), (key, (sure_bound, _))) = 
  18.103 +	    case sure_bound of
  18.104 +		None => (r1, r2)
  18.105 +	      | Some bound => 
  18.106 +		(case key of
  18.107 +		     (u, UPPER) => (r1, Inttab.update ((u, bound), r2))
  18.108 +		   | (u, LOWER) => (Inttab.update ((u, bound), r1), r2))
  18.109 +    in
  18.110 +	VarGraph.foldl split ((Inttab.empty, Inttab.empty), g)
  18.111 +    end
  18.112 +
  18.113 +fun it2list t = 
  18.114 +    let
  18.115 +	fun tolist (l, a) = a::l
  18.116 +    in
  18.117 +	Inttab.foldl tolist ([], t)
  18.118 +    end
  18.119 +
  18.120 +(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm).
  18.121 +   If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *)
  18.122 +fun propagate_sure_bounds safe names g = 
  18.123 +    let		 	    	
  18.124 +	(* returns None if no new sure bound could be calculated, otherwise the new sure bound is returned *)
  18.125 +	fun calc_sure_bound_from_sources g (key as (_, btype)) = 
  18.126 +	    let		
  18.127 +		fun mult_upper x (lower, upper) = 
  18.128 +		    if FloatArith.is_negative x then
  18.129 +			FloatArith.mul x lower
  18.130 +		    else
  18.131 +			FloatArith.mul x upper
  18.132 +			
  18.133 +		fun mult_lower x (lower, upper) = 
  18.134 +		    if FloatArith.is_negative x then
  18.135 +			FloatArith.mul x upper
  18.136 +		    else
  18.137 +			FloatArith.mul x lower
  18.138 +
  18.139 +		val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower
  18.140 +
  18.141 +		fun calc_sure_bound (sure_bound, (row_index, (row_bound, sources))) = 
  18.142 +		    let
  18.143 +			fun add_src_bound (sum, (coeff, src_key)) = 
  18.144 +			    case sum of 
  18.145 +				None => None
  18.146 +			      | Some x => 
  18.147 +				(case get_sure_bound g src_key of
  18.148 +				     None => None
  18.149 +				   | Some src_sure_bound => Some (FloatArith.add x (mult_btype src_sure_bound coeff)))
  18.150 +		    in
  18.151 +			case foldl add_src_bound (Some row_bound, sources) of
  18.152 +			    None => sure_bound
  18.153 +			  | new_sure_bound as (Some new_bound) => 
  18.154 +			    (case sure_bound of 
  18.155 +				 None => new_sure_bound
  18.156 +			       | Some old_bound => 
  18.157 +				 Some (case btype of 
  18.158 +					   UPPER => FloatArith.min old_bound new_bound
  18.159 +					 | LOWER => FloatArith.max old_bound new_bound))				 
  18.160 +		    end		
  18.161 +	    in
  18.162 +		case VarGraph.lookup (g, key) of
  18.163 +		    None => None
  18.164 +		  | Some (sure_bound, f) =>
  18.165 +		    let
  18.166 +			val x = Inttab.foldl calc_sure_bound (sure_bound, f) 
  18.167 +		    in
  18.168 +			if x = sure_bound then None else x
  18.169 +		    end		
  18.170 +    	    end
  18.171 +
  18.172 +	fun propagate ((g, b), (key, _)) = 
  18.173 +	    case calc_sure_bound_from_sources g key of 
  18.174 +		None => (g,b)
  18.175 +	      | Some bound => (update_sure_bound g key bound, 
  18.176 +			       if safe then 
  18.177 +				   case get_sure_bound g key of
  18.178 +				       None => true
  18.179 +				     | _ => b
  18.180 +			       else
  18.181 +				   true)
  18.182 +
  18.183 +	val (g, b) = VarGraph.foldl propagate ((g, false), g) 
  18.184 +    in
  18.185 +	if b then propagate_sure_bounds safe names g else g	
  18.186 +    end	    
  18.187 +    		
  18.188 +exception Load of string;
  18.189 +
  18.190 +fun calcr safe_propagation xlen names prec A b = 
  18.191 +    let
  18.192 +	val empty = Inttab.empty
  18.193 +
  18.194 +	fun instab t i x = Inttab.update ((i,x), t)
  18.195 +
  18.196 +	fun isnegstr x = String.isPrefix "-" x
  18.197 +	fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else "-"^x
  18.198 +
  18.199 +	fun test_1 (lower, upper) = 
  18.200 +	    if lower = upper then
  18.201 +		(if FloatArith.is_equal lower (IntInf.fromInt ~1, FloatArith.izero) then ~1 
  18.202 +		 else if FloatArith.is_equal lower (IntInf.fromInt 1, FloatArith.izero) then 1
  18.203 +		 else 0)
  18.204 +	    else 0	
  18.205 +
  18.206 +	fun calcr (g, (row_index, a)) = 
  18.207 +	    let				
  18.208 +		val b =  FloatSparseMatrixBuilder.v_elem_at b row_index
  18.209 +		val (_, b2) = ExactFloatingPoint.approx_decstr_by_bin prec (case b of None => "0" | Some b => b)
  18.210 +		val approx_a = FloatSparseMatrixBuilder.v_fold (fn (l, (i,s)) => 
  18.211 +								   (i, ExactFloatingPoint.approx_decstr_by_bin prec s)::l) [] a
  18.212 +			       
  18.213 +		fun fold_dest_nodes (g, (dest_index, dest_value)) = 
  18.214 +		    let		
  18.215 +			val dest_test = test_1 dest_value
  18.216 +		    in
  18.217 +			if dest_test = 0 then
  18.218 +			    g
  18.219 +			else let
  18.220 +				val (dest_key as (_, dest_btype), row_bound) = 
  18.221 +				    if dest_test = ~1 then 
  18.222 +					((dest_index, LOWER), FloatArith.neg b2)
  18.223 +				    else
  18.224 +					((dest_index, UPPER), b2)
  18.225 +					
  18.226 +				fun fold_src_nodes (g, (src_index, src_value as (src_lower, src_upper))) = 
  18.227 +				    if src_index = dest_index then g
  18.228 +				    else
  18.229 +					let
  18.230 +					    val coeff = case dest_btype of 
  18.231 +							    UPPER => (FloatArith.neg src_upper, FloatArith.neg src_lower)
  18.232 +							  | LOWER => src_value
  18.233 +					in
  18.234 +					    if FloatArith.is_negative src_lower then
  18.235 +						add_edge g (src_index, UPPER) dest_key row_index coeff
  18.236 +					    else
  18.237 +						add_edge g (src_index, LOWER) dest_key row_index coeff
  18.238 +					end
  18.239 +			    in	    
  18.240 +				foldl fold_src_nodes ((add_row_bound g dest_key row_index row_bound), approx_a)
  18.241 +			    end
  18.242 +		    end
  18.243 +	    in
  18.244 +		case approx_a of
  18.245 +		    [] => g
  18.246 +		  | [(u, a)] => 
  18.247 +		    let
  18.248 +			val atest = test_1 a
  18.249 +		    in
  18.250 +			if atest = ~1 then 			  
  18.251 +			    update_sure_bound g (u, LOWER) (FloatArith.neg b2)
  18.252 +			else if atest = 1 then
  18.253 +			    update_sure_bound g (u, UPPER) b2
  18.254 +			else
  18.255 +			    g
  18.256 +		    end
  18.257 +		  | _ => foldl fold_dest_nodes (g, approx_a)
  18.258 +	    end
  18.259 +	
  18.260 +	val g = FloatSparseMatrixBuilder.m_fold calcr VarGraph.empty A
  18.261 +	val g = propagate_sure_bounds safe_propagation names g
  18.262 +
  18.263 +	val (r1, r2) = split_graph g
  18.264 +	
  18.265 +	fun abs_estimate i r1 r2 = 
  18.266 +	    if i = 0 then FloatSparseMatrixBuilder.empty_spmat
  18.267 +	    else
  18.268 +		let
  18.269 +		    val index = xlen-i
  18.270 +		    val r = abs_estimate (i-1) r1 r2 
  18.271 +		    val b1 = case Inttab.lookup (r1, index) of None => raise (Load ("x-value not bounded from below: "^(names index))) | Some x => x
  18.272 +		    val b2 = case Inttab.lookup (r2, index) of None => raise (Load ("x-value not bounded from above: "^(names index))) | Some x => x
  18.273 +		    val abs_max = FloatArith.max (FloatArith.neg (FloatArith.negative_part b1)) (FloatArith.positive_part b2)    
  18.274 +		    val vec = FloatSparseMatrixBuilder.cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 abs_max) FloatSparseMatrixBuilder.empty_spvec
  18.275 +		in
  18.276 +		    FloatSparseMatrixBuilder.cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) r
  18.277 +		end		    		    
  18.278 +    in
  18.279 +	FloatSparseMatrixBuilder.sign_term (abs_estimate xlen r1 r2)
  18.280 +    end
  18.281 +	    
  18.282 +fun load filename prec safe_propagation =
  18.283 +    let
  18.284 +	val prog = Cplex.load_cplexFile filename
  18.285 +	val prog = Cplex.elim_nonfree_bounds prog
  18.286 +	val prog = Cplex.relax_strict_ineqs prog
  18.287 +	val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog
  18.288 +	val r = calcr safe_propagation xlen names prec A b
  18.289 +	val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems"			
  18.290 +	val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b
  18.291 +	val results = Cplex.solve dualprog
  18.292 +	val (optimal,v) = CplexFloatSparseMatrixConverter.convert_results results indexof
  18.293 +	val A = FloatSparseMatrixBuilder.cut_matrix v None A
  18.294 +	fun id x = x
  18.295 +	val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v
  18.296 +	val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b)
  18.297 +	val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c
  18.298 +	val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec FloatArith.positive_part v
  18.299 +	val A = FloatSparseMatrixBuilder.approx_matrix prec id A
  18.300 +	val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b
  18.301 +	val c = FloatSparseMatrixBuilder.approx_matrix prec id c
  18.302 +    in
  18.303 +	(y1, A, b2, c, r)
  18.304 +    end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s)))
  18.305 +
  18.306 +end
  18.307 \ No newline at end of file
    19.1 --- a/src/HOL/OrderedGroup.thy	Fri Sep 03 10:27:05 2004 +0200
    19.2 +++ b/src/HOL/OrderedGroup.thy	Fri Sep 03 17:10:36 2004 +0200
    19.3 @@ -842,6 +842,42 @@
    19.4  lemma minus_add_cancel: "-(a::'a::ab_group_add) + (a + b) = b"
    19.5  by (simp add: add_assoc[symmetric])
    19.6  
    19.7 +lemma  le_add_right_mono: 
    19.8 +  assumes 
    19.9 +  "a <= b + (c::'a::pordered_ab_group_add)"
   19.10 +  "c <= d"    
   19.11 +  shows "a <= b + d"
   19.12 +  apply (rule_tac order_trans[where y = "b+c"])
   19.13 +  apply (simp_all add: prems)
   19.14 +  done
   19.15 +
   19.16 +lemmas group_eq_simps =
   19.17 +  mult_ac
   19.18 +  add_ac
   19.19 +  add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2
   19.20 +  diff_eq_eq eq_diff_eq
   19.21 +
   19.22 +lemma estimate_by_abs:
   19.23 +"a + b <= (c::'a::lordered_ab_group_abs) \<Longrightarrow> a <= c + abs b" 
   19.24 +proof -
   19.25 +  assume 1: "a+b <= c"
   19.26 +  have 2: "a <= c+(-b)"
   19.27 +    apply (insert 1)
   19.28 +    apply (drule_tac add_right_mono[where c="-b"])
   19.29 +    apply (simp add: group_eq_simps)
   19.30 +    done
   19.31 +  have 3: "(-b) <= abs b" by (rule abs_ge_minus_self)
   19.32 +  show ?thesis by (rule le_add_right_mono[OF 2 3])
   19.33 +qed
   19.34 +
   19.35 +lemma abs_of_ge_0: "0 <= (y::'a::lordered_ab_group_abs) \<Longrightarrow> abs y = y"
   19.36 +proof -
   19.37 +  assume 1:"0 <= y"
   19.38 +  have 2:"-y <= 0" by (simp add: 1)
   19.39 +  from 1 2 have 3:"-y <= y" by (simp only:)
   19.40 +  show ?thesis by (simp add: abs_lattice ge_imp_join_eq[OF 3])
   19.41 +qed
   19.42 +
   19.43  ML {*
   19.44  val add_zero_left = thm"add_0";
   19.45  val add_zero_right = thm"add_0_right";
    20.1 --- a/src/HOL/Ring_and_Field.thy	Fri Sep 03 10:27:05 2004 +0200
    20.2 +++ b/src/HOL/Ring_and_Field.thy	Fri Sep 03 17:10:36 2004 +0200
    20.3 @@ -546,14 +546,14 @@
    20.4  done
    20.5  
    20.6  text{*This list of rewrites decides ring equalities by ordered rewriting.*}
    20.7 -lemmas ring_eq_simps =
    20.8 -  mult_ac
    20.9 +lemmas ring_eq_simps =  
   20.10 +(*  mult_ac*)
   20.11    left_distrib right_distrib left_diff_distrib right_diff_distrib
   20.12 -  add_ac
   20.13 +  group_eq_simps
   20.14 +(*  add_ac
   20.15    add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2
   20.16 -  diff_eq_eq eq_diff_eq
   20.17 +  diff_eq_eq eq_diff_eq *)
   20.18      
   20.19 -
   20.20  subsection {* Fields *}
   20.21  
   20.22  lemma right_inverse [simp]:
   20.23 @@ -1503,6 +1503,90 @@
   20.24  text{*Moving this up spoils many proofs using @{text mult_le_cancel_right}*}
   20.25  declare times_divide_eq_left [simp]
   20.26  
   20.27 +lemma linprog_dual_estimate:
   20.28 +  assumes
   20.29 +  "A * x \<le> (b::'a::lordered_ring)"
   20.30 +  "0 \<le> y"
   20.31 +  "abs (A - A') \<le> \<delta>A"
   20.32 +  "b \<le> b'"
   20.33 +  "abs (c - c') \<le> \<delta>c"
   20.34 +  "abs x \<le> r"
   20.35 +  shows
   20.36 +  "c * x \<le> y * b' + (y * \<delta>A + abs (y * A' - c') + \<delta>c) * r"
   20.37 +proof -
   20.38 +  from prems have 1: "y * b <= y * b'" by (simp add: mult_left_mono)
   20.39 +  from prems have 2: "y * (A * x) <= y * b" by (simp add: mult_left_mono) 
   20.40 +  have 3: "y * (A * x) = c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x" by (simp add: ring_eq_simps)  
   20.41 +  from 1 2 3 have 4: "c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x <= y * b'" by simp
   20.42 +  have 5: "c * x <= y * b' + abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)"
   20.43 +    by (simp only: 4 estimate_by_abs)  
   20.44 +  have 6: "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x) <= abs (y * (A - A') + (y * A' - c') + (c'-c)) * abs x"
   20.45 +    by (simp add: abs_le_mult)
   20.46 +  have 7: "(abs (y * (A - A') + (y * A' - c') + (c'-c))) * abs x <= (abs (y * (A-A') + (y*A'-c')) + abs(c'-c)) * abs x"
   20.47 +    by (simp add: abs_triangle_ineq mult_right_mono)
   20.48 +  have 8: " (abs (y * (A-A') + (y*A'-c')) + abs(c'-c)) * abs x <=  (abs (y * (A-A')) + abs (y*A'-c') + abs(c'-c)) * abs x"
   20.49 +    by (simp add: abs_triangle_ineq mult_right_mono)    
   20.50 +  have 9: "(abs (y * (A-A')) + abs (y*A'-c') + abs(c'-c)) * abs x <= (abs y * abs (A-A') + abs (y*A'-c') + abs (c'-c)) * abs x"
   20.51 +    by (simp add: abs_le_mult mult_right_mono)  
   20.52 +  have 10: "c'-c = -(c-c')" by (simp add: ring_eq_simps)
   20.53 +  have 11: "abs (c'-c) = abs (c-c')" 
   20.54 +    by (subst 10, subst abs_minus_cancel, simp)
   20.55 +  have 12: "(abs y * abs (A-A') + abs (y*A'-c') + abs (c'-c)) * abs x <= (abs y * abs (A-A') + abs (y*A'-c') + \<delta>c) * abs x"
   20.56 +    by (simp add: 11 prems mult_right_mono)
   20.57 +  have 13: "(abs y * abs (A-A') + abs (y*A'-c') + \<delta>c) * abs x <= (abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * abs x"
   20.58 +    by (simp add: prems mult_right_mono mult_left_mono)  
   20.59 +  have r: "(abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * abs x <=  (abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * r"
   20.60 +    apply (rule mult_left_mono)
   20.61 +    apply (simp add: prems)
   20.62 +    apply (rule_tac add_mono[of "0::'a" _ "0", simplified])+
   20.63 +    apply (rule mult_left_mono[of "0" "\<delta>A", simplified])
   20.64 +    apply (simp_all)
   20.65 +    apply (rule order_trans[where y="abs (A-A')"], simp_all add: prems)
   20.66 +    apply (rule order_trans[where y="abs (c-c')"], simp_all add: prems)
   20.67 +    done    
   20.68 +  from 6 7 8 9 12 13 r have 14:" abs((y * (A - A') + (y * A' - c') + (c'-c)) * x) <=(abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * r"     
   20.69 +    by (simp)
   20.70 +  show ?thesis 
   20.71 +    apply (rule_tac le_add_right_mono[of _ _ "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)"])
   20.72 +    apply (simp_all add: 5 14[simplified abs_of_ge_0[of y, simplified prems]])
   20.73 +    done
   20.74 +qed
   20.75 +
   20.76 +lemma le_ge_imp_abs_diff_1:
   20.77 +  assumes
   20.78 +  "A1 <= (A::'a::lordered_ring)"
   20.79 +  "A <= A2" 
   20.80 +  shows "abs (A-A1) <= A2-A1"
   20.81 +proof -
   20.82 +  have "0 <= A - A1"    
   20.83 +  proof -
   20.84 +    have 1: "A - A1 = A + (- A1)" by simp
   20.85 +    show ?thesis by (simp only: 1 add_right_mono[of A1 A "-A1", simplified, simplified prems])
   20.86 +  qed
   20.87 +  then have "abs (A-A1) = A-A1" by (rule abs_of_ge_0)
   20.88 +  with prems show "abs (A-A1) <= (A2-A1)" by simp
   20.89 +qed
   20.90 +
   20.91 +lemma linprog_dual_estimate_1:
   20.92 +  assumes
   20.93 +  "A * x \<le> (b::'a::lordered_ring)"
   20.94 +  "0 \<le> y"
   20.95 +  "A1 <= A"
   20.96 +  "A <= A2"
   20.97 +  "c1 <= c"
   20.98 +  "c <= c2"
   20.99 +  "abs x \<le> r"
  20.100 +  shows
  20.101 +  "c * x \<le> y * b + (y * (A2 - A1) + abs (y * A1 - c1) + (c2 - c1)) * r"
  20.102 +proof -
  20.103 +  from prems have delta_A: "abs (A-A1) <= (A2-A1)" by (simp add: le_ge_imp_abs_diff_1)
  20.104 +  from prems have delta_c: "abs (c-c1) <= (c2-c1)" by (simp add: le_ge_imp_abs_diff_1)
  20.105 +  show ?thesis
  20.106 +    apply (rule_tac linprog_dual_estimate)
  20.107 +    apply (auto intro: delta_A delta_c simp add: prems)
  20.108 +    done
  20.109 +qed
  20.110 +
  20.111  ML {*
  20.112  val left_distrib = thm "left_distrib";
  20.113  val right_distrib = thm "right_distrib";