--- a/src/HOL/IsaMakefile Fri Sep 03 10:27:05 2004 +0200
+++ b/src/HOL/IsaMakefile Fri Sep 03 17:10:36 2004 +0200
@@ -1,5 +1,3 @@
-#
-# $Id$
#
# IsaMakefile for HOL
#
@@ -80,8 +78,7 @@
$(SRC)/Provers/quasi.ML \
$(SRC)/Provers/simplifier.ML $(SRC)/Provers/splitter.ML \
$(SRC)/Provers/trancl.ML \
- $(SRC)/TFL/isand.ML $(SRC)/TFL/casesplit.ML $(SRC)/TFL/dcterm.ML\
- $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \
+ $(SRC)/TFL/dcterm.ML $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \
$(SRC)/TFL/tfl.ML $(SRC)/TFL/thms.ML $(SRC)/TFL/thry.ML \
$(SRC)/TFL/usyntax.ML $(SRC)/TFL/utils.ML \
Datatype.thy Datatype_Universe.ML Datatype_Universe.thy \
@@ -97,8 +94,7 @@
Nat.thy NatArith.ML NatArith.thy \
Power.thy PreList.thy Product_Type.ML Product_Type.thy \
Refute.thy ROOT.ML \
- Recdef.thy Reconstruction.thy Record.thy\
- Relation.ML Relation.thy Relation_Power.ML \
+ Recdef.thy Record.thy Relation.ML Relation.thy Relation_Power.ML \
Relation_Power.thy LOrder.thy OrderedGroup.thy OrderedGroup.ML Ring_and_Field.thy\
Set.ML Set.thy SetInterval.ML SetInterval.thy \
Sum_Type.ML Sum_Type.thy Tools/datatype_abs_proofs.ML Tools/datatype_aux.ML \
@@ -109,7 +105,7 @@
Tools/primrec_package.ML \
Tools/prop_logic.ML \
Tools/recdef_package.ML Tools/recfun_codegen.ML \
- Tools/record_package.ML Tools/reconstruction.ML\
+ Tools/record_package.ML \
Tools/refute.ML Tools/refute_isar.ML \
Tools/rewrite_hol_proof.ML \
Tools/sat_solver.ML \
@@ -635,14 +631,18 @@
## HOL-Matrix
-HOL-Matrix: HOL $(LOG)/HOL-Matrix.gz
+HOL-Matrix: HOL HOL-Complex $(LOG)/HOL-Matrix.gz
-$(LOG)/HOL-Matrix.gz: $(OUT)/HOL Matrix/ROOT.ML \
- Matrix/Matrix.thy Matrix/LinProg.thy Matrix/MatrixGeneral.thy \
- Matrix/SparseMatrix.thy Matrix/document/root.tex
- @$(ISATOOL) usedir $(OUT)/HOL Matrix
-
-
+$(LOG)/HOL-Matrix.gz: $(OUT)/HOL-Complex \
+ Matrix/MatrixGeneral.thy Matrix/Matrix.thy Matrix/SparseMatrix.thy \
+ Matrix/Float.thy Matrix/FloatArith.ML Matrix/ExactFloatingPoint.ML \
+ Matrix/Cplex.ML Matrix/CplexMatrixConverter.ML \
+ Matrix/FloatSparseMatrixBuilder.ML \
+ Matrix/conv.ML Matrix/eq_codegen.ML Matrix/codegen_prep.ML \
+ Matrix/fspmlp.ML \
+ Matrix/MatrixLP_gensimp.ML Matrix/MatrixLP.ML Matrix/MatrixLP.thy \
+ Matrix/document/root.tex Matrix/ROOT.ML
+ @cd Matrix; $(ISATOOL) usedir -b $(OUT)/HOL-Complex HOL-Matrix
## TLA
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/Cplex.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,1027 @@
+signature CPLEX =
+sig
+
+ datatype cplexTerm = cplexVar of string | cplexNum of string | cplexInf
+ | cplexNeg of cplexTerm
+ | cplexProd of cplexTerm * cplexTerm
+ | cplexSum of (cplexTerm list)
+
+ datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq
+
+ datatype cplexGoal = cplexMinimize of cplexTerm
+ | cplexMaximize of cplexTerm
+
+ datatype cplexConstr = cplexConstr of cplexComp *
+ (cplexTerm * cplexTerm)
+
+ datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm
+ * cplexComp * cplexTerm
+ | cplexBound of cplexTerm * cplexComp * cplexTerm
+
+ datatype cplexProg = cplexProg of string
+ * cplexGoal
+ * ((string option * cplexConstr)
+ list)
+ * cplexBounds list
+
+ datatype cplexResult = Unbounded
+ | Infeasible
+ | Undefined
+ | Optimal of string *
+ (((* name *) string *
+ (* value *) string) list)
+
+ exception Load_cplexFile of string
+ exception Load_cplexResult of string
+ exception Save_cplexFile of string
+ exception Execute of string
+
+ val load_cplexFile : string -> cplexProg
+
+ val save_cplexFile : string -> cplexProg -> unit
+
+ val elim_nonfree_bounds : cplexProg -> cplexProg
+
+ val relax_strict_ineqs : cplexProg -> cplexProg
+
+ val is_normed_cplexProg : cplexProg -> bool
+
+ val load_cplexResult : string -> cplexResult
+
+ val solve : cplexProg -> cplexResult
+end;
+
+structure Cplex : CPLEX =
+struct
+
+exception Load_cplexFile of string;
+exception Load_cplexResult of string;
+exception Save_cplexFile of string;
+
+datatype cplexTerm = cplexVar of string
+ | cplexNum of string
+ | cplexInf
+ | cplexNeg of cplexTerm
+ | cplexProd of cplexTerm * cplexTerm
+ | cplexSum of (cplexTerm list)
+datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq
+datatype cplexGoal = cplexMinimize of cplexTerm | cplexMaximize of cplexTerm
+datatype cplexConstr = cplexConstr of cplexComp * (cplexTerm * cplexTerm)
+datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm
+ * cplexComp * cplexTerm
+ | cplexBound of cplexTerm * cplexComp * cplexTerm
+datatype cplexProg = cplexProg of string
+ * cplexGoal
+ * ((string option * cplexConstr) list)
+ * cplexBounds list
+
+fun rev_cmp cplexLe = cplexGe
+ | rev_cmp cplexLeq = cplexGeq
+ | rev_cmp cplexGe = cplexLe
+ | rev_cmp cplexGeq = cplexLeq
+ | rev_cmp cplexEq = cplexEq
+
+fun the None = raise (Load_cplexFile "Some expected")
+ | the (Some x) = x;
+
+fun modulo_signed is_something (cplexNeg u) = is_something u
+ | modulo_signed is_something u = is_something u
+
+fun is_Num (cplexNum _) = true
+ | is_Num _ = false
+
+fun is_Inf cplexInf = true
+ | is_Inf _ = false
+
+fun is_Var (cplexVar _) = true
+ | is_Var _ = false
+
+fun is_Neg (cplexNeg x ) = true
+ | is_Neg _ = false
+
+fun is_normed_Prod (cplexProd (t1, t2)) =
+ (is_Num t1) andalso (is_Var t2)
+ | is_normed_Prod x = is_Var x
+
+fun is_normed_Sum (cplexSum ts) =
+ (ts <> []) andalso forall (modulo_signed is_normed_Prod) ts
+ | is_normed_Sum x = modulo_signed is_normed_Prod x
+
+fun is_normed_Constr (cplexConstr (c, (t1, t2))) =
+ (is_normed_Sum t1) andalso (modulo_signed is_Num t2)
+
+fun is_Num_or_Inf x = is_Inf x orelse is_Num x
+
+fun is_normed_Bounds (cplexBounds (t1, c1, t2, c2, t3)) =
+ (c1 = cplexLe orelse c1 = cplexLeq) andalso
+ (c2 = cplexLe orelse c2 = cplexLeq) andalso
+ is_Var t2 andalso
+ modulo_signed is_Num_or_Inf t1 andalso
+ modulo_signed is_Num_or_Inf t3
+ | is_normed_Bounds (cplexBound (t1, c, t2)) =
+ (is_Var t1 andalso (modulo_signed is_Num_or_Inf t2))
+ orelse
+ (c <> cplexEq andalso
+ is_Var t2 andalso (modulo_signed is_Num_or_Inf t1))
+
+fun term_of_goal (cplexMinimize x) = x
+ | term_of_goal (cplexMaximize x) = x
+
+fun is_normed_cplexProg (cplexProg (name, goal, constraints, bounds)) =
+ is_normed_Sum (term_of_goal goal) andalso
+ forall (fn (_,x) => is_normed_Constr x) constraints andalso
+ forall is_normed_Bounds bounds
+
+fun is_NL s = s = "\n"
+
+fun is_blank s = forall (fn c => c <> #"\n" andalso Char.isSpace c) (String.explode s)
+
+fun is_num a =
+ let
+ val b = String.explode a
+ fun num4 cs = forall Char.isDigit cs
+ fun num3 [] = true
+ | num3 (ds as (c::cs)) =
+ if c = #"+" orelse c = #"-" then
+ num4 cs
+ else
+ num4 ds
+ fun num2 [] = true
+ | num2 (c::cs) =
+ if c = #"e" orelse c = #"E" then num3 cs
+ else (Char.isDigit c) andalso num2 cs
+ fun num1 [] = true
+ | num1 (c::cs) =
+ if c = #"." then num2 cs
+ else if c = #"e" orelse c = #"E" then num3 cs
+ else (Char.isDigit c) andalso num1 cs
+ fun num [] = true
+ | num (c::cs) =
+ if c = #"." then num2 cs
+ else (Char.isDigit c) andalso num1 cs
+ in
+ num b
+ end
+
+fun is_delimiter s = s = "+" orelse s = "-" orelse s = ":"
+
+fun is_cmp s = s = "<" orelse s = ">" orelse s = "<="
+ orelse s = ">=" orelse s = "="
+
+fun is_symbol a =
+ let
+ val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~"
+ fun is_symbol_char c = Char.isAlphaNum c orelse
+ exists (fn d => d=c) symbol_char
+ fun is_symbol_start c = is_symbol_char c andalso
+ not (Char.isDigit c) andalso
+ not (c= #".")
+ val b = String.explode a
+ in
+ b <> [] andalso is_symbol_start (hd b) andalso
+ forall is_symbol_char b
+ end
+
+fun to_upper s = String.implode (map Char.toUpper (String.explode s))
+
+fun keyword x =
+ let
+ val a = to_upper x
+ in
+ if a = "BOUNDS" orelse a = "BOUND" then
+ Some "BOUNDS"
+ else if a = "MINIMIZE" orelse a = "MINIMUM" orelse a = "MIN" then
+ Some "MINIMIZE"
+ else if a = "MAXIMIZE" orelse a = "MAXIMUM" orelse a = "MAX" then
+ Some "MAXIMIZE"
+ else if a = "ST" orelse a = "S.T." orelse a = "ST." then
+ Some "ST"
+ else if a = "FREE" orelse a = "END" then
+ Some a
+ else if a = "GENERAL" orelse a = "GENERALS" orelse a = "GEN" then
+ Some "GENERAL"
+ else if a = "INTEGER" orelse a = "INTEGERS" orelse a = "INT" then
+ Some "INTEGER"
+ else if a = "BINARY" orelse a = "BINARIES" orelse a = "BIN" then
+ Some "BINARY"
+ else if a = "INF" orelse a = "INFINITY" then
+ Some "INF"
+ else
+ None
+ end
+
+val TOKEN_ERROR = ~1
+val TOKEN_BLANK = 0
+val TOKEN_NUM = 1
+val TOKEN_DELIMITER = 2
+val TOKEN_SYMBOL = 3
+val TOKEN_LABEL = 4
+val TOKEN_CMP = 5
+val TOKEN_KEYWORD = 6
+val TOKEN_NL = 7
+
+(* tokenize takes a list of chars as argument and returns a list of
+ int * string pairs, each string representing a "cplex token",
+ and each int being one of TOKEN_NUM, TOKEN_DELIMITER, TOKEN_CMP
+ or TOKEN_SYMBOL *)
+fun tokenize s =
+ let
+ val flist = [(is_NL, TOKEN_NL),
+ (is_blank, TOKEN_BLANK),
+ (is_num, TOKEN_NUM),
+ (is_delimiter, TOKEN_DELIMITER),
+ (is_cmp, TOKEN_CMP),
+ (is_symbol, TOKEN_SYMBOL)]
+ fun match_helper [] s = (fn x => false, TOKEN_ERROR)
+ | match_helper (f::fs) s =
+ if ((fst f) s) then f else match_helper fs s
+ fun match s = match_helper flist s
+ fun tok s =
+ if s = "" then [] else
+ let
+ val h = String.substring (s,0,1)
+ val (f, j) = match h
+ fun len i =
+ if size s = i then i
+ else if f (String.substring (s,0,i+1)) then
+ len (i+1)
+ else i
+ in
+ if j < 0 then
+ (if h = "\\" then []
+ else raise (Load_cplexFile ("token expected, found: "
+ ^s)))
+ else
+ let
+ val l = len 1
+ val u = String.substring (s,0,l)
+ val v = String.extract (s,l,NONE)
+ in
+ if j = 0 then tok v else (j, u) :: tok v
+ end
+ end
+ in
+ tok s
+ end
+
+exception Tokenize of string;
+
+fun tokenize_general flist s =
+ let
+ fun match_helper [] s = raise (Tokenize s)
+ | match_helper (f::fs) s =
+ if ((fst f) s) then f else match_helper fs s
+ fun match s = match_helper flist s
+ fun tok s =
+ if s = "" then [] else
+ let
+ val h = String.substring (s,0,1)
+ val (f, j) = match h
+ fun len i =
+ if size s = i then i
+ else if f (String.substring (s,0,i+1)) then
+ len (i+1)
+ else i
+ val l = len 1
+ in
+ (j, String.substring (s,0,l)) :: tok (String.extract (s,l,NONE))
+ end
+ in
+ tok s
+ end
+
+fun load_cplexFile name =
+ let
+ val f = TextIO.openIn name
+ val ignore_NL = ref true
+ val rest = ref []
+
+ fun is_symbol s c = (fst c) = TOKEN_SYMBOL andalso (to_upper (snd c)) = s
+
+ fun readToken_helper () =
+ if length (!rest) > 0 then
+ let val u = hd (!rest) in
+ (
+ rest := tl (!rest);
+ Some u
+ )
+ end
+ else
+ let val s = TextIO.inputLine f in
+ if s = "" then None else
+ let val t = tokenize s in
+ if (length t >= 2 andalso
+ snd(hd (tl t)) = ":")
+ then
+ rest := (TOKEN_LABEL, snd (hd t)) :: (tl (tl t))
+ else if (length t >= 2) andalso is_symbol "SUBJECT" (hd (t))
+ andalso is_symbol "TO" (hd (tl t))
+ then
+ rest := (TOKEN_SYMBOL, "ST") :: (tl (tl t))
+ else
+ rest := t;
+ readToken_helper ()
+ end
+ end
+
+ fun readToken_helper2 () =
+ let val c = readToken_helper () in
+ if c = None then None
+ else if !ignore_NL andalso fst (the c) = TOKEN_NL then
+ readToken_helper2 ()
+ else if fst (the c) = TOKEN_SYMBOL
+ andalso keyword (snd (the c)) <> None
+ then Some (TOKEN_KEYWORD, the (keyword (snd (the c))))
+ else c
+ end
+
+ fun readToken () = readToken_helper2 ()
+
+ fun pushToken a = rest := (a::(!rest))
+
+ fun is_value token =
+ fst token = TOKEN_NUM orelse (fst token = TOKEN_KEYWORD
+ andalso snd token = "INF")
+
+ fun get_value token =
+ if fst token = TOKEN_NUM then
+ cplexNum (snd token)
+ else if fst token = TOKEN_KEYWORD andalso snd token = "INF"
+ then
+ cplexInf
+ else
+ raise (Load_cplexFile "num expected")
+
+ fun readTerm_Product only_num =
+ let val c = readToken () in
+ if c = None then None
+ else if fst (the c) = TOKEN_SYMBOL
+ then (
+ if only_num then (pushToken (the c); None)
+ else Some (cplexVar (snd (the c)))
+ )
+ else if only_num andalso is_value (the c) then
+ Some (get_value (the c))
+ else if is_value (the c) then
+ let val t1 = get_value (the c)
+ val d = readToken ()
+ in
+ if d = None then Some t1
+ else if fst (the d) = TOKEN_SYMBOL then
+ Some (cplexProd (t1, cplexVar (snd (the d))))
+ else
+ (pushToken (the d); Some t1)
+ end
+ else (pushToken (the c); None)
+ end
+
+ fun readTerm_Signed only_signed only_num =
+ let
+ val c = readToken ()
+ in
+ if c = None then None
+ else
+ let val d = the c in
+ if d = (TOKEN_DELIMITER, "+") then
+ readTerm_Product only_num
+ else if d = (TOKEN_DELIMITER, "-") then
+ Some (cplexNeg (the (readTerm_Product
+ only_num)))
+ else (pushToken d;
+ if only_signed then None
+ else readTerm_Product only_num)
+ end
+ end
+
+ fun readTerm_Sum first_signed =
+ let val c = readTerm_Signed first_signed false in
+ if c = None then [] else (the c)::(readTerm_Sum true)
+ end
+
+ fun readTerm () =
+ let val c = readTerm_Sum false in
+ if c = [] then None
+ else if tl c = [] then Some (hd c)
+ else Some (cplexSum c)
+ end
+
+ fun readLabeledTerm () =
+ let val c = readToken () in
+ if c = None then (None, None)
+ else if fst (the c) = TOKEN_LABEL then
+ let val t = readTerm () in
+ if t = None then
+ raise (Load_cplexFile ("term after label "^
+ (snd (the c))^
+ " expected"))
+ else (Some (snd (the c)), t)
+ end
+ else (pushToken (the c); (None, readTerm ()))
+ end
+
+ fun readGoal () =
+ let
+ val g = readToken ()
+ in
+ if g = Some (TOKEN_KEYWORD, "MAXIMIZE") then
+ cplexMaximize (the (snd (readLabeledTerm ())))
+ else if g = Some (TOKEN_KEYWORD, "MINIMIZE") then
+ cplexMinimize (the (snd (readLabeledTerm ())))
+ else raise (Load_cplexFile "MAXIMIZE or MINIMIZE expected")
+ end
+
+ fun str2cmp b =
+ (case b of
+ "<" => cplexLe
+ | "<=" => cplexLeq
+ | ">" => cplexGe
+ | ">=" => cplexGeq
+ | "=" => cplexEq
+ | _ => raise (Load_cplexFile (b^" is no TOKEN_CMP")))
+
+ fun readConstraint () =
+ let
+ val t = readLabeledTerm ()
+ fun make_constraint b t1 t2 =
+ cplexConstr
+ (str2cmp b,
+ (t1, t2))
+ in
+ if snd t = None then None
+ else
+ let val c = readToken () in
+ if c = None orelse fst (the c) <> TOKEN_CMP
+ then raise (Load_cplexFile "TOKEN_CMP expected")
+ else
+ let val n = readTerm_Signed false true in
+ if n = None then
+ raise (Load_cplexFile "num expected")
+ else
+ Some (fst t,
+ make_constraint (snd (the c))
+ (the (snd t))
+ (the n))
+ end
+ end
+ end
+
+ fun readST () =
+ let
+ fun readbody () =
+ let val t = readConstraint () in
+ if t = None then []
+ else if (is_normed_Constr (snd (the t))) then
+ (the t)::(readbody ())
+ else if (fst (the t) <> None) then
+ raise (Load_cplexFile
+ ("constraint '"^(the (fst (the t)))^
+ "'is not normed"))
+ else
+ raise (Load_cplexFile
+ "constraint is not normed")
+ end
+ in
+ if readToken () = Some (TOKEN_KEYWORD, "ST")
+ then
+ readbody ()
+ else
+ raise (Load_cplexFile "ST expected")
+ end
+
+ fun readCmp () =
+ let val c = readToken () in
+ if c = None then None
+ else if fst (the c) = TOKEN_CMP then
+ Some (str2cmp (snd (the c)))
+ else (pushToken (the c); None)
+ end
+
+ fun skip_NL () =
+ let val c = readToken () in
+ if c <> None andalso fst (the c) = TOKEN_NL then
+ skip_NL ()
+ else
+ (pushToken (the c); ())
+ end
+
+ fun is_var (cplexVar _) = true
+ | is_var _ = false
+
+ fun make_bounds c t1 t2 =
+ cplexBound (t1, c, t2)
+
+ fun readBound () =
+ let
+ val _ = skip_NL ()
+ val t1 = readTerm ()
+ in
+ if t1 = None then None
+ else
+ let
+ val c1 = readCmp ()
+ in
+ if c1 = None then
+ let
+ val c = readToken ()
+ in
+ if c = Some (TOKEN_KEYWORD, "FREE") then
+ Some (
+ cplexBounds (cplexNeg cplexInf,
+ cplexLeq,
+ the t1,
+ cplexLeq,
+ cplexInf))
+ else
+ raise (Load_cplexFile "FREE expected")
+ end
+ else
+ let
+ val t2 = readTerm ()
+ in
+ if t2 = None then
+ raise (Load_cplexFile "term expected")
+ else
+ let val c2 = readCmp () in
+ if c2 = None then
+ Some (make_bounds (the c1)
+ (the t1)
+ (the t2))
+ else
+ Some (
+ cplexBounds (the t1,
+ the c1,
+ the t2,
+ the c2,
+ the (readTerm())))
+ end
+ end
+ end
+ end
+
+ fun readBounds () =
+ let
+ fun makestring b = "?"
+ fun readbody () =
+ let
+ val b = readBound ()
+ in
+ if b = None then []
+ else if (is_normed_Bounds (the b)) then
+ (the b)::(readbody())
+ else (
+ raise (Load_cplexFile
+ ("bounds are not normed in: "^
+ (makestring (the b)))))
+ end
+ in
+ if readToken () = Some (TOKEN_KEYWORD, "BOUNDS") then
+ readbody ()
+ else raise (Load_cplexFile "BOUNDS expected")
+ end
+
+ fun readEnd () =
+ if readToken () = Some (TOKEN_KEYWORD, "END") then ()
+ else raise (Load_cplexFile "END expected")
+
+ val result_Goal = readGoal ()
+ val result_ST = readST ()
+ val _ = ignore_NL := false
+ val result_Bounds = readBounds ()
+ val _ = ignore_NL := true
+ val _ = readEnd ()
+ val _ = TextIO.closeIn f
+ in
+ cplexProg (name, result_Goal, result_ST, result_Bounds)
+ end
+
+fun save_cplexFile filename (cplexProg (name, goal, constraints, bounds)) =
+ let
+ val f = TextIO.openOut filename
+
+ fun basic_write s = TextIO.output(f, s)
+
+ val linebuf = ref ""
+ fun buf_flushline s =
+ (basic_write (!linebuf);
+ basic_write "\n";
+ linebuf := s)
+ fun buf_add s = linebuf := (!linebuf) ^ s
+
+ fun write s =
+ if (String.size s) + (String.size (!linebuf)) >= 250 then
+ buf_flushline (" "^s)
+ else
+ buf_add s
+
+ fun writeln s = (buf_add s; buf_flushline "")
+
+ fun write_term (cplexVar x) = write x
+ | write_term (cplexNum x) = write x
+ | write_term cplexInf = write "inf"
+ | write_term (cplexProd (cplexNum "1", b)) = write_term b
+ | write_term (cplexProd (a, b)) =
+ (write_term a; write " "; write_term b)
+ | write_term (cplexNeg x) = (write " - "; write_term x)
+ | write_term (cplexSum ts) = write_terms ts
+ and write_terms [] = ()
+ | write_terms (t::ts) =
+ (if (not (is_Neg t)) then write " + " else ();
+ write_term t; write_terms ts)
+
+ fun write_goal (cplexMaximize term) =
+ (writeln "MAXIMIZE"; write_term term; writeln "")
+ | write_goal (cplexMinimize term) =
+ (writeln "MINIMIZE"; write_term term; writeln "")
+
+ fun write_cmp cplexLe = write "<"
+ | write_cmp cplexLeq = write "<="
+ | write_cmp cplexEq = write "="
+ | write_cmp cplexGe = write ">"
+ | write_cmp cplexGeq = write ">="
+
+ fun write_constr (cplexConstr (cmp, (a,b))) =
+ (write_term a;
+ write " ";
+ write_cmp cmp;
+ write " ";
+ write_term b)
+
+ fun write_constraints [] = ()
+ | write_constraints (c::cs) =
+ (if (fst c <> None)
+ then
+ (write (the (fst c)); write ": ")
+ else
+ ();
+ write_constr (snd c);
+ writeln "";
+ write_constraints cs)
+
+ fun write_bounds [] = ()
+ | write_bounds ((cplexBounds (t1,c1,t2,c2,t3))::bs) =
+ ((if t1 = cplexNeg cplexInf andalso t3 = cplexInf
+ andalso (c1 = cplexLeq orelse c1 = cplexLe)
+ andalso (c2 = cplexLeq orelse c2 = cplexLe)
+ then
+ (write_term t2; write " free")
+ else
+ (write_term t1; write " "; write_cmp c1; write " ";
+ write_term t2; write " "; write_cmp c2; write " ";
+ write_term t3)
+ ); writeln ""; write_bounds bs)
+ | write_bounds ((cplexBound (t1, c, t2)) :: bs) =
+ (write_term t1; write " ";
+ write_cmp c; write " ";
+ write_term t2; writeln ""; write_bounds bs)
+
+ val _ = write_goal goal
+ val _ = (writeln ""; writeln "ST")
+ val _ = write_constraints constraints
+ val _ = (writeln ""; writeln "BOUNDS")
+ val _ = write_bounds bounds
+ val _ = (writeln ""; writeln "END")
+ val _ = TextIO.closeOut f
+ in
+ ()
+ end
+
+fun norm_Constr (constr as cplexConstr (c, (t1, t2))) =
+ if not (modulo_signed is_Num t2) andalso
+ modulo_signed is_Num t1
+ then
+ [cplexConstr (rev_cmp c, (t2, t1))]
+ else if (c = cplexLe orelse c = cplexLeq) andalso
+ (t1 = (cplexNeg cplexInf) orelse t2 = cplexInf)
+ then
+ []
+ else if (c = cplexGe orelse c = cplexGeq) andalso
+ (t1 = cplexInf orelse t2 = cplexNeg cplexInf)
+ then
+ []
+ else
+ [constr]
+
+fun bound2constr (cplexBounds (t1,c1,t2,c2,t3)) =
+ (norm_Constr(cplexConstr (c1, (t1, t2))))
+ @ (norm_Constr(cplexConstr (c2, (t2, t3))))
+ | bound2constr (cplexBound (t1, cplexEq, t2)) =
+ (norm_Constr(cplexConstr (cplexLeq, (t1, t2))))
+ @ (norm_Constr(cplexConstr (cplexLeq, (t2, t1))))
+ | bound2constr (cplexBound (t1, c1, t2)) =
+ norm_Constr(cplexConstr (c1, (t1,t2)))
+
+val emptyset = Symtab.empty
+
+fun singleton v = Symtab.update ((v, ()), emptyset)
+
+fun merge a b = Symtab.merge (op =) (a, b)
+
+fun mergemap f ts = foldl
+ (fn (table, x) => merge table (f x))
+ (Symtab.empty, ts)
+
+fun diff a b = Symtab.foldl (fn (a, (k,v)) =>
+ (Symtab.delete k a) handle UNDEF => a)
+ (a, b)
+
+fun collect_vars (cplexVar v) = singleton v
+ | collect_vars (cplexNeg t) = collect_vars t
+ | collect_vars (cplexProd (t1, t2)) =
+ merge (collect_vars t1) (collect_vars t2)
+ | collect_vars (cplexSum ts) = mergemap collect_vars ts
+ | collect_vars _ = emptyset
+
+(* Eliminates all nonfree bounds from the linear program and produces an
+ equivalent program with only free bounds
+ IF for the input program P holds: is_normed_cplexProg P *)
+fun elim_nonfree_bounds (cplexProg (name, goal, constraints, bounds)) =
+ let
+ fun collect_constr_vars (_, cplexConstr (c, (t1,_))) =
+ (collect_vars t1)
+
+ val cvars = merge (collect_vars (term_of_goal goal))
+ (mergemap collect_constr_vars constraints)
+
+ fun collect_lower_bounded_vars
+ (cplexBounds (t1, c1, cplexVar v, c2, t3)) =
+ singleton v
+ | collect_lower_bounded_vars
+ (cplexBound (_, cplexLe, cplexVar v)) =
+ singleton v
+ | collect_lower_bounded_vars
+ (cplexBound (_, cplexLeq, cplexVar v)) =
+ singleton v
+ | collect_lower_bounded_vars
+ (cplexBound (cplexVar v, cplexGe,_)) =
+ singleton v
+ | collect_lower_bounded_vars
+ (cplexBound (cplexVar v, cplexGeq, _)) =
+ singleton v
+ | collect_lower_bounded_vars
+ (cplexBound (cplexVar v, cplexEq, _)) =
+ singleton v
+ | collect_lower_bounded_vars _ = emptyset
+
+ val lvars = mergemap collect_lower_bounded_vars bounds
+ val positive_vars = diff cvars lvars
+ val zero = cplexNum "0"
+
+ fun make_pos_constr v =
+ (None, cplexConstr (cplexGeq, ((cplexVar v), zero)))
+
+ fun make_free_bound v =
+ cplexBounds (cplexNeg cplexInf, cplexLeq,
+ cplexVar v, cplexLeq,
+ cplexInf)
+
+ val pos_constrs = rev(Symtab.foldl
+ (fn (l, (k,v)) => (make_pos_constr k) :: l)
+ ([], positive_vars))
+ val bound_constrs = map (fn x => (None, x))
+ (flat (map bound2constr bounds))
+ val constraints' = constraints @ pos_constrs @ bound_constrs
+ val bounds' = rev(Symtab.foldl (fn (l, (v,_)) =>
+ (make_free_bound v)::l)
+ ([], cvars))
+ in
+ cplexProg (name, goal, constraints', bounds')
+ end
+
+fun relax_strict_ineqs (cplexProg (name, goals, constrs, bounds)) =
+ let
+ fun relax cplexLe = cplexLeq
+ | relax cplexGe = cplexGeq
+ | relax x = x
+
+ fun relax_constr (n, cplexConstr(c, (t1, t2))) =
+ (n, cplexConstr(relax c, (t1, t2)))
+
+ fun relax_bounds (cplexBounds (t1, c1, t2, c2, t3)) =
+ cplexBounds (t1, relax c1, t2, relax c2, t3)
+ | relax_bounds (cplexBound (t1, c, t2)) =
+ cplexBound (t1, relax c, t2)
+ in
+ cplexProg (name,
+ goals,
+ map relax_constr constrs,
+ map relax_bounds bounds)
+ end
+
+datatype cplexResult = Unbounded
+ | Infeasible
+ | Undefined
+ | Optimal of string * ((string * string) list)
+
+fun is_separator x = forall (fn c => c = #"-") (String.explode x)
+
+fun is_sign x = (x = "+" orelse x = "-")
+
+fun is_colon x = (x = ":")
+
+fun is_resultsymbol a =
+ let
+ val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~-"
+ fun is_symbol_char c = Char.isAlphaNum c orelse
+ exists (fn d => d=c) symbol_char
+ fun is_symbol_start c = is_symbol_char c andalso
+ not (Char.isDigit c) andalso
+ not (c= #".") andalso
+ not (c= #"-")
+ val b = String.explode a
+ in
+ b <> [] andalso is_symbol_start (hd b) andalso
+ forall is_symbol_char b
+ end
+
+val TOKEN_SIGN = 100
+val TOKEN_COLON = 101
+val TOKEN_SEPARATOR = 102
+
+fun load_cplexResult name =
+ let
+ val flist = [(is_NL, TOKEN_NL),
+ (is_blank, TOKEN_BLANK),
+ (is_num, TOKEN_NUM),
+ (is_sign, TOKEN_SIGN),
+ (is_colon, TOKEN_COLON),
+ (is_cmp, TOKEN_CMP),
+ (is_resultsymbol, TOKEN_SYMBOL),
+ (is_separator, TOKEN_SEPARATOR)]
+
+ val tokenize = tokenize_general flist
+
+ val f = TextIO.openIn name
+
+ val rest = ref []
+
+ fun readToken_helper () =
+ if length (!rest) > 0 then
+ let val u = hd (!rest) in
+ (
+ rest := tl (!rest);
+ Some u
+ )
+ end
+ else
+ let val s = TextIO.inputLine f in
+ if s = "" then None else (rest := tokenize s; readToken_helper())
+ end
+
+ fun is_tt tok ty = (tok <> None andalso (fst (the tok)) = ty)
+
+ fun pushToken a = if a = None then () else (rest := ((the a)::(!rest)))
+
+ fun readToken () =
+ let val t = readToken_helper () in
+ if is_tt t TOKEN_BLANK then
+ readToken ()
+ else if is_tt t TOKEN_NL then
+ let val t2 = readToken_helper () in
+ if is_tt t2 TOKEN_SIGN then
+ (pushToken (Some (TOKEN_SEPARATOR, "-")); t)
+ else
+ (pushToken t2; t)
+ end
+ else if is_tt t TOKEN_SIGN then
+ let val t2 = readToken_helper () in
+ if is_tt t2 TOKEN_NUM then
+ (Some (TOKEN_NUM, (snd (the t))^(snd (the t2))))
+ else
+ (pushToken t2; t)
+ end
+ else
+ t
+ end
+
+ fun readRestOfLine P =
+ let
+ val t = readToken ()
+ in
+ if is_tt t TOKEN_NL orelse t = None
+ then P
+ else readRestOfLine P
+ end
+
+ fun readHeader () =
+ let
+ fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ())))
+ fun readObjective () = readRestOfLine ("OBJECTIVE", snd (the (readToken (); readToken (); readToken ())))
+ val t1 = readToken ()
+ val t2 = readToken ()
+ in
+ if is_tt t1 TOKEN_SYMBOL andalso is_tt t2 TOKEN_COLON
+ then
+ case to_upper (snd (the t1)) of
+ "STATUS" => (readStatus ())::(readHeader ())
+ | "OBJECTIVE" => (readObjective())::(readHeader ())
+ | _ => (readRestOfLine (); readHeader ())
+ else
+ (pushToken t2; pushToken t1; [])
+ end
+
+ fun skip_until_sep () =
+ let val x = readToken () in
+ if is_tt x TOKEN_SEPARATOR then
+ readRestOfLine ()
+ else
+ skip_until_sep ()
+ end
+
+ fun load_value () =
+ let
+ val t1 = readToken ()
+ val t2 = readToken ()
+ in
+ if is_tt t1 TOKEN_NUM andalso is_tt t2 TOKEN_SYMBOL then
+ let
+ val t = readToken ()
+ val state = if is_tt t TOKEN_NL then readToken () else t
+ val _ = if is_tt state TOKEN_SYMBOL then () else raise (Load_cplexResult "state expected")
+ val k = readToken ()
+ in
+ if is_tt k TOKEN_NUM then
+ readRestOfLine (Some (snd (the t2), snd (the k)))
+ else
+ raise (Load_cplexResult "number expected")
+ end
+ else
+ (pushToken t2; pushToken t1; None)
+ end
+
+ fun load_values () =
+ let val v = load_value () in
+ if v = None then [] else (the v)::(load_values ())
+ end
+
+ val header = readHeader ()
+
+ val result =
+ case assoc (header, "STATUS") of
+ Some "INFEASIBLE" => Infeasible
+ | Some "UNBOUNDED" => Unbounded
+ | Some "OPTIMAL" => Optimal (the (assoc (header, "OBJECTIVE")),
+ (skip_until_sep ();
+ skip_until_sep ();
+ load_values ()))
+ | _ => Undefined
+
+ val _ = TextIO.closeIn f
+ in
+ result
+ end
+ handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s))
+ | OPTION => raise (Load_cplexResult "OPTION")
+ | x => raise x
+
+exception Execute of string;
+
+fun execute_cplex lpname resultname =
+let
+ fun wrap s = "\""^s^"\""
+ val cplex_path = getenv "CPLEX"
+ val cplex = if cplex_path = "" then "glpsol" else cplex_path
+ val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname)
+in
+ execute command
+end
+
+fun tmp_file s = Path.pack (Path.expand (File.tmp_path (Path.make [s])))
+
+fun solve prog =
+ let
+ val name = LargeInt.toString (Time.toMicroseconds (Time.now ()))
+ val lpname = tmp_file (name^".lp")
+ val resultname = tmp_file (name^".result")
+ val _ = save_cplexFile lpname prog
+ val answer = execute_cplex lpname resultname
+ in
+ let
+ val result = load_cplexResult resultname
+ val _ = OS.FileSys.remove lpname
+ val _ = OS.FileSys.remove resultname
+ in
+ result
+ end
+ handle (Load_cplexResult s) => raise (Execute ("Load_cplexResult: "^s^"\nExecute: "^answer))
+ | _ => raise (Execute answer)
+ end
+end;
+
+val demofile = "/home/obua/flyspeck/kepler/LP/cplexPent2.lp45"
+val demoout = "/home/obua/flyspeck/kepler/LP/test.out"
+val demoresult = "/home/obua/flyspeck/kepler/LP/try/test2.sol"
+
+fun loadcplex () = Cplex.relax_strict_ineqs
+ (Cplex.load_cplexFile demofile)
+
+fun writecplex lp = Cplex.save_cplexFile demoout lp
+
+fun test () =
+ let
+ val lp = loadcplex ()
+ val lp2 = Cplex.elim_nonfree_bounds lp
+ in
+ writecplex lp2
+ end
+
+fun loadresult () = Cplex.load_cplexResult demoresult;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/CplexMatrixConverter.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,132 @@
+signature MATRIX_BUILDER =
+sig
+ type vector
+ type matrix
+
+ val empty_vector : vector
+ val empty_matrix : matrix
+
+ exception Nat_expected of int
+ val set_elem : vector -> int -> string -> vector
+ val set_vector : matrix -> int -> vector -> matrix
+end;
+
+signature CPLEX_MATRIX_CONVERTER =
+sig
+ structure cplex : CPLEX
+ structure matrix_builder : MATRIX_BUILDER
+ type vector = matrix_builder.vector
+ type matrix = matrix_builder.matrix
+ type naming = int * (int -> string) * (string -> int)
+
+ exception Converter of string
+
+ (* program must fulfill is_normed_cplexProg and must be an element of the image of elim_nonfree_bounds *)
+ (* convert_prog maximize c A b naming *)
+ val convert_prog : cplex.cplexProg -> bool * vector * matrix * vector * naming
+
+ (* results must be optimal, converts_results returns the optimal value as string and the solution as vector *)
+ (* convert_results results name2index *)
+ val convert_results : cplex.cplexResult -> (string -> int) -> string * vector
+end;
+
+functor MAKE_CPLEX_MATRIX_CONVERTER (structure cplex: CPLEX and matrix_builder: MATRIX_BUILDER) : CPLEX_MATRIX_CONVERTER =
+struct
+
+structure cplex = cplex
+structure matrix_builder = matrix_builder
+type matrix = matrix_builder.matrix
+type vector = matrix_builder.vector
+type naming = int * (int -> string) * (string -> int)
+
+open matrix_builder
+open cplex
+
+exception Converter of string;
+
+structure Inttab = TableFun(type key = int val ord = int_ord);
+
+fun neg_term (cplexNeg t) = t
+ | neg_term (cplexSum ts) = cplexSum (map neg_term ts)
+ | neg_term t = cplexNeg t
+
+fun convert_prog (cplexProg (s, goal, constrs, bounds)) =
+ let
+ fun build_naming index i2s s2i [] = (index, i2s, s2i)
+ | build_naming index i2s s2i (cplexBounds (cplexNeg cplexInf, cplexLeq, cplexVar v, cplexLeq, cplexInf)::bounds)
+ = build_naming (index+1) (Inttab.update ((index, v), i2s)) (Symtab.update_new ((v, index), s2i)) bounds
+ | build_naming _ _ _ _ = raise (Converter "nonfree bound")
+
+ val (varcount, i2s_tab, s2i_tab) = build_naming 0 Inttab.empty Symtab.empty bounds
+
+ fun i2s i = case Inttab.lookup (i2s_tab, i) of None => raise (Converter "index not found")
+ | Some n => n
+ fun s2i s = case Symtab.lookup (s2i_tab, s) of None => raise (Converter ("name not found: "^s))
+ | Some i => i
+ fun num2str positive (cplexNeg t) = num2str (not positive) t
+ | num2str positive (cplexNum num) = if positive then num else "-"^num
+ | num2str _ _ = raise (Converter "term is not a (possibly signed) number")
+
+ fun setprod vec positive (cplexNeg t) = setprod vec (not positive) t
+ | setprod vec positive (cplexVar v) = set_elem vec (s2i v) (if positive then "1" else "-1")
+ | setprod vec positive (cplexProd (cplexNum num, cplexVar v)) =
+ set_elem vec (s2i v) (if positive then num else "-"^num)
+ | setprod _ _ _ = raise (Converter "term is not a normed product")
+
+ fun sum2vec (cplexSum ts) = foldl (fn (vec, t) => setprod vec true t) (empty_vector, ts)
+ | sum2vec t = setprod empty_vector true t
+
+ fun constrs2Ab j A b [] = (A, b)
+ | constrs2Ab j A b ((_, cplexConstr (cplexLeq, (t1,t2)))::cs) =
+ constrs2Ab (j+1) (set_vector A j (sum2vec t1)) (set_elem b j (num2str true t2)) cs
+ | constrs2Ab j A b ((_, cplexConstr (cplexGeq, (t1,t2)))::cs) =
+ constrs2Ab (j+1) (set_vector A j (sum2vec (neg_term t1))) (set_elem b j (num2str true (neg_term t2))) cs
+ | constrs2Ab j A b ((_, cplexConstr (cplexEq, (t1,t2)))::cs) =
+ constrs2Ab j A b ((None, cplexConstr (cplexLeq, (t1,t2)))::
+ (None, cplexConstr (cplexGeq, (t1, t2)))::cs)
+ | constrs2Ab _ _ _ _ = raise (Converter "no strict constraints allowed")
+
+ val (A, b) = constrs2Ab 0 empty_matrix empty_vector constrs
+
+ val (goal_maximize, goal_term) =
+ case goal of
+ (cplexMaximize t) => (true, t)
+ | (cplexMinimize t) => (false, t)
+ in
+ (goal_maximize, sum2vec goal_term, A, b, (varcount, i2s, s2i))
+ end
+
+fun convert_results (cplex.Optimal (opt, entries)) name2index =
+ let
+ fun setv (v, (name, value)) = (matrix_builder.set_elem v (name2index name) value)
+ in
+ (opt, foldl setv (matrix_builder.empty_vector, entries))
+ end
+ | convert_results _ _ = raise (Converter "No optimal result")
+
+
+end;
+
+structure SimpleMatrixBuilder : MATRIX_BUILDER =
+struct
+type vector = (int * string) list
+type matrix = (int * vector) list
+
+val empty_matrix = []
+val empty_vector = []
+
+exception Nat_expected of int;
+
+fun set_elem v i s = v @ [(i, s)]
+
+fun set_vector m i v = m @ [(i, v)]
+
+end;
+
+
+
+structure SimpleCplexMatrixConverter = MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = SimpleMatrixBuilder);
+
+
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/ExactFloatingPoint.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,215 @@
+structure ExactFloatingPoint :
+sig
+ exception Destruct_floatstr of string
+
+ val destruct_floatstr : (char -> bool) -> (char -> bool) -> string -> bool * string * string * bool * string
+
+ exception Floating_point of string
+
+ type floatrep = IntInf.int * IntInf.int
+ val approx_dec_by_bin : IntInf.int -> floatrep -> floatrep * floatrep
+ val approx_decstr_by_bin : int -> string -> floatrep * floatrep
+end
+=
+struct
+
+fun fst (a,b) = a
+fun snd (a,b) = b
+
+val filter = List.filter;
+
+exception Destruct_floatstr of string;
+
+fun destruct_floatstr isDigit isExp number =
+ let
+ val numlist = filter (not o Char.isSpace) (String.explode number)
+
+ fun countsigns ((#"+")::cs) = countsigns cs
+ | countsigns ((#"-")::cs) =
+ let
+ val (positive, rest) = countsigns cs
+ in
+ (not positive, rest)
+ end
+ | countsigns cs = (true, cs)
+
+ fun readdigits [] = ([], [])
+ | readdigits (q as c::cs) =
+ if (isDigit c) then
+ let
+ val (digits, rest) = readdigits cs
+ in
+ (c::digits, rest)
+ end
+ else
+ ([], q)
+
+ fun readfromexp_helper cs =
+ let
+ val (positive, rest) = countsigns cs
+ val (digits, rest') = readdigits rest
+ in
+ case rest' of
+ [] => (positive, digits)
+ | _ => raise (Destruct_floatstr number)
+ end
+
+ fun readfromexp [] = (true, [])
+ | readfromexp (c::cs) =
+ if isExp c then
+ readfromexp_helper cs
+ else
+ raise (Destruct_floatstr number)
+
+ fun readfromdot [] = ([], readfromexp [])
+ | readfromdot ((#".")::cs) =
+ let
+ val (digits, rest) = readdigits cs
+ val exp = readfromexp rest
+ in
+ (digits, exp)
+ end
+ | readfromdot cs = readfromdot ((#".")::cs)
+
+ val (positive, numlist) = countsigns numlist
+ val (digits1, numlist) = readdigits numlist
+ val (digits2, exp) = readfromdot numlist
+ in
+ (positive, String.implode digits1, String.implode digits2, fst exp, String.implode (snd exp))
+ end
+
+type floatrep = IntInf.int * IntInf.int
+
+exception Floating_point of string;
+
+val ln2_10 = (Math.ln 10.0)/(Math.ln 2.0)
+
+fun intmul a b = IntInf.* (a,b)
+fun intsub a b = IntInf.- (a,b)
+fun intadd a b = IntInf.+ (a,b)
+fun intpow a b = IntInf.pow (a, IntInf.toInt b);
+fun intle a b = IntInf.<= (a, b);
+fun intless a b = IntInf.< (a, b);
+fun intneg a = IntInf.~ a;
+val zero = IntInf.fromInt 0;
+val one = IntInf.fromInt 1;
+val two = IntInf.fromInt 2;
+val ten = IntInf.fromInt 10;
+val five = IntInf.fromInt 5;
+
+fun find_most_significant q r =
+ let
+ fun int2real i =
+ case Real.fromString (IntInf.toString i) of
+ SOME r => r
+ | NONE => raise (Floating_point "int2real")
+ fun subtract (q, r) (q', r') =
+ if intle r r' then
+ (intsub q (intmul q' (intpow ten (intsub r' r))), r)
+ else
+ (intsub (intmul q (intpow ten (intsub r r'))) q', r')
+ fun bin2dec d =
+ if intle zero d then
+ (intpow two d, zero)
+ else
+ (intpow five (intneg d), d)
+
+ val L = IntInf.fromInt (Real.floor (int2real (IntInf.fromInt (IntInf.log2 q)) + (int2real r) * ln2_10))
+ val L1 = intadd L one
+
+ val (q1, r1) = subtract (q, r) (bin2dec L1)
+ in
+ if intle zero q1 then
+ let
+ val (q2, r2) = subtract (q, r) (bin2dec (intadd L1 one))
+ in
+ if intle zero q2 then
+ raise (Floating_point "find_most_significant")
+ else
+ (L1, (q1, r1))
+ end
+ else
+ let
+ val (q0, r0) = subtract (q, r) (bin2dec L)
+ in
+ if intle zero q0 then
+ (L, (q0, r0))
+ else
+ raise (Floating_point "find_most_significant")
+ end
+ end
+
+fun approx_dec_by_bin n (q,r) =
+ let
+ fun addseq acc d' [] = acc
+ | addseq acc d' (d::ds) = addseq (intadd acc (intpow two (intsub d d'))) d' ds
+
+ fun seq2bin [] = (zero, zero)
+ | seq2bin (d::ds) = (intadd (addseq zero d ds) one, d)
+
+ fun approx d_seq d0 precision (q,r) =
+ if q = zero then
+ let val x = seq2bin d_seq in
+ (x, x)
+ end
+ else
+ let
+ val (d, (q', r')) = find_most_significant q r
+ in
+ if intless precision (intsub d0 d) then
+ let
+ val d' = intsub d0 precision
+ val x1 = seq2bin (d_seq)
+ val x2 = (intadd (intmul (fst x1) (intpow two (intsub (snd x1) d'))) one, d') (* = seq2bin (d'::d_seq) *)
+ in
+ (x1, x2)
+ end
+ else
+ approx (d::d_seq) d0 precision (q', r')
+ end
+
+ fun approx_start precision (q, r) =
+ if q = zero then
+ ((zero, zero), (zero, zero))
+ else
+ let
+ val (d, (q', r')) = find_most_significant q r
+ in
+ if intle precision zero then
+ let
+ val x1 = seq2bin [d]
+ in
+ if q' = zero then
+ (x1, x1)
+ else
+ (x1, seq2bin [intadd d one])
+ end
+ else
+ approx [d] d precision (q', r')
+ end
+ in
+ if intle zero q then
+ approx_start n (q,r)
+ else
+ let
+ val ((a1,b1), (a2, b2)) = approx_start n (intneg q, r)
+ in
+ ((intneg a2, b2), (intneg a1, b1))
+ end
+ end
+
+fun approx_decstr_by_bin n decstr =
+ let
+ fun str2int s = case IntInf.fromString s of SOME x => x | NONE => zero
+ fun signint p x = if p then x else intneg x
+
+ val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr
+ val s = IntInf.fromInt (size d2)
+
+ val q = signint p (intadd (intmul (str2int d1) (intpow ten s)) (str2int d2))
+ val r = intsub (signint ep (str2int e)) s
+ in
+ approx_dec_by_bin (IntInf.fromInt n) (q,r)
+ end
+
+end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/Float.thy Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,508 @@
+theory Float = Hyperreal:
+
+constdefs
+ pow2 :: "int \<Rightarrow> real"
+ "pow2 a == if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a))))"
+ float :: "int * int \<Rightarrow> real"
+ "float x == (real (fst x)) * (pow2 (snd x))"
+
+lemma pow2_0[simp]: "pow2 0 = 1"
+by (simp add: pow2_def)
+
+lemma pow2_1[simp]: "pow2 1 = 2"
+by (simp add: pow2_def)
+
+lemma pow2_neg: "pow2 x = inverse (pow2 (-x))"
+by (simp add: pow2_def)
+
+lemma pow2_add1: "pow2 (1 + a) = 2 * (pow2 a)"
+proof -
+ have h: "! n. nat (2 + int n) - Suc 0 = nat (1 + int n)" by arith
+ have g: "! a b. a - -1 = a + (1::int)" by arith
+ have pos: "! n. pow2 (int n + 1) = 2 * pow2 (int n)"
+ apply (auto, induct_tac n)
+ apply (simp_all add: pow2_def)
+ apply (rule_tac m1="2" and n1="nat (2 + int na)" in ssubst[OF realpow_num_eq_if])
+ apply (auto simp add: h)
+ apply arith
+ done
+ show ?thesis
+ proof (induct a)
+ case (1 n)
+ from pos show ?case by (simp add: ring_eq_simps)
+ next
+ case (2 n)
+ show ?case
+ apply (auto)
+ apply (subst pow2_neg[of "- int n"])
+ apply (subst pow2_neg[of "-1 - int n"])
+ apply (auto simp add: g pos)
+ done
+ qed
+qed
+
+lemma pow2_add: "pow2 (a+b) = (pow2 a) * (pow2 b)"
+proof (induct b)
+ case (1 n)
+ show ?case
+ proof (induct n)
+ case 0
+ show ?case by simp
+ next
+ case (Suc m)
+ show ?case by (auto simp add: ring_eq_simps pow2_add1 prems)
+ qed
+next
+ case (2 n)
+ show ?case
+ proof (induct n)
+ case 0
+ show ?case
+ apply (auto)
+ apply (subst pow2_neg[of "a + -1"])
+ apply (subst pow2_neg[of "-1"])
+ apply (simp)
+ apply (insert pow2_add1[of "-a"])
+ apply (simp add: ring_eq_simps)
+ apply (subst pow2_neg[of "-a"])
+ apply (simp)
+ done
+ case (Suc m)
+ have a: "int m - (a + -2) = 1 + (int m - a + 1)" by arith
+ have b: "int m - -2 = 1 + (int m + 1)" by arith
+ show ?case
+ apply (auto)
+ apply (subst pow2_neg[of "a + (-2 - int m)"])
+ apply (subst pow2_neg[of "-2 - int m"])
+ apply (auto simp add: ring_eq_simps)
+ apply (subst a)
+ apply (subst b)
+ apply (simp only: pow2_add1)
+ apply (subst pow2_neg[of "int m - a + 1"])
+ apply (subst pow2_neg[of "int m + 1"])
+ apply auto
+ apply (insert prems)
+ apply (auto simp add: ring_eq_simps)
+ done
+ qed
+qed
+
+lemma "float (a, e) + float (b, e) = float (a + b, e)"
+by (simp add: float_def ring_eq_simps)
+
+constdefs
+ int_of_real :: "real \<Rightarrow> int"
+ "int_of_real x == SOME y. real y = x"
+ real_is_int :: "real \<Rightarrow> bool"
+ "real_is_int x == ? (u::int). x = real u"
+
+lemma real_is_int_def2: "real_is_int x = (x = real (int_of_real x))"
+by (auto simp add: real_is_int_def int_of_real_def)
+
+lemma float_transfer: "real_is_int ((real a)*(pow2 c)) \<Longrightarrow> float (a, b) = float (int_of_real ((real a)*(pow2 c)), b - c)"
+by (simp add: float_def real_is_int_def2 pow2_add[symmetric])
+
+lemma pow2_int: "pow2 (int c) = (2::real)^c"
+by (simp add: pow2_def)
+
+lemma float_transfer_nat: "float (a, b) = float (a * 2^c, b - int c)"
+by (simp add: float_def pow2_int[symmetric] pow2_add[symmetric])
+
+lemma real_is_int_real[simp]: "real_is_int (real (x::int))"
+by (auto simp add: real_is_int_def int_of_real_def)
+
+lemma int_of_real_real[simp]: "int_of_real (real x) = x"
+by (simp add: int_of_real_def)
+
+lemma real_int_of_real[simp]: "real_is_int x \<Longrightarrow> real (int_of_real x) = x"
+by (auto simp add: int_of_real_def real_is_int_def)
+
+lemma real_is_int_add_int_of_real: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> (int_of_real (a+b)) = (int_of_real a) + (int_of_real b)"
+by (auto simp add: int_of_real_def real_is_int_def)
+
+lemma real_is_int_add[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a+b)"
+apply (subst real_is_int_def2)
+apply (simp add: real_is_int_add_int_of_real real_int_of_real)
+done
+
+lemma int_of_real_sub: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> (int_of_real (a-b)) = (int_of_real a) - (int_of_real b)"
+by (auto simp add: int_of_real_def real_is_int_def)
+
+lemma real_is_int_sub[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a-b)"
+apply (subst real_is_int_def2)
+apply (simp add: int_of_real_sub real_int_of_real)
+done
+
+lemma real_is_int_rep: "real_is_int x \<Longrightarrow> ?! (a::int). real a = x"
+by (auto simp add: real_is_int_def)
+
+lemma int_of_real_mult:
+ assumes "real_is_int a" "real_is_int b"
+ shows "(int_of_real (a*b)) = (int_of_real a) * (int_of_real b)"
+proof -
+ from prems have a: "?! (a'::int). real a' = a" by (rule_tac real_is_int_rep, auto)
+ from prems have b: "?! (b'::int). real b' = b" by (rule_tac real_is_int_rep, auto)
+ from a obtain a'::int where a':"a = real a'" by auto
+ from b obtain b'::int where b':"b = real b'" by auto
+ have r: "real a' * real b' = real (a' * b')" by auto
+ show ?thesis
+ apply (simp add: a' b')
+ apply (subst r)
+ apply (simp only: int_of_real_real)
+ done
+qed
+
+lemma real_is_int_mult[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a*b)"
+apply (subst real_is_int_def2)
+apply (simp add: int_of_real_mult)
+done
+
+lemma real_is_int_0[simp]: "real_is_int (0::real)"
+by (simp add: real_is_int_def int_of_real_def)
+
+lemma real_is_int_1[simp]: "real_is_int (1::real)"
+proof -
+ have "real_is_int (1::real) = real_is_int(real (1::int))" by auto
+ also have "\<dots> = True" by (simp only: real_is_int_real)
+ ultimately show ?thesis by auto
+qed
+
+lemma real_is_int_n1: "real_is_int (-1::real)"
+proof -
+ have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto
+ also have "\<dots> = True" by (simp only: real_is_int_real)
+ ultimately show ?thesis by auto
+qed
+
+lemma real_is_int_number_of[simp]: "real_is_int ((number_of::bin\<Rightarrow>real) x)"
+proof -
+ have neg1: "real_is_int (-1::real)"
+ proof -
+ have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto
+ also have "\<dots> = True" by (simp only: real_is_int_real)
+ ultimately show ?thesis by auto
+ qed
+
+ {
+ fix x::int
+ have "!! y. real_is_int ((number_of::bin\<Rightarrow>real) (Abs_Bin x))"
+ apply (simp add: number_of_eq)
+ apply (subst Abs_Bin_inverse)
+ apply (simp add: Bin_def)
+ apply (induct x)
+ apply (induct_tac n)
+ apply (simp)
+ apply (simp)
+ apply (induct_tac n)
+ apply (simp add: neg1)
+ proof -
+ fix n :: nat
+ assume rn: "(real_is_int (of_int (- (int (Suc n)))))"
+ have s: "-(int (Suc (Suc n))) = -1 + - (int (Suc n))" by simp
+ show "real_is_int (of_int (- (int (Suc (Suc n)))))"
+ apply (simp only: s of_int_add)
+ apply (rule real_is_int_add)
+ apply (simp add: neg1)
+ apply (simp only: rn)
+ done
+ qed
+ }
+ note Abs_Bin = this
+ {
+ fix x :: bin
+ have "? u. x = Abs_Bin u"
+ apply (rule exI[where x = "Rep_Bin x"])
+ apply (simp add: Rep_Bin_inverse)
+ done
+ }
+ then obtain u::int where "x = Abs_Bin u" by auto
+ with Abs_Bin show ?thesis by auto
+qed
+
+lemma int_of_real_0[simp]: "int_of_real (0::real) = (0::int)"
+by (simp add: int_of_real_def)
+
+lemma int_of_real_1[simp]: "int_of_real (1::real) = (1::int)"
+proof -
+ have 1: "(1::real) = real (1::int)" by auto
+ show ?thesis by (simp only: 1 int_of_real_real)
+qed
+
+lemma int_of_real_number_of[simp]: "int_of_real (number_of b) = number_of b"
+proof -
+ have "real_is_int (number_of b)" by simp
+ then have uu: "?! u::int. number_of b = real u" by (auto simp add: real_is_int_rep)
+ then obtain u::int where u:"number_of b = real u" by auto
+ have "number_of b = real ((number_of b)::int)"
+ by (simp add: number_of_eq real_of_int_def)
+ have ub: "number_of b = real ((number_of b)::int)"
+ by (simp add: number_of_eq real_of_int_def)
+ from uu u ub have unb: "u = number_of b"
+ by blast
+ have "int_of_real (number_of b) = u" by (simp add: u)
+ with unb show ?thesis by simp
+qed
+
+lemma float_transfer_even: "even a \<Longrightarrow> float (a, b) = float (a div 2, b+1)"
+ apply (subst float_transfer[where a="a" and b="b" and c="-1", simplified])
+ apply (simp_all add: pow2_def even_def real_is_int_def ring_eq_simps)
+ apply (auto)
+proof -
+ fix q::int
+ have a:"b - (-1\<Colon>int) = (1\<Colon>int) + b" by arith
+ show "(float (q, (b - (-1\<Colon>int)))) = (float (q, ((1\<Colon>int) + b)))"
+ by (simp add: a)
+qed
+
+consts
+ norm_float :: "int*int \<Rightarrow> int*int"
+
+lemma int_div_zdiv: "int (a div b) = (int a) div (int b)"
+apply (subst split_div, auto)
+apply (subst split_zdiv, auto)
+apply (rule_tac a="int (b * i) + int j" and b="int b" and r="int j" and r'=ja in IntDiv.unique_quotient)
+apply (auto simp add: IntDiv.quorem_def int_eq_of_nat)
+done
+
+lemma int_mod_zmod: "int (a mod b) = (int a) mod (int b)"
+apply (subst split_mod, auto)
+apply (subst split_zmod, auto)
+apply (rule_tac a="int (b * i) + int j" and b="int b" and q="int i" and q'=ia in IntDiv.unique_remainder)
+apply (auto simp add: IntDiv.quorem_def int_eq_of_nat)
+done
+
+lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> abs((a::int) div 2) < abs a"
+by arith
+
+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>"
+apply (auto)
+apply (rule abs_div_2_less)
+apply (auto)
+done
+
+ML {* simp_depth_limit := 2 *}
+recdef norm_float "measure (% (a,b). nat (abs a))"
+ "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)))"
+(hints simp: terminating_norm_float)
+ML {* simp_depth_limit := 1000 *}
+
+
+lemma norm_float: "float x = float (norm_float x)"
+proof -
+ {
+ fix a b :: int
+ have norm_float_pair: "float (a,b) = float (norm_float (a,b))"
+ proof (induct a b rule: norm_float.induct)
+ case (1 u v)
+ show ?case
+ proof cases
+ assume u: "u \<noteq> 0 \<and> even u"
+ with prems have ind: "float (u div 2, v + 1) = float (norm_float (u div 2, v + 1))" by auto
+ with u have "float (u,v) = float (u div 2, v+1)" by (simp add: float_transfer_even)
+ then show ?thesis
+ apply (subst norm_float.simps)
+ apply (simp add: ind)
+ done
+ next
+ assume "~(u \<noteq> 0 \<and> even u)"
+ then show ?thesis
+ by (simp add: prems float_def)
+ qed
+ qed
+ }
+ note helper = this
+ have "? a b. x = (a,b)" by auto
+ then obtain a b where "x = (a, b)" by blast
+ then show ?thesis by (simp only: helper)
+qed
+
+lemma pow2_int: "pow2 (int n) = 2^n"
+ by (simp add: pow2_def)
+
+lemma float_add:
+ "float (a1, e1) + float (a2, e2) =
+ (if e1<=e2 then float (a1+a2*2^(nat(e2-e1)), e1)
+ else float (a1*2^(nat (e1-e2))+a2, e2))"
+ apply (simp add: float_def ring_eq_simps)
+ apply (auto simp add: pow2_int[symmetric] pow2_add[symmetric])
+ done
+
+lemma float_mult:
+ "float (a1, e1) * float (a2, e2) =
+ (float (a1 * a2, e1 + e2))"
+ by (simp add: float_def pow2_add)
+
+lemma float_minus:
+ "- (float (a,b)) = float (-a, b)"
+ by (simp add: float_def)
+
+lemma zero_less_pow2:
+ "0 < pow2 x"
+proof -
+ {
+ fix y
+ have "0 <= y \<Longrightarrow> 0 < pow2 y"
+ by (induct y, induct_tac n, simp_all add: pow2_add)
+ }
+ note helper=this
+ show ?thesis
+ apply (case_tac "0 <= x")
+ apply (simp add: helper)
+ apply (subst pow2_neg)
+ apply (simp add: helper)
+ done
+qed
+
+lemma zero_le_float:
+ "(0 <= float (a,b)) = (0 <= a)"
+ apply (auto simp add: float_def)
+ apply (auto simp add: zero_le_mult_iff zero_less_pow2)
+ apply (insert zero_less_pow2[of b])
+ apply (simp_all)
+ done
+
+lemma float_abs:
+ "abs (float (a,b)) = (if 0 <= a then (float (a,b)) else (float (-a,b)))"
+ apply (auto simp add: abs_if)
+ apply (simp_all add: zero_le_float[symmetric, of a b] float_minus)
+ done
+
+lemma norm_0_1: "(0::_::number_ring) = Numeral0 & (1::_::number_ring) = Numeral1"
+ by auto
+
+lemma add_left_zero: "0 + a = (a::'a::comm_monoid_add)"
+ by simp
+
+lemma add_right_zero: "a + 0 = (a::'a::comm_monoid_add)"
+ by simp
+
+lemma mult_left_one: "1 * a = (a::'a::semiring_1)"
+ by simp
+
+lemma mult_right_one: "a * 1 = (a::'a::semiring_1)"
+ by simp
+
+lemma int_pow_0: "(a::int)^(Numeral0) = 1"
+ by simp
+
+lemma int_pow_1: "(a::int)^(Numeral1) = a"
+ by simp
+
+lemma zero_eq_Numeral0_nring: "(0::'a::number_ring) = Numeral0"
+ by simp
+
+lemma one_eq_Numeral1_nring: "(1::'a::number_ring) = Numeral1"
+ by simp
+
+lemma zero_eq_Numeral0_nat: "(0::nat) = Numeral0"
+ by simp
+
+lemma one_eq_Numeral1_nat: "(1::nat) = Numeral1"
+ by simp
+
+lemma zpower_Pls: "(z::int)^Numeral0 = Numeral1"
+ by simp
+
+lemma zpower_Min: "(z::int)^((-1)::nat) = Numeral1"
+proof -
+ have 1:"((-1)::nat) = 0"
+ by simp
+ show ?thesis by (simp add: 1)
+qed
+
+lemma fst_cong: "a=a' \<Longrightarrow> fst (a,b) = fst (a',b)"
+ by simp
+
+lemma snd_cong: "b=b' \<Longrightarrow> snd (a,b) = snd (a,b')"
+ by simp
+
+lemma lift_bool: "x \<Longrightarrow> x=True"
+ by simp
+
+lemma nlift_bool: "~x \<Longrightarrow> x=False"
+ by simp
+
+lemma not_false_eq_true: "(~ False) = True" by simp
+
+lemma not_true_eq_false: "(~ True) = False" by simp
+
+
+lemmas binarith =
+ Pls_0_eq Min_1_eq
+ bin_pred_Pls bin_pred_Min bin_pred_1 bin_pred_0
+ bin_succ_Pls bin_succ_Min bin_succ_1 bin_succ_0
+ bin_add_Pls bin_add_Min bin_add_BIT_0 bin_add_BIT_10
+ bin_add_BIT_11 bin_minus_Pls bin_minus_Min bin_minus_1
+ bin_minus_0 bin_mult_Pls bin_mult_Min bin_mult_1 bin_mult_0
+ bin_add_Pls_right bin_add_Min_right
+
+thm eq_number_of_eq
+
+lemma int_eq_number_of_eq: "(((number_of v)::int)=(number_of w)) = iszero ((number_of (bin_add v (bin_minus w)))::int)"
+ by simp
+
+lemma int_iszero_number_of_Pls: "iszero (Numeral0::int)"
+ by (simp only: iszero_number_of_Pls)
+
+lemma int_nonzero_number_of_Min: "~(iszero ((-1)::int))"
+ by simp
+thm iszero_number_of_1
+
+lemma int_iszero_number_of_0: "iszero ((number_of (w BIT False))::int) = iszero ((number_of w)::int)"
+ by simp
+
+lemma int_iszero_number_of_1: "\<not> iszero ((number_of (w BIT True))::int)"
+ by simp
+
+lemma int_less_number_of_eq_neg: "(((number_of x)::int) < number_of y) = neg ((number_of (bin_add x (bin_minus y)))::int)"
+ by simp
+
+lemma int_not_neg_number_of_Pls: "\<not> (neg (Numeral0::int))"
+ by simp
+
+lemma int_neg_number_of_Min: "neg (-1::int)"
+ by simp
+
+lemma int_neg_number_of_BIT: "neg ((number_of (w BIT x))::int) = neg ((number_of w)::int)"
+ by simp
+
+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))"
+ by simp
+
+lemmas intarithrel =
+ (*abs_zero abs_one*)
+ int_eq_number_of_eq
+ lift_bool[OF int_iszero_number_of_Pls] nlift_bool[OF int_nonzero_number_of_Min] int_iszero_number_of_0
+ 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]
+ int_neg_number_of_BIT int_le_number_of_eq
+
+lemma int_number_of_add_sym: "((number_of v)::int) + number_of w = number_of (bin_add v w)"
+ by simp
+
+lemma int_number_of_diff_sym: "((number_of v)::int) - number_of w = number_of (bin_add v (bin_minus w))"
+ by simp
+
+lemma int_number_of_mult_sym: "((number_of v)::int) * number_of w = number_of (bin_mult v w)"
+ by simp
+
+lemma int_number_of_minus_sym: "- ((number_of v)::int) = number_of (bin_minus v)"
+ by simp
+
+lemmas intarith = int_number_of_add_sym int_number_of_minus_sym int_number_of_diff_sym int_number_of_mult_sym
+
+lemmas natarith = (*zero_eq_Numeral0_nat one_eq_Numeral1_nat*) add_nat_number_of
+ diff_nat_number_of mult_nat_number_of eq_nat_number_of less_nat_number_of
+
+lemmas powerarith = nat_number_of zpower_number_of_even
+ zpower_number_of_odd[simplified zero_eq_Numeral0_nring one_eq_Numeral1_nring]
+ zpower_Pls zpower_Min
+
+lemmas floatarith[simplified norm_0_1] = float_add float_mult float_minus float_abs zero_le_float
+
+lemmas arith = binarith intarith intarithrel natarith powerarith floatarith not_false_eq_true not_true_eq_false
+
+(* needed for the verifying code generator *)
+lemmas arith_no_let = arith[simplified Let_def]
+
+end
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/FloatArith.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,56 @@
+structure FloatArith =
+struct
+
+type float = IntInf.int * IntInf.int
+
+val izero = IntInf.fromInt 0
+fun imul a b = IntInf.* (a,b)
+fun isub a b = IntInf.- (a,b)
+fun iadd a b = IntInf.+ (a,b)
+
+val floatzero = (izero, izero)
+
+fun positive_part (a,b) =
+ (if IntInf.< (a,izero) then izero else a, b)
+
+fun negative_part (a,b) =
+ (if IntInf.< (a,izero) then a else izero, b)
+
+fun is_negative (a,b) =
+ if IntInf.< (a, izero) then true else false
+
+fun is_positive (a,b) =
+ if IntInf.< (izero, a) then true else false
+
+fun is_zero (a,b) =
+ if a = izero then true else false
+
+fun ipow2 a = IntInf.pow ((IntInf.fromInt 2), IntInf.toInt a)
+
+fun add (a1, b1) (a2, b2) =
+ if b1 < b2 then
+ (iadd a1 (imul a2 (ipow2 (isub b2 b1))), b1)
+ else
+ (iadd (imul a1 (ipow2 (isub b1 b2))) a2, b2)
+
+fun sub (a1, b1) (a2, b2) =
+ if b1 < b2 then
+ (isub a1 (imul a2 (ipow2 (isub b2 b1))), b1)
+ else
+ (isub (imul a1 (ipow2 (isub b1 b2))) a2, b2)
+
+fun neg (a, b) = (IntInf.~ a, b)
+
+fun is_equal a b = is_zero (sub a b)
+
+fun is_less a b = is_negative (sub a b)
+
+fun max a b = if is_less a b then b else a
+
+fun min a b = if is_less a b then a else b
+
+fun abs a = if is_negative a then neg a else a
+
+fun mul (a1, b1) (a2, b2) = (imul a1 a2, iadd b1 b2)
+
+end;
--- a/src/HOL/Matrix/LinProg.thy Fri Sep 03 10:27:05 2004 +0200
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,79 +0,0 @@
-(* Title: HOL/Matrix/LinProg.thy
- ID: $Id$
- Author: Steven Obua
-*)
-
-theory LinProg = Matrix:
-
-lemma linprog_by_duality_approx_general:
- assumes
- fmuladdprops:
- "! a b c d. a <= b & c <= d \<longrightarrow> fadd a c <= fadd b d"
- "! c a b. 0 <= c & a <= b \<longrightarrow> fmul c a <= fmul c b"
- "! a. fmul 0 a = 0"
- "! a. fmul a 0 = 0"
- "fadd 0 0 = 0"
- "associative fadd"
- "commutative fadd"
- "associative fmul"
- "distributive fmul fadd"
- and specificprops:
- "mult_matrix fmul fadd (combine_matrix fadd A dA) x <= (b::('a::{order, zero}) matrix)"
- "mult_matrix fmul fadd y A = c"
- "0 <= y"
- shows
- "combine_matrix fadd (mult_matrix fmul fadd c x) (mult_matrix fmul fadd (mult_matrix fmul fadd y dA) x)
- <= mult_matrix fmul fadd y b"
-proof -
- let ?mul = "mult_matrix fmul fadd"
- let ?add = "combine_matrix fadd"
- let ?t1 = "?mul y (?mul (?add A dA) x)"
- have a: "?t1 <= ?mul y b" by (rule le_left_mult, simp_all!)
- have b: "?t1 = ?mul (?mul y (?add A dA)) x" by (simp! add: mult_matrix_assoc_simple[THEN sym])
- have assoc: "associative ?add" by (simp! add: combine_matrix_assoc)
- have r_distr: "r_distributive ?mul ?add"
- apply (rule r_distributive_matrix)
- by (simp! add: distributive_def)+
- have l_distr: "l_distributive ?mul ?add"
- apply (rule l_distributive_matrix)
- by (simp! add: distributive_def)+
- have c:"?mul y (?add A dA) = ?add (?mul y A) (?mul y dA)"
- by (insert r_distr, simp add: r_distributive_def)
- have d:"?mul (?add (?mul y A) (?mul y dA)) x = ?add (?mul (?mul y A) x) (?mul (?mul y dA) x)"
- by (insert l_distr, simp add: l_distributive_def)
- have e:"\<dots> = ?add (?mul c x) (?mul (?mul y dA) x)" by (simp!)
- from a b c d e show "?add (?mul c x) (?mul (?mul y dA) x) <= ?mul y b" by simp
-qed
-
-lemma linprog_by_duality_approx:
- assumes
- "(A + dA) * x <= (b::('a::lordered_ring) matrix)"
- "y * A = c"
- "0 <= y"
- shows
- "c * x + (y * dA) * x <= y * b"
-apply (simp add: times_matrix_def plus_matrix_def)
-apply (rule linprog_by_duality_approx_general)
-apply (simp_all)
-apply (simp_all add: associative_def add_assoc mult_assoc)
-apply (simp_all add: commutative_def add_commute)
-apply (simp_all add: distributive_def l_distributive_def r_distributive_def left_distrib right_distrib)
-apply (simp_all! add: plus_matrix_def times_matrix_def)
-apply (simp add: add_mono)
-apply (simp add: mult_left_mono)
-done
-
-lemma linprog_by_duality:
- assumes
- "A * x <= (b::('a::lordered_ring) matrix)"
- "y * A = c"
- "0 <= y"
- shows
- "c * x <= y * b"
-proof -
- have a:"(A + 0) * x <= b" by (simp!)
- have b:"c * x + (y*0)*x <= y * b" by (rule linprog_by_duality_approx, simp_all!)
- show "c * x <= y*b" by (insert b, simp)
-qed
-
-end
--- a/src/HOL/Matrix/Matrix.thy Fri Sep 03 10:27:05 2004 +0200
+++ b/src/HOL/Matrix/Matrix.thy Fri Sep 03 17:10:36 2004 +0200
@@ -21,7 +21,10 @@
times_matrix_def: "A * B == mult_matrix (op *) (op +) A B"
lemma is_meet_combine_matrix_meet: "is_meet (combine_matrix meet)"
-by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le)
+ by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le)
+
+lemma is_join_combine_matrix_join: "is_join (combine_matrix join)"
+ by (simp_all add: is_join_def le_matrix_def join_left_le join_right_le join_imp_le)
instance matrix :: (lordered_ab_group) lordered_ab_group_meet
proof
@@ -284,7 +287,20 @@
apply (auto)
done
+lemma Rep_minus[simp]: "Rep_matrix (-(A::_::lordered_ab_group)) x y = - (Rep_matrix A x y)"
+ by (simp add: minus_matrix_def)
+lemma join_matrix: "join (A::('a::lordered_ring) matrix) B = combine_matrix join A B"
+ apply (insert join_unique[of "(combine_matrix join)::('a matrix \<Rightarrow> 'a matrix \<Rightarrow> 'a matrix)", simplified is_join_combine_matrix_join])
+ apply (simp)
+ done
+lemma meet_matrix: "meet (A::('a::lordered_ring) matrix) B = combine_matrix meet A B"
+ apply (insert meet_unique[of "(combine_matrix meet)::('a matrix \<Rightarrow> 'a matrix \<Rightarrow> 'a matrix)", simplified is_meet_combine_matrix_meet])
+ apply (simp)
+ done
+
+lemma Rep_abs[simp]: "Rep_matrix (abs (A::_::lordered_ring)) x y = abs (Rep_matrix A x y)"
+ by (simp add: abs_lattice join_matrix)
end
--- a/src/HOL/Matrix/MatrixGeneral.thy Fri Sep 03 10:27:05 2004 +0200
+++ b/src/HOL/Matrix/MatrixGeneral.thy Fri Sep 03 17:10:36 2004 +0200
@@ -1392,7 +1392,44 @@
apply (rule le_foldseq)
by (auto)
+lemma spec2: "! j i. P j i \<Longrightarrow> P j i" by blast
+lemma neg_imp: "(\<not> Q \<longrightarrow> \<not> P) \<Longrightarrow> P \<longrightarrow> Q" by blast
+
lemma singleton_matrix_le[simp]: "(singleton_matrix j i a <= singleton_matrix j i b) = (a <= (b::_::order))"
by (auto simp add: le_matrix_def)
+lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))"
+ apply (auto)
+ apply (simp add: le_matrix_def)
+ apply (drule_tac j=j and i=i in spec2)
+ apply (simp)
+ apply (simp add: le_matrix_def)
+ done
+
+lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)"
+ apply (auto)
+ apply (simp add: le_matrix_def)
+ apply (drule_tac j=j and i=i in spec2)
+ apply (simp)
+ apply (simp add: le_matrix_def)
+ done
+
+lemma move_matrix_le_zero[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))"
+ apply (auto simp add: le_matrix_def neg_def)
+ apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
+ apply (auto)
+ done
+
+lemma move_matrix_zero_le[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)"
+ apply (auto simp add: le_matrix_def neg_def)
+ apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
+ apply (auto)
+ done
+
+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))"
+ apply (auto simp add: le_matrix_def neg_def)
+ apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
+ apply (auto)
+ done
+
end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/MatrixLP.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,83 @@
+use "MatrixLP_gensimp.ML";
+
+signature MATRIX_LP =
+sig
+ val lp_dual_estimate : string -> int -> thm
+ val simplify : thm -> thm
+end
+
+structure MatrixLP : MATRIX_LP =
+struct
+
+val _ = SimprocsCodegen.use_simprocs_code sg geninput
+
+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'
+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'''''
+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'''''
+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'''''
+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'''''
+val simp_sorted_sparse_matrix = simp_SparseMatrix_sorted_sparse_matrix
+val simp_sorted_spvec = simp_SparseMatrix_sorted_spvec
+
+fun single_arg simp ct = snd (simp (snd (Thm.dest_comb ct)))
+
+fun double_arg simp ct =
+ let
+ val (a, y) = Thm.dest_comb ct
+ val (_, x) = Thm.dest_comb a
+ in
+ snd (simp x y)
+ end
+
+fun spmc ct =
+ (case term_of ct of
+ ((Const ("SparseMatrix.le_spmat", _)) $ _) =>
+ single_arg simp_le_spmat ct
+ | ((Const ("SparseMatrix.add_spmat", _)) $ _) =>
+ single_arg simp_add_spmat ct
+ | ((Const ("SparseMatrix.diff_spmat", _)) $ _ $ _) =>
+ double_arg simp_diff_spmat ct
+ | ((Const ("SparseMatrix.abs_spmat", _)) $ _) =>
+ single_arg simp_abs_spmat ct
+ | ((Const ("SparseMatrix.mult_spmat", _)) $ _ $ _) =>
+ double_arg simp_mult_spmat ct
+ | ((Const ("SparseMatrix.sorted_sparse_matrix", _)) $ _) =>
+ single_arg simp_sorted_sparse_matrix ct
+ | ((Const ("SparseMatrix.sorted_spvec", ty)) $ _) =>
+ single_arg simp_sorted_spvec ct
+ | _ => raise Fail_conv)
+
+fun lp_dual_estimate lptfile prec =
+ let
+ val th = inst_real (thm "SparseMatrix.spm_linprog_dual_estimate_1")
+ val sg = sign_of (theory "MatrixLP")
+ val (y, (A1, A2), (c1, c2), b, r) =
+ let
+ open fspmlp
+ val l = load lptfile prec false
+ in
+ (y l, A l, c l, b l, r l)
+ end
+ fun var s x = (cterm_of sg (Var ((s,0), FloatSparseMatrixBuilder.real_spmatT)), x)
+ 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
+ in
+ th
+ end
+
+fun simplify th =
+ let
+ val simp_th = botc spmc (cprop_of th)
+ val th = strip_shyps (equal_elim simp_th th)
+ fun removeTrue th = removeTrue (implies_elim th TrueI) handle _ => th
+ in
+ removeTrue th
+ end
+
+fun float2real (x,y) = (Real.fromInt x) * (Math.pow (2.0, Real.fromInt y))
+
+end
+
+val result = ref ([]:thm list)
+
+fun run c = timeit (fn () => (result := [MatrixLP.simplify c]; ()))
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/MatrixLP.thy Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,4 @@
+theory MatrixLP = Float + SparseMatrix:
+
+end
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/MatrixLP_gensimp.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,35 @@
+open BasisLibrary;
+
+use "Cplex.ML";
+use "CplexMatrixConverter.ML";
+use "ExactFloatingPoint.ML";
+use "FloatArith.ML";
+use "FloatSparseMatrixBuilder.ML";
+use "fspmlp.ML";
+use "codegen_prep.ML";
+use "eq_codegen.ML";
+use "conv.ML";
+
+
+val th = theory "MatrixLP"
+val sg = sign_of th
+
+fun prep x = standard (mk_meta_eq x)
+
+fun inst_real thm = standard (Thm.instantiate ([(fst (hd (term_tvars (prop_of thm))),
+ ctyp_of sg HOLogic.realT)], []) thm);
+
+val spm_lp_dual_estimate_simp =
+ (thms "Float.arith_no_let") @
+ (map inst_real (thms "SparseMatrix.sparse_row_matrix_arith_simps")) @
+ (thms "SparseMatrix.sorted_sp_simps") @
+ [thm "Product_Type.fst_conv", thm "Product_Type.snd_conv"] @
+ (thms "SparseMatrix.boolarith")
+
+val simp_with = Simplifier.simplify (HOL_basic_ss addsimps [thm "SparseMatrix.if_case_eq", thm "Float.norm_0_1"])
+
+val spm_lp_dual_estimate_simp = map simp_with spm_lp_dual_estimate_simp
+
+val geninput = codegen_prep.prepare_thms spm_lp_dual_estimate_simp
+
+val _ = SimprocsCodegen.use_simprocs_code sg geninput
--- a/src/HOL/Matrix/ROOT.ML Fri Sep 03 10:27:05 2004 +0200
+++ b/src/HOL/Matrix/ROOT.ML Fri Sep 03 17:10:36 2004 +0200
@@ -1,10 +1,1 @@
-(* Title: HOL/Matrix/ROOT.ML
- ID: $Id$
- Author: Steven Obua
- License: 2004 Technische Universität München
-
-Theory of matrices with an application of matrix theory to linear
-programming.
-*)
-
-use_thy "LinProg";
+use_thy "MatrixLP"
--- a/src/HOL/Matrix/SparseMatrix.thy Fri Sep 03 10:27:05 2004 +0200
+++ b/src/HOL/Matrix/SparseMatrix.thy Fri Sep 03 17:10:36 2004 +0200
@@ -94,9 +94,62 @@
done
consts
+ abs_spvec :: "('a::lordered_ring) spvec \<Rightarrow> 'a spvec"
+ minus_spvec :: "('a::lordered_ring) spvec \<Rightarrow> 'a spvec"
smult_spvec :: "('a::lordered_ring) \<Rightarrow> 'a spvec \<Rightarrow> 'a spvec"
addmult_spvec :: "('a::lordered_ring) * 'a spvec * 'a spvec \<Rightarrow> 'a spvec"
+primrec
+ "minus_spvec [] = []"
+ "minus_spvec (a#as) = (fst a, -(snd a))#(minus_spvec as)"
+
+primrec
+ "abs_spvec [] = []"
+ "abs_spvec (a#as) = (fst a, abs (snd a))#(abs_spvec as)"
+
+lemma sparse_row_vector_minus:
+ "sparse_row_vector (minus_spvec v) = - (sparse_row_vector v)"
+ apply (induct v)
+ apply (simp_all add: sparse_row_vector_cons)
+ apply (simp add: Rep_matrix_inject[symmetric])
+ apply (rule ext)+
+ apply simp
+ done
+
+lemma sparse_row_vector_abs:
+ "sorted_spvec v \<Longrightarrow> sparse_row_vector (abs_spvec v) = abs (sparse_row_vector v)"
+ apply (induct v)
+ apply (simp_all add: sparse_row_vector_cons)
+ apply (frule_tac sorted_spvec_cons1, simp)
+ apply (simp only: Rep_matrix_inject[symmetric])
+ apply (rule ext)+
+ apply auto
+ apply (subgoal_tac "Rep_matrix (sparse_row_vector list) 0 a = 0")
+ apply (simp)
+ apply (rule sorted_sparse_row_vector_zero)
+ apply auto
+ done
+
+lemma sorted_spvec_minus_spvec:
+ "sorted_spvec v \<Longrightarrow> sorted_spvec (minus_spvec v)"
+ apply (induct v)
+ apply (simp)
+ apply (frule sorted_spvec_cons1, simp)
+ apply (simp add: sorted_spvec.simps)
+ apply (case_tac list)
+ apply (simp_all)
+ done
+
+lemma sorted_spvec_abs_spvec:
+ "sorted_spvec v \<Longrightarrow> sorted_spvec (abs_spvec v)"
+ apply (induct v)
+ apply (simp)
+ apply (frule sorted_spvec_cons1, simp)
+ apply (simp add: sorted_spvec.simps)
+ apply (case_tac list)
+ apply (simp_all)
+ done
+
defs
smult_spvec_def: "smult_spvec y arr == map (% a. (fst a, y * snd a)) arr"
@@ -542,9 +595,6 @@
else
(le_spvec(snd a, snd b) & le_spmat (as, bs))))"
-lemma spec2: "! j i. P j i \<Longrightarrow> P j i" by blast
-lemma neg_imp: "(\<not> Q \<longrightarrow> \<not> P) \<Longrightarrow> P \<longrightarrow> Q" by blast
-
constdefs
disj_matrices :: "('a::zero) matrix \<Rightarrow> 'a matrix \<Rightarrow> bool"
"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))"
@@ -597,22 +647,6 @@
(B + A <= C) = (A <= C & (B::('a::lordered_ab_group) matrix) <= 0)"
by (auto simp add: disj_matrices_add[of B A 0 C,simplified] disj_matrices_commute)
-lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))"
- apply (auto)
- apply (simp add: le_matrix_def)
- apply (drule_tac j=j and i=i in spec2)
- apply (simp)
- apply (simp add: le_matrix_def)
- done
-
-lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)"
- apply (auto)
- apply (simp add: le_matrix_def)
- apply (drule_tac j=j and i=i in spec2)
- apply (simp)
- apply (simp add: le_matrix_def)
- done
-
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)"
apply (simp add: disj_matrices_def)
apply (rule conjI)
@@ -641,27 +675,6 @@
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)"
by (auto simp add: disj_matrices_def)
-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)"
- apply (rule le_spvec.induct)
- 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)
- apply (rule conjI, intro strip)
- apply (simp add: sorted_spvec_cons1)
- apply (subst disj_matrices_add_x_le)
- apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute)
- apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
- apply (simp, blast)
- apply (intro strip, rule conjI, intro strip)
- apply (simp add: sorted_spvec_cons1)
- apply (subst disj_matrices_add_le_x)
- 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)
- apply (blast)
- apply (intro strip)
- apply (simp add: sorted_spvec_cons1)
- apply (case_tac "a=aa", simp_all)
- apply (subst disj_matrices_add)
- apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
- done
-
lemma disj_move_sparse_vec_mat[simplified disj_matrices_commute]:
"j <= a \<Longrightarrow> sorted_spvec((a,c)#as) \<Longrightarrow> disj_matrices (move_matrix (sparse_row_vector b) (int j) i) (sparse_row_matrix as)"
apply (auto simp add: neg_def disj_matrices_def)
@@ -685,24 +698,28 @@
apply (rule nrows, rule order_trans[of _ 1], simp, drule nrows_notzero, drule less_le_trans[OF _ nrows_spvec], arith)+
done
-lemma move_matrix_le_zero[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))"
- apply (auto simp add: le_matrix_def neg_def)
- apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
- apply (auto)
+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)"
+ apply (rule le_spvec.induct)
+ 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)
+ apply (rule conjI, intro strip)
+ apply (simp add: sorted_spvec_cons1)
+ apply (subst disj_matrices_add_x_le)
+ apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute)
+ apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
+ apply (simp, blast)
+ apply (intro strip, rule conjI, intro strip)
+ apply (simp add: sorted_spvec_cons1)
+ apply (subst disj_matrices_add_le_x)
+ 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)
+ apply (blast)
+ apply (intro strip)
+ apply (simp add: sorted_spvec_cons1)
+ apply (case_tac "a=aa", simp_all)
+ apply (subst disj_matrices_add)
+ apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
done
-lemma move_matrix_zero_le[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)"
- apply (auto simp add: le_matrix_def neg_def)
- apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
- apply (auto)
- done
-
-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))"
- apply (auto simp add: le_matrix_def neg_def)
- apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
- apply (auto)
- done
-
lemma le_spvec_empty2_sparse_row[rule_format]: "(sorted_spvec b) \<longrightarrow> (le_spvec (b,[]) = (sparse_row_vector b <= 0))"
apply (induct b)
apply (simp_all add: sorted_spvec_cons1)
@@ -752,39 +769,169 @@
apply (simp add: sorted_spvec_cons1 le_spvec_iff_sparse_row_le)
done
-term smult_spvec
-term addmult_spvec
-term add_spvec
-term mult_spvec_spmat
-term mult_spmat
-term add_spmat
-term le_spvec
-term le_spmat
-term sorted_spvec
-term sorted_spmat
+ML {* simp_depth_limit := 999 *}
+
+consts
+ abs_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat"
+ minus_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat"
+
+primrec
+ "abs_spmat [] = []"
+ "abs_spmat (a#as) = (fst a, abs_spvec (snd a))#(abs_spmat as)"
+
+primrec
+ "minus_spmat [] = []"
+ "minus_spmat (a#as) = (fst a, minus_spvec (snd a))#(minus_spmat as)"
+
+lemma sparse_row_matrix_minus:
+ "sparse_row_matrix (minus_spmat A) = - (sparse_row_matrix A)"
+ apply (induct A)
+ apply (simp_all add: sparse_row_vector_minus sparse_row_matrix_cons)
+ apply (subst Rep_matrix_inject[symmetric])
+ apply (rule ext)+
+ apply simp
+ done
-thm sparse_row_mult_spmat
-thm sparse_row_add_spmat
-thm le_spmat_iff_sparse_row_le
+lemma Rep_sparse_row_vector_zero: "x \<noteq> 0 \<Longrightarrow> Rep_matrix (sparse_row_vector v) x y = 0"
+proof -
+ assume x:"x \<noteq> 0"
+ have r:"nrows (sparse_row_vector v) <= Suc 0" by (rule nrows_spvec)
+ show ?thesis
+ apply (rule nrows)
+ apply (subgoal_tac "Suc 0 <= x")
+ apply (insert r)
+ apply (simp only:)
+ apply (insert x)
+ apply arith
+ done
+qed
+
+lemma sparse_row_matrix_abs:
+ "sorted_spvec A \<Longrightarrow> sorted_spmat A \<Longrightarrow> sparse_row_matrix (abs_spmat A) = abs (sparse_row_matrix A)"
+ apply (induct A)
+ apply (simp_all add: sparse_row_vector_abs sparse_row_matrix_cons)
+ apply (frule_tac sorted_spvec_cons1, simp)
+ apply (subst Rep_matrix_inject[symmetric])
+ apply (rule ext)+
+ apply auto
+ apply (case_tac "x=a")
+ apply (simp)
+ apply (subst sorted_sparse_row_matrix_zero)
+ apply auto
+ apply (subst Rep_sparse_row_vector_zero)
+ apply (simp_all add: neg_def)
+ done
+
+lemma sorted_spvec_minus_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec (minus_spmat A)"
+ apply (induct A)
+ apply (simp)
+ apply (frule sorted_spvec_cons1, simp)
+ apply (simp add: sorted_spvec.simps)
+ apply (case_tac list)
+ apply (simp_all)
+ done
+
+lemma sorted_spvec_abs_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec (abs_spmat A)"
+ apply (induct A)
+ apply (simp)
+ apply (frule sorted_spvec_cons1, simp)
+ apply (simp add: sorted_spvec.simps)
+ apply (case_tac list)
+ apply (simp_all)
+ done
+
+lemma sorted_spmat_minus_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat (minus_spmat A)"
+ apply (induct A)
+ apply (simp_all add: sorted_spvec_minus_spvec)
+ done
+
+lemma sorted_spmat_abs_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat (abs_spmat A)"
+ apply (induct A)
+ apply (simp_all add: sorted_spvec_abs_spvec)
+ done
-thm sorted_spvec_mult_spmat
-thm sorted_spmat_mult_spmat
-thm sorted_spvec_add_spmat
-thm sorted_spmat_add_spmat
+constdefs
+ diff_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat \<Rightarrow> 'a spmat"
+ "diff_spmat A B == add_spmat (A, minus_spmat B)"
+
+lemma sorted_spmat_diff_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat B \<Longrightarrow> sorted_spmat (diff_spmat A B)"
+ by (simp add: diff_spmat_def sorted_spmat_minus_spmat sorted_spmat_add_spmat)
+
+lemma sorted_spvec_diff_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec B \<Longrightarrow> sorted_spvec (diff_spmat A B)"
+ by (simp add: diff_spmat_def sorted_spvec_minus_spmat sorted_spvec_add_spmat)
+
+lemma sparse_row_diff_spmat: "sparse_row_matrix (diff_spmat A B ) = (sparse_row_matrix A) - (sparse_row_matrix B)"
+ by (simp add: diff_spmat_def sparse_row_add_spmat sparse_row_matrix_minus)
+
+constdefs
+ sorted_sparse_matrix :: "'a spmat \<Rightarrow> bool"
+ "sorted_sparse_matrix A == (sorted_spvec A) & (sorted_spmat A)"
+
+lemma sorted_sparse_matrix_imp_spvec: "sorted_sparse_matrix A \<Longrightarrow> sorted_spvec A"
+ by (simp add: sorted_sparse_matrix_def)
+
+lemma sorted_sparse_matrix_imp_spmat: "sorted_sparse_matrix A \<Longrightarrow> sorted_spmat A"
+ by (simp add: sorted_sparse_matrix_def)
+
+lemmas sparse_row_matrix_op_simps =
+ sorted_sparse_matrix_imp_spmat sorted_sparse_matrix_imp_spvec
+ sparse_row_add_spmat sorted_spvec_add_spmat sorted_spmat_add_spmat
+ sparse_row_diff_spmat sorted_spvec_diff_spmat sorted_spmat_diff_spmat
+ sparse_row_matrix_minus sorted_spvec_minus_spmat sorted_spmat_minus_spmat
+ sparse_row_mult_spmat sorted_spvec_mult_spmat sorted_spmat_mult_spmat
+ sparse_row_matrix_abs sorted_spvec_abs_spmat sorted_spmat_abs_spmat
+ le_spmat_iff_sparse_row_le
+
+lemma zero_eq_Numeral0: "(0::_::number_ring) = Numeral0" by simp
-thm smult_spvec_empty
-thm smult_spvec_cons
-thm addmult_spvec.simps
-thm add_spvec.simps
-thm add_spmat.simps
-thm mult_spvec_spmat.simps
-thm mult_spmat.simps
-thm le_spvec.simps
-thm le_spmat.simps
-thm sorted_spvec.simps
-thm sorted_spmat.simps
+lemmas sparse_row_matrix_arith_simps[simplified zero_eq_Numeral0] =
+ mult_spmat.simps mult_spvec_spmat.simps
+ addmult_spvec.simps
+ smult_spvec_empty smult_spvec_cons
+ add_spmat.simps add_spvec.simps
+ minus_spmat.simps minus_spvec.simps
+ abs_spmat.simps abs_spvec.simps
+ diff_spmat_def
+ le_spmat.simps le_spvec.simps
+
+lemmas sorted_sp_simps =
+ sorted_spvec.simps
+ sorted_spmat.simps
+ sorted_sparse_matrix_def
+
+lemma bool1: "(\<not> True) = False" by blast
+lemma bool2: "(\<not> False) = True" by blast
+lemma bool3: "((P\<Colon>bool) \<and> True) = P" by blast
+lemma bool4: "(True \<and> (P\<Colon>bool)) = P" by blast
+lemma bool5: "((P\<Colon>bool) \<and> False) = False" by blast
+lemma bool6: "(False \<and> (P\<Colon>bool)) = False" by blast
+lemma bool7: "((P\<Colon>bool) \<or> True) = True" by blast
+lemma bool8: "(True \<or> (P\<Colon>bool)) = True" by blast
+lemma bool9: "((P\<Colon>bool) \<or> False) = P" by blast
+lemma bool10: "(False \<or> (P\<Colon>bool)) = P" by blast
+lemmas boolarith = bool1 bool2 bool3 bool4 bool5 bool6 bool7 bool8 bool9 bool10
+
+lemma if_case_eq: "(if b then x else y) = (case b of True => x | False => y)" by simp
+
+lemma spm_linprog_dual_estimate_1:
+ assumes
+ "sorted_sparse_matrix A1"
+ "sorted_sparse_matrix A2"
+ "sorted_sparse_matrix c1"
+ "sorted_sparse_matrix c2"
+ "sorted_sparse_matrix y"
+ "sorted_spvec b"
+ "sorted_spvec r"
+ "le_spmat ([], y)"
+ "A * x \<le> sparse_row_matrix (b::('a::lordered_ring) spmat)"
+ "sparse_row_matrix A1 <= A"
+ "A <= sparse_row_matrix A2"
+ "sparse_row_matrix c1 <= c"
+ "c <= sparse_row_matrix c2"
+ "abs x \<le> sparse_row_matrix r"
+ shows
+ "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),
+ abs_spmat (diff_spmat (mult_spmat y A1) c1)), diff_spmat c2 c1)) r))"
+ by (insert prems, simp add: sparse_row_matrix_op_simps linprog_dual_estimate_1[where A=A])
end
-
-
-
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/codegen_prep.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,110 @@
+structure codegen_prep =
+struct
+
+exception Prepare_thms of string;
+
+fun is_meta_eq th =
+ let
+ fun check ((Const ("==", _)) $ _ $ _) = true
+ | check _ = false
+ in
+ check (concl_of th)
+ end
+
+fun printty ty =
+ let
+ fun immerse s [] = []
+ | immerse s [x] = [x]
+ | immerse s (x::xs) = x::s::(immerse s xs)
+
+ fun py t =
+ let
+ val (name, args) = dest_Type t
+ val args' = map printty args
+ in
+ concat (immerse "_" (name::args'))
+ end
+
+ val (args, res) = strip_type ty
+ val tystrs = map py (args@[res])
+
+ val s = "'"^(concat (immerse "_" tystrs))^"'"
+ in
+ s
+ end
+
+fun head_name th =
+ let val s =
+ let
+ val h = fst (strip_comb (hd (snd (strip_comb (concl_of th)))))
+ val n = fst (dest_Const h)
+ val ty = snd (dest_Const h)
+ val hn = fst (dest_Const h)
+ in
+ hn^"_"^(printty ty) handle _ => (writeln ("warning: polymorphic "^hn); hn)
+ end
+ in
+ s
+ end
+
+val mangle_id =
+ let
+ fun tr #"=" = "_eq_"
+ | tr #"." = "_"
+ | tr #" " = "_"
+ | tr #"<" = "_le_"
+ | tr #">" = "_ge_"
+ | tr #"(" = "_bro_"
+ | tr #")" = "_brc_"
+ | tr #"+" = "_plus_"
+ | tr #"*" = "_mult_"
+ | tr #"-" = "_minus_"
+ | tr #"|" = "_or_"
+ | tr #"&" = "_and_"
+ | tr x = Char.toString x
+ fun m s = "simp_"^(String.translate tr s)
+ in
+ m
+ end
+
+fun group f [] = []
+ | group f (x::xs) =
+ let
+ val g = group f xs
+ val key = f x
+ in
+ case assoc (g, key) of
+ None => (key, [x])::g
+ | Some l => overwrite (g, (key, x::l))
+ end
+
+fun prepare_thms ths =
+ let
+ val ths = (filter is_meta_eq ths) @ (map (standard o mk_meta_eq) (filter (not o is_meta_eq) ths))
+ val _ = if forall Thm.no_prems ths then () else raise (Prepare_thms "premisse found")
+ val thmgroups = group head_name ths
+ val test_clashgroups = group (fn (a,b) => mangle_id a) thmgroups
+ val _ = if (length thmgroups <> length test_clashgroups) then raise (Prepare_thms "clash after name mangling") else ()
+
+ fun prep (name, ths) =
+ let
+ val m = mangle_id name
+
+ in
+ (m, true, ths)
+ end
+
+ val thmgroups = map prep thmgroups
+ in
+ (thmgroups)
+ end
+
+fun writestr filename s =
+ let
+ val f = TextIO.openOut filename
+ val _ = TextIO.output(f, s)
+ val _ = TextIO.closeOut f
+ in
+ ()
+ end
+end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/conv.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,26 @@
+exception Fail_conv;
+
+fun orelsec conv1 conv2 ct = conv1 ct handle Fail_conv => conv2 ct
+
+val allc = Thm.reflexive
+
+fun thenc conv1 conv2 ct =
+ let
+ fun rhs_of t = snd (Thm.dest_comb (strip_imp_concl (cprop_of t)))
+
+ val eq = conv1 ct
+ in
+ Thm.transitive eq (conv2 (rhs_of eq))
+ end
+
+fun subc conv ct =
+ case term_of ct of
+ _ $ _ =>
+ let
+ val (ct1, ct2) = Thm.dest_comb ct
+ in
+ Thm.combination (conv ct1) (conv ct2)
+ end
+ | _ => allc ct
+
+fun botc conv ct = thenc (subc (botc conv)) (orelsec conv allc) ct
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/eq_codegen.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,493 @@
+fun inst_cterm inst ct = fst (Drule.dest_equals
+ (Thm.cprop_of (Thm.instantiate inst (reflexive ct))));
+fun tyinst_cterm tyinst = inst_cterm (tyinst, []);
+
+val bla = ref ([] : term list);
+
+(******************************************************)
+(* Code generator for equational proofs *)
+(******************************************************)
+fun my_mk_meta_eq thm =
+ let
+ val (_, eq) = Thm.dest_comb (cprop_of thm);
+ val (ct, rhs) = Thm.dest_comb eq;
+ val (_, lhs) = Thm.dest_comb ct
+ in Thm.implies_elim (Drule.instantiate' [Some (ctyp_of_term lhs)]
+ [Some lhs, Some rhs] eq_reflection) thm
+ end;
+
+structure SimprocsCodegen =
+struct
+
+val simp_thms = ref ([] : thm list);
+
+fun parens b = if b then Pretty.enclose "(" ")" else Pretty.block;
+
+fun gen_mk_val f xs ps = Pretty.block ([Pretty.str "val ",
+ f (length xs > 1) (flat
+ (separate [Pretty.str ",", Pretty.brk 1] (map (single o Pretty.str) xs))),
+ Pretty.str " =", Pretty.brk 1] @ ps @ [Pretty.str ";"]);
+
+val mk_val = gen_mk_val parens;
+val mk_vall = gen_mk_val (K (Pretty.enclose "[" "]"));
+
+fun rename s = if s mem ThmDatabase.ml_reserved then s ^ "'" else s;
+
+fun mk_decomp_name (Var ((s, i), _)) = rename (if i=0 then s else s ^ string_of_int i)
+ | mk_decomp_name (Const (s, _)) = rename (Codegen.mk_id (Sign.base_name s))
+ | mk_decomp_name _ = "ct";
+
+fun decomp_term_code cn ((vs, bs, ps), (v, t)) =
+ if exists (equal t o fst) bs then (vs, bs, ps)
+ else (case t of
+ Var _ => (vs, bs @ [(t, v)], ps)
+ | Const _ => (vs, if cn then bs @ [(t, v)] else bs, ps)
+ | Bound _ => (vs, bs, ps)
+ | Abs (s, T, t) =>
+ let
+ val v1 = variant vs s;
+ val v2 = variant (v1 :: vs) (mk_decomp_name t)
+ in
+ decomp_term_code cn ((v1 :: v2 :: vs,
+ bs @ [(Free (s, T), v1)],
+ ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_abs", Pretty.brk 1,
+ Pretty.str "None", Pretty.brk 1, Pretty.str v]]), (v2, t))
+ end
+ | t $ u =>
+ let
+ val v1 = variant vs (mk_decomp_name t);
+ val v2 = variant (v1 :: vs) (mk_decomp_name u);
+ val (vs', bs', ps') = decomp_term_code cn ((v1 :: v2 :: vs, bs,
+ ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_comb", Pretty.brk 1,
+ Pretty.str v]]), (v1, t));
+ val (vs'', bs'', ps'') = decomp_term_code cn ((vs', bs', ps'), (v2, u))
+ in
+ if bs'' = bs then (vs, bs, ps) else (vs'', bs'', ps'')
+ end);
+
+val strip_tv = implode o tl o explode;
+
+fun mk_decomp_tname (TVar ((s, i), _)) =
+ strip_tv ((if i=0 then s else s ^ string_of_int i) ^ "T")
+ | mk_decomp_tname (Type (s, _)) = Codegen.mk_id (Sign.base_name s) ^ "T"
+ | mk_decomp_tname _ = "cT";
+
+fun decomp_type_code ((vs, bs, ps), (v, TVar (ixn, _))) =
+ if exists (equal ixn o fst) bs then (vs, bs, ps)
+ else (vs, bs @ [(ixn, v)], ps)
+ | decomp_type_code ((vs, bs, ps), (v, Type (_, Ts))) =
+ let
+ val vs' = variantlist (map mk_decomp_tname Ts, vs);
+ val (vs'', bs', ps') =
+ foldl decomp_type_code ((vs @ vs', bs, ps @
+ [mk_vall vs' [Pretty.str "Thm.dest_ctyp", Pretty.brk 1,
+ Pretty.str v]]), vs' ~~ Ts)
+ in
+ if bs' = bs then (vs, bs, ps) else (vs'', bs', ps')
+ end;
+
+fun gen_mk_bindings s dest decomp ((vs, bs, ps), (v, x)) =
+ let
+ val s' = variant vs s;
+ val (vs', bs', ps') = decomp ((s' :: vs, bs, ps @
+ [mk_val [s'] (dest v)]), (s', x))
+ in
+ if bs' = bs then (vs, bs, ps) else (vs', bs', ps')
+ end;
+
+val mk_term_bindings = gen_mk_bindings "ct"
+ (fn s => [Pretty.str "cprop_of", Pretty.brk 1, Pretty.str s])
+ (decomp_term_code true);
+
+val mk_type_bindings = gen_mk_bindings "cT"
+ (fn s => [Pretty.str "Thm.ctyp_of_term", Pretty.brk 1, Pretty.str s])
+ decomp_type_code;
+
+fun pretty_pattern b (Const (s, _)) = Pretty.block [Pretty.str "Const",
+ Pretty.brk 1, Pretty.str ("(\"" ^ s ^ "\", _)")]
+ | pretty_pattern b (t as _ $ _) = parens b
+ (flat (separate [Pretty.str " $", Pretty.brk 1]
+ (map (single o pretty_pattern true) (op :: (strip_comb t)))))
+ | pretty_pattern b _ = Pretty.str "_";
+
+fun term_consts' t = foldl_aterms
+ (fn (cs, c as Const _) => c ins cs | (cs, _) => cs) ([], t);
+
+fun mk_apps s b p [] = p
+ | mk_apps s b p (q :: qs) =
+ mk_apps s b (parens (b orelse not (null qs))
+ [Pretty.str s, Pretty.brk 1, p, Pretty.brk 1, q]) qs;
+
+fun mk_refleq eq ct = mk_val [eq] [Pretty.str ("Thm.reflexive " ^ ct)];
+
+fun mk_tyinst ((s, i), s') =
+ Pretty.block [Pretty.str ("((" ^ quote s ^ ","), Pretty.brk 1,
+ Pretty.str (string_of_int i ^ "),"), Pretty.brk 1,
+ Pretty.str (s' ^ ")")];
+
+fun inst_ty b ty_bs t s = (case term_tvars t of
+ [] => Pretty.str s
+ | Ts => parens b [Pretty.str "tyinst_cterm", Pretty.brk 1,
+ Pretty.list "[" "]" (map (fn (ixn, _) => mk_tyinst
+ (ixn, the (assoc (ty_bs, ixn)))) Ts),
+ Pretty.brk 1, Pretty.str s]);
+
+fun mk_cterm_code b ty_bs ts xs (vals, t $ u) =
+ let
+ val (vals', p1) = mk_cterm_code true ty_bs ts xs (vals, t);
+ val (vals'', p2) = mk_cterm_code true ty_bs ts xs (vals', u)
+ in
+ (vals'', parens b [Pretty.str "Thm.capply", Pretty.brk 1,
+ p1, Pretty.brk 1, p2])
+ end
+ | mk_cterm_code b ty_bs ts xs (vals, Abs (s, T, t)) =
+ let
+ val u = Free (s, T);
+ val Some s' = assoc (ts, u);
+ val p = Pretty.str s';
+ val (vals', p') = mk_cterm_code true ty_bs ts (p :: xs)
+ (if null (typ_tvars T) then vals
+ else vals @ [(u, (("", s'), [mk_val [s'] [inst_ty true ty_bs u s']]))], t)
+ in (vals',
+ parens b [Pretty.str "Thm.cabs", Pretty.brk 1, p, Pretty.brk 1, p'])
+ end
+ | mk_cterm_code b ty_bs ts xs (vals, Bound i) = (vals, nth_elem (i, xs))
+ | mk_cterm_code b ty_bs ts xs (vals, t) = (case assoc (vals, t) of
+ None =>
+ let val Some s = assoc (ts, t)
+ in (if is_Const t andalso not (null (term_tvars t)) then
+ vals @ [(t, (("", s), [mk_val [s] [inst_ty true ty_bs t s]]))]
+ else vals, Pretty.str s)
+ end
+ | Some ((_, s), _) => (vals, Pretty.str s));
+
+fun get_cases sg =
+ Symtab.foldl (fn (tab, (k, {case_rewrites, ...})) => Symtab.update_new
+ ((fst (dest_Const (head_of (fst (HOLogic.dest_eq (HOLogic.dest_Trueprop
+ (prop_of (hd case_rewrites))))))), map my_mk_meta_eq case_rewrites), tab))
+ (Symtab.empty, DatatypePackage.get_datatypes_sg sg);
+
+fun decomp_case th =
+ let
+ val (lhs, _) = Logic.dest_equals (prop_of th);
+ val (f, ts) = strip_comb lhs;
+ val (us, u) = split_last ts;
+ val (Const (s, _), vs) = strip_comb u
+ in (us, s, vs, u) end;
+
+fun rename vs t =
+ let
+ fun mk_subst ((vs, subs), Var ((s, i), T)) =
+ let val s' = variant vs s
+ in if s = s' then (vs, subs)
+ else (s' :: vs, ((s, i), Var ((s', i), T)) :: subs)
+ end;
+ val (vs', subs) = foldl mk_subst ((vs, []), term_vars t)
+ in (vs', subst_Vars subs t) end;
+
+fun is_instance sg t u = t = subst_TVars_Vartab
+ (Type.typ_match (Sign.tsig_of sg) (Vartab.empty,
+ (fastype_of u, fastype_of t))) u handle Type.TYPE_MATCH => false;
+
+(*
+fun lookup sg fs t = apsome snd (Library.find_first
+ (is_instance sg t o fst) fs);
+*)
+
+fun lookup sg fs t = (case Library.find_first (is_instance sg t o fst) fs of
+ None => (bla := (t ins !bla); None)
+ | Some (_, x) => Some x);
+
+fun unint sg fs t = forall (is_none o lookup sg fs) (term_consts' t);
+
+fun mk_let s i xs ys =
+ Pretty.blk (0, [Pretty.blk (i, separate Pretty.fbrk (Pretty.str s :: xs)),
+ Pretty.fbrk,
+ Pretty.blk (i, ([Pretty.str "in", Pretty.fbrk] @ ys)),
+ Pretty.fbrk, Pretty.str "end"]);
+
+(*****************************************************************************)
+(* Generate bindings for simplifying term t *)
+(* mkeq: whether to generate reflexivity theorem for uninterpreted terms *)
+(* fs: interpreted functions *)
+(* ts: atomic terms *)
+(* vs: used identifiers *)
+(* vals: list of bindings of the form ((eq, ct), ps) where *)
+(* eq: name of equational theorem *)
+(* ct: name of simplified cterm *)
+(* ps: ML code for creating the above two items *)
+(*****************************************************************************)
+
+fun mk_simpl_code sg case_tab mkeq fs ts ty_bs thm_bs ((vs, vals), t) =
+ (case assoc (vals, t) of
+ Some ((eq, ct), ps) => (* binding already generated *)
+ if mkeq andalso eq="" then
+ let val eq' = variant vs "eq"
+ in ((eq' :: vs, overwrite (vals,
+ (t, ((eq', ct), ps @ [mk_refleq eq' ct])))), (eq', ct))
+ end
+ else ((vs, vals), (eq, ct))
+ | None => (case assoc (ts, t) of
+ Some v => (* atomic term *)
+ let val xs = if not (null (term_tvars t)) andalso is_Const t then
+ [mk_val [v] [inst_ty false ty_bs t v]] else []
+ in
+ if mkeq then
+ let val eq = variant vs "eq"
+ in ((eq :: vs, vals @
+ [(t, ((eq, v), xs @ [mk_refleq eq v]))]), (eq, v))
+ end
+ else ((vs, if null xs then vals else vals @
+ [(t, (("", v), xs))]), ("", v))
+ end
+ | None => (* complex term *)
+ let val (f as Const (cname, _), us) = strip_comb t
+ in case Symtab.lookup (case_tab, cname) of
+ Some cases => (* case expression *)
+ let
+ val (us', u) = split_last us;
+ val b = unint sg fs u;
+ val ((vs1, vals1), (eq, ct)) =
+ mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs ((vs, vals), u);
+ val xs = variantlist (replicate (length us') "f", vs1);
+ val (vals2, ps) = foldl_map
+ (mk_cterm_code false ty_bs ts []) (vals1, us');
+ val fvals = map (fn (x, p) => mk_val [x] [p]) (xs ~~ ps);
+ val uT = fastype_of u;
+ val (us'', _, _, u') = decomp_case (hd cases);
+ val (vs2, ty_bs', ty_vals) = mk_type_bindings
+ (mk_type_bindings ((vs1 @ xs, [], []),
+ (hd xs, fastype_of (hd us''))), (ct, fastype_of u'));
+ val insts1 = map mk_tyinst ty_bs';
+ val i = length vals2;
+
+ fun mk_case_code ((vs, vals), (f, (name, eqn))) =
+ let
+ val (fvs, cname, cvs, _) = decomp_case eqn;
+ val Ts = binder_types (fastype_of f);
+ val ys = variantlist (map (fst o fst o dest_Var) cvs, vs);
+ val cvs' = map Var (map (rpair 0) ys ~~ Ts);
+ val rs = cvs' ~~ cvs;
+ val lhs = list_comb (Const (cname, Ts ---> uT), cvs');
+ val rhs = foldl betapply (f, cvs');
+ val (vs', tm_bs, tm_vals) = decomp_term_code false
+ ((vs @ ys, [], []), (ct, lhs));
+ val ((vs'', all_vals), (eq', ct')) = mk_simpl_code sg case_tab
+ false fs (tm_bs @ ts) ty_bs thm_bs ((vs', vals), rhs);
+ val (old_vals, eq_vals) = splitAt (i, all_vals);
+ val vs''' = vs @ filter (fn v => exists
+ (fn (_, ((v', _), _)) => v = v') old_vals) (vs'' \\ vs');
+ val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(",
+ inst_ty false ty_bs' t (the (assoc (thm_bs, t))), Pretty.str ",",
+ Pretty.brk 1, Pretty.str (s ^ ")")]) ((fvs ~~ xs) @
+ (map (fn (v, s) => (the (assoc (rs, v)), s)) tm_bs));
+ val eq'' = if null insts1 andalso null insts2 then Pretty.str name
+ else parens (eq' <> "") [Pretty.str
+ (if null cvs then "Thm.instantiate" else "Drule.instantiate"),
+ Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1,
+ Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2,
+ Pretty.str ")", Pretty.brk 1, Pretty.str name];
+ val eq''' = if eq' = "" then eq'' else
+ Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1,
+ eq'', Pretty.brk 1, Pretty.str eq']
+ in
+ ((vs''', old_vals), Pretty.block [pretty_pattern false lhs,
+ Pretty.str " =>",
+ Pretty.brk 1, mk_let "let" 2 (tm_vals @ flat (map (snd o snd) eq_vals))
+ [Pretty.str ("(" ^ ct' ^ ","), Pretty.brk 1, eq''', Pretty.str ")"]])
+ end;
+
+ val case_names = map (fn i => Sign.base_name cname ^ "_" ^
+ string_of_int i) (1 upto length cases);
+ val ((vs3, vals3), case_ps) = foldl_map mk_case_code
+ ((vs2, vals2), us' ~~ (case_names ~~ cases));
+ val eq' = variant vs3 "eq";
+ val ct' = variant (eq' :: vs3) "ct";
+ val eq'' = variant (eq' :: ct' :: vs3) "eq";
+ val case_vals =
+ fvals @ ty_vals @
+ [mk_val [ct', eq'] ([Pretty.str "(case", Pretty.brk 1,
+ Pretty.str ("term_of " ^ ct ^ " of"), Pretty.brk 1] @
+ flat (separate [Pretty.brk 1, Pretty.str "| "]
+ (map single case_ps)) @ [Pretty.str ")"])]
+ in
+ if b then
+ ((eq' :: ct' :: vs3, vals3 @
+ [(t, ((eq', ct'), case_vals))]), (eq', ct'))
+ else
+ let val ((vs4, vals4), (_, ctcase)) = mk_simpl_code sg case_tab false
+ fs ts ty_bs thm_bs ((eq' :: eq'' :: ct' :: vs3, vals3), f)
+ in
+ ((vs4, vals4 @ [(t, ((eq'', ct'), case_vals @
+ [mk_val [eq''] [Pretty.str "Thm.transitive", Pretty.brk 1,
+ Pretty.str "(Thm.combination", Pretty.brk 1,
+ Pretty.str "(Thm.reflexive", Pretty.brk 1,
+ mk_apps "Thm.capply" true (Pretty.str ctcase)
+ (map Pretty.str xs),
+ Pretty.str ")", Pretty.brk 1, Pretty.str (eq ^ ")"),
+ Pretty.brk 1, Pretty.str eq']]))]), (eq'', ct'))
+ end
+ end
+
+ | None =>
+ let
+ val b = forall (unint sg fs) us;
+ val (q, eqs) = foldl_map
+ (mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs) ((vs, vals), us);
+ val ((vs', vals'), (eqf, ctf)) = if is_some (lookup sg fs f) andalso b
+ then (q, ("", ""))
+ else mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs (q, f);
+ val ct = variant vs' "ct";
+ val eq = variant (ct :: vs') "eq";
+ val ctv = mk_val [ct] [mk_apps "Thm.capply" false
+ (Pretty.str ctf) (map (Pretty.str o snd) eqs)];
+ fun combp b = mk_apps "Thm.combination" b
+ (Pretty.str eqf) (map (Pretty.str o fst) eqs)
+ in
+ case (lookup sg fs f, b) of
+ (None, true) => (* completely uninterpreted *)
+ if mkeq then ((ct :: eq :: vs', vals' @
+ [(t, ((eq, ct), [ctv, mk_refleq eq ct]))]), (eq, ct))
+ else ((ct :: vs', vals' @ [(t, (("", ct), [ctv]))]), ("", ct))
+ | (None, false) => (* function uninterpreted *)
+ ((eq :: ct :: vs', vals' @
+ [(t, ((eq, ct), [ctv, mk_val [eq] [combp false]]))]), (eq, ct))
+ | (Some (s, _, _), true) => (* arguments uninterpreted *)
+ ((eq :: ct :: vs', vals' @
+ [(t, ((eq, ct), [mk_val [ct, eq] (separate (Pretty.brk 1)
+ (Pretty.str s :: map (Pretty.str o snd) eqs))]))]), (eq, ct))
+ | (Some (s, _, _), false) => (* function and arguments interpreted *)
+ let val eq' = variant (eq :: ct :: vs') "eq"
+ in ((eq' :: eq :: ct :: vs', vals' @ [(t, ((eq', ct),
+ [mk_val [ct, eq] (separate (Pretty.brk 1)
+ (Pretty.str s :: map (Pretty.str o snd) eqs)),
+ mk_val [eq'] [Pretty.str "Thm.transitive", Pretty.brk 1,
+ combp true, Pretty.brk 1, Pretty.str eq]]))]), (eq', ct))
+ end
+ end
+ end));
+
+fun lhs_of thm = fst (Logic.dest_equals (prop_of thm));
+fun rhs_of thm = snd (Logic.dest_equals (prop_of thm));
+
+fun mk_funs_code sg case_tab fs fs' =
+ let
+ val case_thms = mapfilter (fn s => (case Symtab.lookup (case_tab, s) of
+ None => None
+ | Some thms => Some (unsuffix "_case" (Sign.base_name s) ^ ".cases",
+ map (fn i => Sign.base_name s ^ "_" ^ string_of_int i)
+ (1 upto length thms) ~~ thms)))
+ (foldr add_term_consts (map (prop_of o snd)
+ (flat (map (#3 o snd) fs')), []));
+ val case_vals = map (fn (s, cs) => mk_vall (map fst cs)
+ [Pretty.str "map my_mk_meta_eq", Pretty.brk 1,
+ Pretty.str ("(thms \"" ^ s ^ "\")")]) case_thms;
+ val (vs, thm_bs, thm_vals) = foldl mk_term_bindings (([], [], []),
+ flat (map (map (apsnd prop_of) o #3 o snd) fs') @
+ map (apsnd prop_of) (flat (map snd case_thms)));
+
+ fun mk_fun_code (prfx, (fname, d, eqns)) =
+ let
+ val (f, ts) = strip_comb (lhs_of (snd (hd eqns)));
+ val args = variantlist (replicate (length ts) "ct", vs);
+ val (vs', ty_bs, ty_vals) = foldl mk_type_bindings
+ ((vs @ args, [], []), args ~~ map fastype_of ts);
+ val insts1 = map mk_tyinst ty_bs;
+
+ fun mk_eqn_code (name, eqn) =
+ let
+ val (_, argts) = strip_comb (lhs_of eqn);
+ val (vs'', tm_bs, tm_vals) = foldl (decomp_term_code false)
+ ((vs', [], []), args ~~ argts);
+ val ((vs''', eq_vals), (eq, ct)) = mk_simpl_code sg case_tab false fs
+ (tm_bs @ filter_out (is_Var o fst) thm_bs) ty_bs thm_bs
+ ((vs'', []), rhs_of eqn);
+ val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(",
+ inst_ty false ty_bs t (the (assoc (thm_bs, t))), Pretty.str ",", Pretty.brk 1,
+ Pretty.str (s ^ ")")]) tm_bs
+ val eq' = if null insts1 andalso null insts2 then Pretty.str name
+ else parens (eq <> "") [Pretty.str "Thm.instantiate",
+ Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1,
+ Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2,
+ Pretty.str ")", Pretty.brk 1, Pretty.str name];
+ val eq'' = if eq = "" then eq' else
+ Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1,
+ eq', Pretty.brk 1, Pretty.str eq]
+ in
+ Pretty.block [parens (length argts > 1)
+ (Pretty.commas (map (pretty_pattern false) argts)),
+ Pretty.str " =>",
+ Pretty.brk 1, mk_let "let" 2 (ty_vals @ tm_vals @ flat (map (snd o snd) eq_vals))
+ [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1, eq'', Pretty.str ")"]]
+ end;
+
+ val default = if d then
+ let
+ val Some s = assoc (thm_bs, f);
+ val ct = variant vs' "ct"
+ in [Pretty.brk 1, Pretty.str "handle", Pretty.brk 1,
+ Pretty.str "Match =>", Pretty.brk 1, mk_let "let" 2
+ (ty_vals @ (if null (term_tvars f) then [] else
+ [mk_val [s] [inst_ty false ty_bs f s]]) @
+ [mk_val [ct] [mk_apps "Thm.capply" false (Pretty.str s)
+ (map Pretty.str args)]])
+ [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1,
+ Pretty.str "Thm.reflexive", Pretty.brk 1, Pretty.str (ct ^ ")")]]
+ end
+ else []
+ in
+ ("and ", Pretty.block (separate (Pretty.brk 1)
+ (Pretty.str (prfx ^ fname) :: map Pretty.str args) @
+ [Pretty.str " =", Pretty.brk 1, Pretty.str "(case", Pretty.brk 1,
+ Pretty.list "(" ")" (map (fn s => Pretty.str ("term_of " ^ s)) args),
+ Pretty.str " of", Pretty.brk 1] @
+ flat (separate [Pretty.brk 1, Pretty.str "| "]
+ (map (single o mk_eqn_code) eqns)) @ [Pretty.str ")"] @ default))
+ end;
+
+ val (_, decls) = foldl_map mk_fun_code ("fun ", map snd fs')
+ in
+ mk_let "local" 2 (case_vals @ thm_vals) (separate Pretty.fbrk decls)
+ end;
+
+fun mk_simprocs_code sg eqns =
+ let
+ val case_tab = get_cases sg;
+ fun get_head th = head_of (fst (Logic.dest_equals (prop_of th)));
+ fun attach_term (x as (_, _, (_, th) :: _)) = (get_head th, x);
+ val eqns' = map attach_term eqns;
+ fun mk_node (s, _, (_, th) :: _) = (s, get_head th);
+ fun mk_edges (s, _, ths) = map (pair s) (distinct
+ (mapfilter (fn t => apsome #1 (lookup sg eqns' t))
+ (flat (map (term_consts' o prop_of o snd) ths))));
+ val gr = foldr (uncurry Graph.add_edge)
+ (map (pair "" o #1) eqns @ flat (map mk_edges eqns),
+ foldr (uncurry Graph.new_node)
+ (("", Bound 0) :: map mk_node eqns, Graph.empty));
+ val keys = rev (Graph.all_succs gr [""] \ "");
+ fun gr_ord (x :: _, y :: _) =
+ int_ord (find_index (equal x) keys, find_index (equal y) keys);
+ val scc = map (fn xs => filter (fn (_, (s, _, _)) => s mem xs) eqns')
+ (sort gr_ord (Graph.strong_conn gr \ [""]));
+ in
+ flat (separate [Pretty.str ";", Pretty.fbrk, Pretty.str " ", Pretty.fbrk]
+ (map (fn eqns'' => [mk_funs_code sg case_tab eqns' eqns'']) scc)) @
+ [Pretty.str ";", Pretty.fbrk]
+ end;
+
+fun use_simprocs_code sg eqns =
+ let
+ fun attach_name (i, x) = (i+1, ("simp_thm_" ^ string_of_int i, x));
+ fun attach_names (i, (s, b, eqs)) =
+ let val (i', eqs') = foldl_map attach_name (i, eqs)
+ in (i', (s, b, eqs')) end;
+ val (_, eqns') = foldl_map attach_names (1, eqns);
+ val (names, thms) = split_list (flat (map #3 eqns'));
+ val s = setmp print_mode [] Pretty.string_of
+ (mk_let "local" 2 [mk_vall names [Pretty.str "!SimprocsCodegen.simp_thms"]]
+ (mk_simprocs_code sg eqns'))
+ in
+ (simp_thms := thms; use_text Context.ml_output false s)
+ end;
+
+end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/fspmlp.ML Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,303 @@
+signature FSPMLP =
+sig
+ type linprog
+
+ val y : linprog -> cterm
+ val A : linprog -> cterm * cterm
+ val b : linprog -> cterm
+ val c : linprog -> cterm * cterm
+ val r : linprog -> cterm
+
+ exception Load of string
+
+ val load : string -> int -> bool -> linprog
+end
+
+structure fspmlp : FSPMLP =
+struct
+
+type linprog = cterm * (cterm * cterm) * cterm * (cterm * cterm) * cterm
+
+fun y (c1, c2, c3, c4, c5) = c1
+fun A (c1, c2, c3, c4, c5) = c2
+fun b (c1, c2, c3, c4, c5) = c3
+fun c (c1, c2, c3, c4, c5) = c4
+fun r (c1, c2, c3, c4, c5) = c5
+
+structure CplexFloatSparseMatrixConverter =
+MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder);
+
+datatype bound_type = LOWER | UPPER
+
+fun intbound_ord ((i1, b1),(i2,b2)) =
+ if i1 < i2 then LESS
+ else if i1 = i2 then
+ (if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER)
+ else GREATER
+
+structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord));
+
+structure VarGraph = TableFun(type key = int*bound_type val ord = intbound_ord);
+(* key -> (float option) * (int -> (float * (((float * float) * key) list)))) *)
+(* dest_key -> (sure_bound * (row_index -> (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *)
+
+exception Internal of string;
+
+fun add_row_bound g dest_key row_index row_bound =
+ let
+ val x =
+ case VarGraph.lookup (g, dest_key) of
+ None => (None, Inttab.update ((row_index, (row_bound, [])), Inttab.empty))
+ | Some (sure_bound, f) =>
+ (sure_bound,
+ case Inttab.lookup (f, row_index) of
+ None => Inttab.update ((row_index, (row_bound, [])), f)
+ | Some _ => raise (Internal "add_row_bound"))
+ in
+ VarGraph.update ((dest_key, x), g)
+ end
+
+fun update_sure_bound g (key as (_, btype)) bound =
+ let
+ val x =
+ case VarGraph.lookup (g, key) of
+ None => (Some bound, Inttab.empty)
+ | Some (None, f) => (Some bound, f)
+ | Some (Some old_bound, f) =>
+ (Some ((case btype of
+ UPPER => FloatArith.min
+ | LOWER => FloatArith.max)
+ old_bound bound), f)
+ in
+ VarGraph.update ((key, x), g)
+ end
+
+fun get_sure_bound g key =
+ case VarGraph.lookup (g, key) of
+ None => None
+ | Some (sure_bound, _) => sure_bound
+
+(*fun get_row_bound g key row_index =
+ case VarGraph.lookup (g, key) of
+ None => None
+ | Some (sure_bound, f) =>
+ (case Inttab.lookup (f, row_index) of
+ None => None
+ | Some (row_bound, _) => (sure_bound, row_bound))*)
+
+fun add_edge g src_key dest_key row_index coeff =
+ case VarGraph.lookup (g, dest_key) of
+ None => raise (Internal "add_edge: dest_key not found")
+ | Some (sure_bound, f) =>
+ (case Inttab.lookup (f, row_index) of
+ None => raise (Internal "add_edge: row_index not found")
+ | Some (row_bound, sources) =>
+ VarGraph.update ((dest_key, (sure_bound, Inttab.update ((row_index, (row_bound, (coeff, src_key) :: sources)), f))), g))
+
+fun split_graph g =
+ let
+ fun split ((r1, r2), (key, (sure_bound, _))) =
+ case sure_bound of
+ None => (r1, r2)
+ | Some bound =>
+ (case key of
+ (u, UPPER) => (r1, Inttab.update ((u, bound), r2))
+ | (u, LOWER) => (Inttab.update ((u, bound), r1), r2))
+ in
+ VarGraph.foldl split ((Inttab.empty, Inttab.empty), g)
+ end
+
+fun it2list t =
+ let
+ fun tolist (l, a) = a::l
+ in
+ Inttab.foldl tolist ([], t)
+ end
+
+(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm).
+ If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *)
+fun propagate_sure_bounds safe names g =
+ let
+ (* returns None if no new sure bound could be calculated, otherwise the new sure bound is returned *)
+ fun calc_sure_bound_from_sources g (key as (_, btype)) =
+ let
+ fun mult_upper x (lower, upper) =
+ if FloatArith.is_negative x then
+ FloatArith.mul x lower
+ else
+ FloatArith.mul x upper
+
+ fun mult_lower x (lower, upper) =
+ if FloatArith.is_negative x then
+ FloatArith.mul x upper
+ else
+ FloatArith.mul x lower
+
+ val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower
+
+ fun calc_sure_bound (sure_bound, (row_index, (row_bound, sources))) =
+ let
+ fun add_src_bound (sum, (coeff, src_key)) =
+ case sum of
+ None => None
+ | Some x =>
+ (case get_sure_bound g src_key of
+ None => None
+ | Some src_sure_bound => Some (FloatArith.add x (mult_btype src_sure_bound coeff)))
+ in
+ case foldl add_src_bound (Some row_bound, sources) of
+ None => sure_bound
+ | new_sure_bound as (Some new_bound) =>
+ (case sure_bound of
+ None => new_sure_bound
+ | Some old_bound =>
+ Some (case btype of
+ UPPER => FloatArith.min old_bound new_bound
+ | LOWER => FloatArith.max old_bound new_bound))
+ end
+ in
+ case VarGraph.lookup (g, key) of
+ None => None
+ | Some (sure_bound, f) =>
+ let
+ val x = Inttab.foldl calc_sure_bound (sure_bound, f)
+ in
+ if x = sure_bound then None else x
+ end
+ end
+
+ fun propagate ((g, b), (key, _)) =
+ case calc_sure_bound_from_sources g key of
+ None => (g,b)
+ | Some bound => (update_sure_bound g key bound,
+ if safe then
+ case get_sure_bound g key of
+ None => true
+ | _ => b
+ else
+ true)
+
+ val (g, b) = VarGraph.foldl propagate ((g, false), g)
+ in
+ if b then propagate_sure_bounds safe names g else g
+ end
+
+exception Load of string;
+
+fun calcr safe_propagation xlen names prec A b =
+ let
+ val empty = Inttab.empty
+
+ fun instab t i x = Inttab.update ((i,x), t)
+
+ fun isnegstr x = String.isPrefix "-" x
+ fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else "-"^x
+
+ fun test_1 (lower, upper) =
+ if lower = upper then
+ (if FloatArith.is_equal lower (IntInf.fromInt ~1, FloatArith.izero) then ~1
+ else if FloatArith.is_equal lower (IntInf.fromInt 1, FloatArith.izero) then 1
+ else 0)
+ else 0
+
+ fun calcr (g, (row_index, a)) =
+ let
+ val b = FloatSparseMatrixBuilder.v_elem_at b row_index
+ val (_, b2) = ExactFloatingPoint.approx_decstr_by_bin prec (case b of None => "0" | Some b => b)
+ val approx_a = FloatSparseMatrixBuilder.v_fold (fn (l, (i,s)) =>
+ (i, ExactFloatingPoint.approx_decstr_by_bin prec s)::l) [] a
+
+ fun fold_dest_nodes (g, (dest_index, dest_value)) =
+ let
+ val dest_test = test_1 dest_value
+ in
+ if dest_test = 0 then
+ g
+ else let
+ val (dest_key as (_, dest_btype), row_bound) =
+ if dest_test = ~1 then
+ ((dest_index, LOWER), FloatArith.neg b2)
+ else
+ ((dest_index, UPPER), b2)
+
+ fun fold_src_nodes (g, (src_index, src_value as (src_lower, src_upper))) =
+ if src_index = dest_index then g
+ else
+ let
+ val coeff = case dest_btype of
+ UPPER => (FloatArith.neg src_upper, FloatArith.neg src_lower)
+ | LOWER => src_value
+ in
+ if FloatArith.is_negative src_lower then
+ add_edge g (src_index, UPPER) dest_key row_index coeff
+ else
+ add_edge g (src_index, LOWER) dest_key row_index coeff
+ end
+ in
+ foldl fold_src_nodes ((add_row_bound g dest_key row_index row_bound), approx_a)
+ end
+ end
+ in
+ case approx_a of
+ [] => g
+ | [(u, a)] =>
+ let
+ val atest = test_1 a
+ in
+ if atest = ~1 then
+ update_sure_bound g (u, LOWER) (FloatArith.neg b2)
+ else if atest = 1 then
+ update_sure_bound g (u, UPPER) b2
+ else
+ g
+ end
+ | _ => foldl fold_dest_nodes (g, approx_a)
+ end
+
+ val g = FloatSparseMatrixBuilder.m_fold calcr VarGraph.empty A
+ val g = propagate_sure_bounds safe_propagation names g
+
+ val (r1, r2) = split_graph g
+
+ fun abs_estimate i r1 r2 =
+ if i = 0 then FloatSparseMatrixBuilder.empty_spmat
+ else
+ let
+ val index = xlen-i
+ val r = abs_estimate (i-1) r1 r2
+ val b1 = case Inttab.lookup (r1, index) of None => raise (Load ("x-value not bounded from below: "^(names index))) | Some x => x
+ val b2 = case Inttab.lookup (r2, index) of None => raise (Load ("x-value not bounded from above: "^(names index))) | Some x => x
+ val abs_max = FloatArith.max (FloatArith.neg (FloatArith.negative_part b1)) (FloatArith.positive_part b2)
+ val vec = FloatSparseMatrixBuilder.cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 abs_max) FloatSparseMatrixBuilder.empty_spvec
+ in
+ FloatSparseMatrixBuilder.cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) r
+ end
+ in
+ FloatSparseMatrixBuilder.sign_term (abs_estimate xlen r1 r2)
+ end
+
+fun load filename prec safe_propagation =
+ let
+ val prog = Cplex.load_cplexFile filename
+ val prog = Cplex.elim_nonfree_bounds prog
+ val prog = Cplex.relax_strict_ineqs prog
+ val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog
+ val r = calcr safe_propagation xlen names prec A b
+ val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems"
+ val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b
+ val results = Cplex.solve dualprog
+ val (optimal,v) = CplexFloatSparseMatrixConverter.convert_results results indexof
+ val A = FloatSparseMatrixBuilder.cut_matrix v None A
+ fun id x = x
+ val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v
+ val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b)
+ val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c
+ val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec FloatArith.positive_part v
+ val A = FloatSparseMatrixBuilder.approx_matrix prec id A
+ val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b
+ val c = FloatSparseMatrixBuilder.approx_matrix prec id c
+ in
+ (y1, A, b2, c, r)
+ end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s)))
+
+end
\ No newline at end of file
--- a/src/HOL/OrderedGroup.thy Fri Sep 03 10:27:05 2004 +0200
+++ b/src/HOL/OrderedGroup.thy Fri Sep 03 17:10:36 2004 +0200
@@ -842,6 +842,42 @@
lemma minus_add_cancel: "-(a::'a::ab_group_add) + (a + b) = b"
by (simp add: add_assoc[symmetric])
+lemma le_add_right_mono:
+ assumes
+ "a <= b + (c::'a::pordered_ab_group_add)"
+ "c <= d"
+ shows "a <= b + d"
+ apply (rule_tac order_trans[where y = "b+c"])
+ apply (simp_all add: prems)
+ done
+
+lemmas group_eq_simps =
+ mult_ac
+ add_ac
+ add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2
+ diff_eq_eq eq_diff_eq
+
+lemma estimate_by_abs:
+"a + b <= (c::'a::lordered_ab_group_abs) \<Longrightarrow> a <= c + abs b"
+proof -
+ assume 1: "a+b <= c"
+ have 2: "a <= c+(-b)"
+ apply (insert 1)
+ apply (drule_tac add_right_mono[where c="-b"])
+ apply (simp add: group_eq_simps)
+ done
+ have 3: "(-b) <= abs b" by (rule abs_ge_minus_self)
+ show ?thesis by (rule le_add_right_mono[OF 2 3])
+qed
+
+lemma abs_of_ge_0: "0 <= (y::'a::lordered_ab_group_abs) \<Longrightarrow> abs y = y"
+proof -
+ assume 1:"0 <= y"
+ have 2:"-y <= 0" by (simp add: 1)
+ from 1 2 have 3:"-y <= y" by (simp only:)
+ show ?thesis by (simp add: abs_lattice ge_imp_join_eq[OF 3])
+qed
+
ML {*
val add_zero_left = thm"add_0";
val add_zero_right = thm"add_0_right";
--- a/src/HOL/Ring_and_Field.thy Fri Sep 03 10:27:05 2004 +0200
+++ b/src/HOL/Ring_and_Field.thy Fri Sep 03 17:10:36 2004 +0200
@@ -546,14 +546,14 @@
done
text{*This list of rewrites decides ring equalities by ordered rewriting.*}
-lemmas ring_eq_simps =
- mult_ac
+lemmas ring_eq_simps =
+(* mult_ac*)
left_distrib right_distrib left_diff_distrib right_diff_distrib
- add_ac
+ group_eq_simps
+(* add_ac
add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2
- diff_eq_eq eq_diff_eq
+ diff_eq_eq eq_diff_eq *)
-
subsection {* Fields *}
lemma right_inverse [simp]:
@@ -1503,6 +1503,90 @@
text{*Moving this up spoils many proofs using @{text mult_le_cancel_right}*}
declare times_divide_eq_left [simp]
+lemma linprog_dual_estimate:
+ assumes
+ "A * x \<le> (b::'a::lordered_ring)"
+ "0 \<le> y"
+ "abs (A - A') \<le> \<delta>A"
+ "b \<le> b'"
+ "abs (c - c') \<le> \<delta>c"
+ "abs x \<le> r"
+ shows
+ "c * x \<le> y * b' + (y * \<delta>A + abs (y * A' - c') + \<delta>c) * r"
+proof -
+ from prems have 1: "y * b <= y * b'" by (simp add: mult_left_mono)
+ from prems have 2: "y * (A * x) <= y * b" by (simp add: mult_left_mono)
+ have 3: "y * (A * x) = c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x" by (simp add: ring_eq_simps)
+ from 1 2 3 have 4: "c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x <= y * b'" by simp
+ have 5: "c * x <= y * b' + abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)"
+ by (simp only: 4 estimate_by_abs)
+ have 6: "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x) <= abs (y * (A - A') + (y * A' - c') + (c'-c)) * abs x"
+ by (simp add: abs_le_mult)
+ 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"
+ by (simp add: abs_triangle_ineq mult_right_mono)
+ 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"
+ by (simp add: abs_triangle_ineq mult_right_mono)
+ 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"
+ by (simp add: abs_le_mult mult_right_mono)
+ have 10: "c'-c = -(c-c')" by (simp add: ring_eq_simps)
+ have 11: "abs (c'-c) = abs (c-c')"
+ by (subst 10, subst abs_minus_cancel, simp)
+ 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"
+ by (simp add: 11 prems mult_right_mono)
+ 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"
+ by (simp add: prems mult_right_mono mult_left_mono)
+ 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"
+ apply (rule mult_left_mono)
+ apply (simp add: prems)
+ apply (rule_tac add_mono[of "0::'a" _ "0", simplified])+
+ apply (rule mult_left_mono[of "0" "\<delta>A", simplified])
+ apply (simp_all)
+ apply (rule order_trans[where y="abs (A-A')"], simp_all add: prems)
+ apply (rule order_trans[where y="abs (c-c')"], simp_all add: prems)
+ done
+ 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"
+ by (simp)
+ show ?thesis
+ apply (rule_tac le_add_right_mono[of _ _ "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)"])
+ apply (simp_all add: 5 14[simplified abs_of_ge_0[of y, simplified prems]])
+ done
+qed
+
+lemma le_ge_imp_abs_diff_1:
+ assumes
+ "A1 <= (A::'a::lordered_ring)"
+ "A <= A2"
+ shows "abs (A-A1) <= A2-A1"
+proof -
+ have "0 <= A - A1"
+ proof -
+ have 1: "A - A1 = A + (- A1)" by simp
+ show ?thesis by (simp only: 1 add_right_mono[of A1 A "-A1", simplified, simplified prems])
+ qed
+ then have "abs (A-A1) = A-A1" by (rule abs_of_ge_0)
+ with prems show "abs (A-A1) <= (A2-A1)" by simp
+qed
+
+lemma linprog_dual_estimate_1:
+ assumes
+ "A * x \<le> (b::'a::lordered_ring)"
+ "0 \<le> y"
+ "A1 <= A"
+ "A <= A2"
+ "c1 <= c"
+ "c <= c2"
+ "abs x \<le> r"
+ shows
+ "c * x \<le> y * b + (y * (A2 - A1) + abs (y * A1 - c1) + (c2 - c1)) * r"
+proof -
+ from prems have delta_A: "abs (A-A1) <= (A2-A1)" by (simp add: le_ge_imp_abs_diff_1)
+ from prems have delta_c: "abs (c-c1) <= (c2-c1)" by (simp add: le_ge_imp_abs_diff_1)
+ show ?thesis
+ apply (rule_tac linprog_dual_estimate)
+ apply (auto intro: delta_A delta_c simp add: prems)
+ done
+qed
+
ML {*
val left_distrib = thm "left_distrib";
val right_distrib = thm "right_distrib";