summaryrefslogtreecommitdiff
path: root/contrib/micromega/vector.ml
diff options
context:
space:
mode:
Diffstat (limited to 'contrib/micromega/vector.ml')
-rw-r--r--contrib/micromega/vector.ml674
1 files changed, 674 insertions, 0 deletions
diff --git a/contrib/micromega/vector.ml b/contrib/micromega/vector.ml
new file mode 100644
index 00000000..fee4ebfc
--- /dev/null
+++ b/contrib/micromega/vector.ml
@@ -0,0 +1,674 @@
+(************************************************************************)
+(* v * The Coq Proof Assistant / The Coq Development Team *)
+(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *)
+(* \VV/ **************************************************************)
+(* // * This file is distributed under the terms of the *)
+(* * GNU Lesser General Public License Version 2.1 *)
+(************************************************************************)
+(* *)
+(* Micromega: A reflexive tactic using the Positivstellensatz *)
+(* *)
+(* Frédéric Besson (Irisa/Inria) 2006-2008 *)
+(* *)
+(************************************************************************)
+
+open Num
+
+module type S =
+sig
+ type t
+
+ val fresh : t -> int
+
+ val null : t
+
+ val is_null : t -> bool
+
+ val get : int -> t -> num
+
+ val update : int -> (num -> num) -> t -> t
+ (* behaviour is undef if index < 0 -- might loop*)
+
+ val set : int -> num -> t -> t
+
+ (*
+ For efficiency...
+
+ val get_update : int -> (num -> num) -> t -> num * t
+ *)
+
+ val mul : num -> t -> t
+
+ val uminus : t -> t
+
+ val add : t -> t -> t
+
+ val dotp : t -> t -> num
+
+ val lin_comb : num -> t -> num -> t -> t
+ (* lin_comb n1 t1 n2 t2 = (n1 * t1) + (n2 * t2) *)
+
+ val gcd : t -> Big_int.big_int
+
+ val normalise : t -> num * t
+
+ val hash : t -> int
+
+ val compare : t -> t -> int
+
+ type it
+
+ val iterator : t -> it
+ val element : it -> (num*it) option
+
+ val string : t -> string
+
+ type status = Pos | Neg
+
+ (* the result list is ordered by fst *)
+ val status : t -> (int * status) list
+
+ val from_list : num list -> t
+ val to_list : t -> num list
+
+end
+
+
+module type SystemS =
+sig
+
+ module Vect : S
+
+ module Cstr :
+ sig
+ type kind = Eq | Ge
+ val string_of_kind : kind -> string
+ type cstr = {coeffs : Vect.t ; op : kind ; cst : num}
+ val string_of_cstr : cstr -> string
+ val compare : cstr -> cstr -> int
+ end
+ open Cstr
+
+
+ module CstrBag :
+ sig
+ type t
+ exception Contradiction
+
+ val empty : t
+
+ val is_empty : t -> bool
+
+ val add : cstr -> t -> t
+ (* c can be deduced from add c t *)
+
+ val find : (cstr -> bool) -> t -> cstr option
+
+ val fold : (cstr -> 'a -> 'a) -> t -> 'a -> 'a
+
+ val status : t -> (int * (int list * int list)) list
+ (* aggregate of vector statuses *)
+
+ val remove : cstr -> t -> t
+
+ (* remove_list the ith element -- it is the ith element visited by 'fold' *)
+
+ val split : (cstr -> int) -> t -> (int -> t)
+
+ type it
+ val iterator : t -> it
+ val element : it -> (cstr*it) option
+
+ end
+
+end
+
+let zero_num = Int 0
+let unit_num = Int 1
+
+
+
+
+module Cstr(V:S) =
+struct
+ type kind = Eq | Ge
+ let string_of_kind = function Eq -> "Eq" | Ge -> "Ge"
+
+ type cstr = {coeffs : V.t ; op : kind ; cst : num}
+
+ let string_of_cstr {coeffs =a ; op = b ; cst =c} =
+ Printf.sprintf "{coeffs = %s;op=%s;cst=%s}" (V.string a) (string_of_kind b) (string_of_num c)
+
+ type t = cstr
+ let compare
+ {coeffs = v1 ; op = op1 ; cst = c1}
+ {coeffs = v2 ; op = op2 ; cst = c2} =
+ Mutils.Cmp.compare_lexical [
+ (fun () -> V.compare v1 v2);
+ (fun () -> Pervasives.compare op1 op2);
+ (fun () -> compare_num c1 c2)
+ ]
+
+
+end
+
+
+
+module VList : S with type t = num list =
+struct
+ type t = num list
+
+ let fresh l = failwith "not implemented"
+
+ let null = []
+
+ let is_null = List.for_all ((=/) zero_num)
+
+ let normalise l = failwith "Not implemented"
+ (*match l with (* Buggy : What if the first num is zero! *)
+ | [] -> (Int 0,[])
+ | [n] -> (n,[Int 1])
+ | n::l -> (n, (Int 1) :: List.map (fun x -> x // n) l)
+ *)
+
+
+ let get i l = try List.nth l i with _ -> zero_num
+
+ (* This is not tail-recursive *)
+ let rec update i f t =
+ match t with
+ | [] -> if i = 0 then [f zero_num] else (zero_num)::(update (i-1) f [])
+ | e::t -> if i = 0 then (f e)::t else e::(update (i-1) f t)
+
+ let rec set i n t =
+ match t with
+ | [] -> if i = 0 then [n] else (zero_num)::(set (i-1) n [])
+ | e::t -> if i = 0 then (n)::t else e::(set (i-1) n t)
+
+
+
+
+ let rec mul z t =
+ match z with
+ | Int 0 -> null
+ | Int 1 -> t
+ | _ -> List.map (mult_num z) t
+
+ let uminus t = mul (Int (-1)) t
+
+ let rec add t1 t2 =
+ match t1,t2 with
+ | [], _ -> t2
+ | _ , [] -> t1
+ | e1::t1,e2::t2 -> (e1 +/ e2 )::(add t1 t2)
+
+ let dotp t1 t2 =
+ let rec _dotp t1 t2 acc =
+ match t1, t2 with
+ | [] , _ -> acc
+ | _ , [] -> acc
+ | e1::t1,e2::t2 -> _dotp t1 t2 (acc +/ (e1 */ e2)) in
+ _dotp t1 t2 zero_num
+
+ let add_mul n t1 t2 =
+ match n with
+ | Int 0 -> t2
+ | Int 1 -> add t1 t2
+ | _ ->
+ let rec _add_mul t1 t2 =
+ match t1,t2 with
+ | [], _ -> t2
+ | _ , [] -> mul n t1
+ | e1::t1,e2::t2 -> ( (n */e1) +/ e2 )::(_add_mul t1 t2) in
+ _add_mul t1 t2
+
+ let lin_comb n1 t1 n2 t2 =
+ match n1,n2 with
+ | Int 0 , _ -> mul n2 t2
+ | Int 1 , _ -> add_mul n2 t2 t1
+ | _ , Int 0 -> mul n1 t1
+ | _ , Int 1 -> add_mul n1 t1 t2
+ | _ ->
+ let rec _lin_comb t1 t2 =
+ match t1,t2 with
+ | [], _ -> mul n2 t2
+ | _ , [] -> mul n1 t1
+ | e1::t1,e2::t2 -> ( (n1 */e1) +/ (n2 */ e2 ))::(_lin_comb t1 t2) in
+ _lin_comb t1 t2
+
+ (* could be computed on the fly *)
+ let gcd t =Mutils.gcd_list t
+
+
+
+
+ let hash = Mutils.Cmp.hash_list int_of_num
+
+ let compare = Mutils.Cmp.compare_list compare_num
+
+ type it = t
+ let iterator (x:t) : it = x
+ let element it =
+ match it with
+ | [] -> None
+ | e::l -> Some (e,l)
+
+ (* TODO: Buffer! *)
+ let string l = List.fold_right (fun n s -> (string_of_num n)^";"^s) l ""
+
+ type status = Pos | Neg
+
+ let status l =
+ let rec xstatus i l =
+ match l with
+ | [] -> []
+ | e::l ->
+ begin
+ match compare_num e (Int 0) with
+ | 1 -> (i,Pos):: (xstatus (i+1) l)
+ | 0 -> xstatus (i+1) l
+ | -1 -> (i,Neg) :: (xstatus (i+1) l)
+ | _ -> assert false
+ end in
+ xstatus 0 l
+
+ let from_list l = l
+ let to_list l = l
+
+end
+
+module VMap : S =
+struct
+ module Map = Map.Make(struct type t = int let compare (x:int) (y:int) = Pervasives.compare x y end)
+
+ type t = num Map.t
+
+ let null = Map.empty
+
+ let fresh m = failwith "not implemented"
+
+ let is_null = Map.is_empty
+
+ let normalise m = failwith "Not implemented"
+
+
+
+ let get i l = try Map.find i l with _ -> zero_num
+
+ let update i f t =
+ try
+ let res = f (Map.find i t) in
+ if res =/ zero_num
+ then Map.remove i t
+ else Map.add i res t
+ with
+ Not_found ->
+ let res = f zero_num in
+ if res =/ zero_num then t else Map.add i res t
+
+ let set i n t =
+ if n =/ zero_num then Map.remove i t
+ else Map.add i n t
+
+
+ let rec mul z t =
+ match z with
+ | Int 0 -> null
+ | Int 1 -> t
+ | _ -> Map.map (mult_num z) t
+
+ let uminus t = mul (Int (-1)) t
+
+
+ let map2 f m1 m2 =
+ let res,m2' =
+ Map.fold (fun k e (res,m2) ->
+ let v = f e (get k m2) in
+ if v =/ zero_num
+ then (res,Map.remove k m2)
+ else (Map.add k v res,Map.remove k m2)) m1 (Map.empty,m2) in
+ Map.fold (fun k e res ->
+ let v = f zero_num e in
+ if v =/ zero_num
+ then res else Map.add k v res) m2' res
+
+ let add t1 t2 = map2 (+/) t1 t2
+
+
+ let dotp t1 t2 =
+ Map.fold (fun k e res ->
+ res +/ (e */ get k t2)) t1 zero_num
+
+
+
+ let add_mul n t1 t2 =
+ match n with
+ | Int 0 -> t2
+ | Int 1 -> add t1 t2
+ | _ -> map2 (fun x y -> (n */ x) +/ y) t1 t2
+
+ let lin_comb n1 t1 n2 t2 =
+ match n1,n2 with
+ | Int 0 , _ -> mul n2 t2
+ | Int 1 , _ -> add_mul n2 t2 t1
+ | _ , Int 0 -> mul n1 t1
+ | _ , Int 1 -> add_mul n1 t1 t2
+ | _ -> map2 (fun x y -> (n1 */ x) +/ (n2 */ y)) t1 t2
+
+
+ let hash map = Map.fold (fun k e res -> k lxor (int_of_num e) lxor res) map 0
+
+ let compare = Map.compare compare_num
+
+ type it = t * int
+
+ let iterator (x:t) : it = (x,0)
+
+ let element (mp,id) =
+ try
+ Some (Map.find id mp, (mp, id+1))
+ with
+ Not_found -> None
+
+ (* TODO: Buffer! *)
+ type status = Pos | Neg
+
+ let status l = Map.fold (fun k e l ->
+ match compare_num e (Int 0) with
+ | 1 -> (k,Pos)::l
+ | 0 -> l
+ | -1 -> (k,Neg) :: l
+ | _ -> assert false) l []
+ let from_list l =
+ let rec from_list i l map =
+ match l with
+ | [] -> map
+ | e::l -> from_list (i+1) l (if e <>/ Int 0 then Map.add i e map else map) in
+ from_list 0 l Map.empty
+
+ let gcd m =
+ let res = Map.fold (fun _ e x -> Big_int.gcd_big_int x (Mutils.numerator e)) m Big_int.zero_big_int in
+ if Big_int.compare_big_int res Big_int.zero_big_int = 0
+ then Big_int.unit_big_int else res
+
+
+ let to_list m =
+ let l = List.rev (Map.fold (fun k e l -> (k,e)::l) m []) in
+ let rec xto_list i l =
+ match l with
+ | [] -> []
+ | (x,v)::l' -> if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in
+ xto_list 0 l
+
+ let string l = VList.string (to_list l)
+
+
+end
+
+
+module VSparse : S =
+struct
+
+ type t = (int*num) list
+
+ let null = []
+
+ let fresh l = List.fold_left (fun acc (i,_) -> max (i+1) acc) 0 l
+
+ let is_null l = l = []
+
+ let rec is_sorted l =
+ match l with
+ | [] -> true
+ | [e] -> true
+ | (i,_)::(j,x)::l -> i < j && is_sorted ((j,x)::l)
+
+
+ let check l = (List.for_all (fun (_,n) -> compare_num n (Int 0) <> 0) l) && (is_sorted l)
+
+ (* let get i t =
+ assert (check t);
+ try List.assoc i t with Not_found -> zero_num *)
+
+ let rec get (i:int) t =
+ match t with
+ | [] -> zero_num
+ | (j,n)::t ->
+ match compare i j with
+ | 0 -> n
+ | 1 -> get i t
+ | _ -> zero_num
+
+ let cons i v rst = if v =/ Int 0 then rst else (i,v)::rst
+
+ let rec update i f t =
+ match t with
+ | [] -> cons i (f zero_num) []
+ | (k,v)::l ->
+ match Pervasives.compare i k with
+ | 0 -> cons k (f v) l
+ | -1 -> cons i (f zero_num) t
+ | 1 -> (k,v) ::(update i f l)
+ | _ -> failwith "compare_num"
+
+ let update i f t =
+ assert (check t);
+ let res = update i f t in
+ assert (check t) ; res
+
+
+ let rec set i n t =
+ match t with
+ | [] -> cons i n []
+ | (k,v)::l ->
+ match Pervasives.compare i k with
+ | 0 -> cons k n l
+ | -1 -> cons i n t
+ | 1 -> (k,v) :: (set i n l)
+ | _ -> failwith "compare_num"
+
+
+ let rec map f l =
+ match l with
+ | [] -> []
+ | (i,e)::l -> cons i (f e) (map f l)
+
+ let rec mul z t =
+ match z with
+ | Int 0 -> null
+ | Int 1 -> t
+ | _ -> List.map (fun (i,n) -> (i, mult_num z n)) t
+
+ let mul z t =
+ assert (check t) ;
+ let res = mul z t in
+ assert (check res) ;
+ res
+
+ let uminus t = mul (Int (-1)) t
+
+
+ let normalise l =
+ match l with
+ | [] -> (Int 0,[])
+ | (i,n)::_ -> (n, mul ((Int 1) // n) l)
+
+
+ let rec map2 f m1 m2 =
+ match m1, m2 with
+ | [] , [] -> []
+ | l , [] -> map (fun x -> f x zero_num) l
+ | [] ,l -> map (f zero_num) l
+ | (i,e)::l1,(i',e')::l2 ->
+ match Pervasives.compare i i' with
+ | 0 -> cons i (f e e') (map2 f l1 l2)
+ | -1 -> cons i (f e zero_num) (map2 f l1 m2)
+ | 1 -> cons i' (f zero_num e') (map2 f m1 l2)
+ | _ -> assert false
+
+ (* let add t1 t2 = map2 (+/) t1 t2*)
+
+ let rec add (m1:t) (m2:t) =
+ match m1, m2 with
+ | [] , [] -> []
+ | l , [] -> l
+ | [] ,l -> l
+ | (i,e)::l1,(i',e')::l2 ->
+ match Pervasives.compare i i' with
+ | 0 -> cons i ( e +/ e') (add l1 l2)
+ | -1 -> (i,e) :: (add l1 m2)
+ | 1 -> (i', e') :: (add m1 l2)
+ | _ -> assert false
+
+
+
+
+ let add t1 t2 =
+ assert (check t1 && check t2);
+ let res = add t1 t2 in
+ assert (check res);
+ res
+
+
+ let rec dotp (t1:t) (t2:t) =
+ match t1, t2 with
+ | [] , _ -> zero_num
+ | _ , [] -> zero_num
+ | (i,e)::l1 , (i',e')::l2 ->
+ match Pervasives.compare i i' with
+ | 0 -> (e */ e') +/ (dotp l1 l2)
+ | -1 -> dotp l1 t2
+ | 1 -> dotp t1 l2
+ | _ -> assert false
+
+ let dotp t1 t2 =
+ assert (check t1 && check t2) ; dotp t1 t2
+
+ let add_mul n t1 t2 =
+ match n with
+ | Int 0 -> t2
+ | Int 1 -> add t1 t2
+ | _ -> map2 (fun x y -> (n */ x) +/ y) t1 t2
+
+ let add_mul n (t1:t) (t2:t) =
+ match n with
+ | Int 0 -> t2
+ | Int 1 -> add t1 t2
+ | _ ->
+ let rec xadd_mul m1 m2 =
+ match m1, m2 with
+ | [] , [] -> []
+ | _ , [] -> mul n m1
+ | [] , _ -> m2
+ | (i,e)::l1,(i',e')::l2 ->
+ match Pervasives.compare i i' with
+ | 0 -> cons i ( n */ e +/ e') (xadd_mul l1 l2)
+ | -1 -> (i,n */ e) :: (xadd_mul l1 m2)
+ | 1 -> (i', e') :: (xadd_mul m1 l2)
+ | _ -> assert false in
+ xadd_mul t1 t2
+
+
+
+
+ let lin_comb n1 t1 n2 t2 =
+ match n1,n2 with
+ | Int 0 , _ -> mul n2 t2
+ | Int 1 , _ -> add_mul n2 t2 t1
+ | _ , Int 0 -> mul n1 t1
+ | _ , Int 1 -> add_mul n1 t1 t2
+ | _ -> (*map2 (fun x y -> (n1 */ x) +/ (n2 */ y)) t1 t2*)
+ let rec xlin_comb m1 m2 =
+ match m1, m2 with
+ | [] , [] -> []
+ | _ , [] -> mul n1 m1
+ | [] , _ -> mul n2 m2
+ | (i,e)::l1,(i',e')::l2 ->
+ match Pervasives.compare i i' with
+ | 0 -> cons i ( n1 */ e +/ n2 */ e') (xlin_comb l1 l2)
+ | -1 -> (i,n1 */ e) :: (xlin_comb l1 m2)
+ | 1 -> (i', n2 */ e') :: (xlin_comb m1 l2)
+ | _ -> assert false in
+ xlin_comb t1 t2
+
+
+
+
+
+ let lin_comb n1 t1 n2 t2 =
+ assert (check t1 && check t2);
+ let res = lin_comb n1 t1 n2 t2 in
+ assert (check res); res
+
+ let hash = Mutils.Cmp.hash_list (fun (x,y) -> (Hashtbl.hash x) lxor (int_of_num y))
+
+
+ let compare = Mutils.Cmp.compare_list (fun x y -> Mutils.Cmp.compare_lexical
+ [
+ (fun () -> Pervasives.compare (fst x) (fst y));
+ (fun () -> compare_num (snd x) (snd y))])
+
+ (*
+ let compare (x:t) (y:t) =
+ let rec xcompare acc1 acc2 x y =
+ match x , y with
+ | [] , [] -> xcomp acc1 acc2
+ | [] , _ -> -1
+ | _ , [] -> 1
+ | (i,n1)::l1 , (j,n2)::l2 ->
+ match Pervasives.compare i j with
+ | 0 -> xcompare (n1::acc1) (n2::acc2) l1 l2
+ | c -> c
+ and xcomp acc1 acc2 = Mutils.Cmp.compare_list compare_num acc1 acc2 in
+ xcompare [] [] x y
+ *)
+
+ type it = t
+
+ let iterator (x:t) : it = x
+
+ let element l = failwith "Not_implemented"
+
+ (* TODO: Buffer! *)
+ type status = Pos | Neg
+
+ let status l = List.map (fun (i,e) ->
+ match compare_num e (Int 0) with
+ | 1 -> i,Pos
+ | -1 -> i,Neg
+ | _ -> assert false) l
+
+ let from_list (l: num list) =
+ let rec xfrom_list i l =
+ match l with
+ | [] -> []
+ | e::l ->
+ if e <>/ Int 0
+ then (i,e)::(xfrom_list (i+1) l)
+ else xfrom_list (i+1) l in
+
+ let res = xfrom_list 0 l in
+ assert (check res) ; res
+
+
+ let gcd m =
+ let res = List.fold_left (fun x (i,e) -> Big_int.gcd_big_int x (Mutils.numerator e)) Big_int.zero_big_int m in
+ if Big_int.compare_big_int res Big_int.zero_big_int = 0
+ then Big_int.unit_big_int else res
+
+ let to_list m =
+ let rec xto_list i l =
+ match l with
+ | [] -> []
+ | (x,v)::l' ->
+ if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in
+ xto_list 0 m
+
+ let to_list l =
+ assert (check l);
+ to_list l
+
+
+ let string l = VList.string (to_list l)
+
+end
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-20 14:45:18 +0000