aboutsummaryrefslogtreecommitdiffhomepage
path: root/plugins/groebner/ideal.ml4
diff options
context:
space:
mode:
Diffstat (limited to 'plugins/groebner/ideal.ml4')
-rw-r--r--plugins/groebner/ideal.ml41335
1 files changed, 1335 insertions, 0 deletions
diff --git a/plugins/groebner/ideal.ml4 b/plugins/groebner/ideal.ml4
new file mode 100644
index 000000000..73db36d46
--- /dev/null
+++ b/plugins/groebner/ideal.ml4
@@ -0,0 +1,1335 @@
+(************************************************************************)
+(* 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 *)
+(************************************************************************)
+
+(*i camlp4deps: "lib/refutpat.cmo" i*)
+(* NB: The above camlp4 extension adds a let* syntax for refutable patterns *)
+
+(*
+ Nullstellensatz par calcul de base de Grobner
+
+ On utilise une representation creuse des polynomes:
+ un monome est un tableau d'exposants (un par variable),
+ avec son degre en tete.
+ un polynome est une liste de (coefficient,monome).
+
+ L'algorithme de Buchberger a proprement parler est tire du code caml
+ extrait du code Coq ecrit par L.Thery.
+
+ *)
+
+open Utile
+
+exception NotInIdeal
+
+module type S = sig
+
+(* Monomials *)
+type mon = int array
+
+val mult_mon : int -> mon -> mon -> mon
+val deg : mon -> int
+val compare_mon : int -> mon -> mon -> int
+val div_mon : int -> mon -> mon -> mon
+val div_mon_test : int -> mon -> mon -> bool
+val ppcm_mon : int -> mon -> mon -> mon
+
+(* Polynomials *)
+
+type deg = int
+type coef
+type poly
+val repr : poly -> (coef * mon) list
+val polconst : deg -> coef -> poly
+val zeroP : poly
+val gen : deg -> int -> poly
+
+val equal : poly -> poly -> bool
+val name_var : string list ref
+val getvar : string list -> int -> string
+val lstringP : poly list -> string
+val printP : poly -> unit
+val lprintP : poly list -> unit
+
+val div_pol_coef : poly -> coef -> poly
+val plusP : deg -> poly -> poly -> poly
+val mult_t_pol : deg -> coef -> mon -> poly -> poly
+val selectdiv : deg -> mon -> poly list -> poly
+val oppP : poly -> poly
+val emultP : coef -> poly -> poly
+val multP : deg -> poly -> poly -> poly
+val puisP : deg -> poly -> int -> poly
+val contentP : poly -> coef
+val contentPlist : poly list -> coef
+val pgcdpos : coef -> coef -> coef
+val div_pol : deg -> poly -> poly -> coef -> coef -> mon -> poly
+val reduce2 : deg -> poly -> poly list -> coef * poly
+
+val poldepcontent : coef list ref
+val coefpoldep_find : poly -> poly -> poly
+val coefpoldep_set : poly -> poly -> poly -> unit
+val initcoefpoldep : deg -> poly list -> unit
+val reduce2_trace : deg -> poly -> poly list -> poly list -> poly list * poly
+val spol : deg -> poly -> poly -> poly
+val etrangers : deg -> poly -> poly -> bool
+val div_ppcm : deg -> poly -> poly -> poly -> bool
+
+val genPcPf : deg -> poly -> poly list -> poly list -> poly list
+val genOCPf : deg -> poly list -> poly list
+
+val is_homogeneous : poly -> bool
+
+type certificate =
+ { coef : coef; power : int;
+ gb_comb : poly list list; last_comb : poly list }
+val test_dans_ideal : deg -> poly -> poly list -> poly list ->
+ poly list * poly * certificate
+val in_ideal : deg -> poly list -> poly -> poly list * poly * certificate
+
+end
+
+(***********************************************************************
+ Global options
+*)
+let lexico = ref false
+let use_hmon = ref false
+
+(***********************************************************************
+ Functor
+*)
+
+module Make (P:Polynom.S) = struct
+
+ type coef = P.t
+ let coef0 = P.of_num (Num.Int 0)
+ let coef1 = P.of_num (Num.Int 1)
+ let coefm1 = P.of_num (Num.Int (-1))
+ let string_of_coef c = "["^(P.to_string c)^"]"
+
+(***********************************************************************
+ Monomes
+ tableau d'entiers
+ le premier coefficient est le degre
+ *)
+
+type mon = int array
+type deg = int
+type poly = (coef * mon) list
+
+(* d représente le nb de variables du monome *)
+
+(* Multiplication de monomes ayant le meme nb de variables = d *)
+let mult_mon d m m' =
+ let m'' = Array.create (d+1) 0 in
+ for i=0 to d do
+ m''.(i)<- (m.(i)+m'.(i));
+ done;
+ m''
+
+(* Degré d'un monome *)
+let deg m = m.(0)
+
+
+let compare_mon d m m' =
+ if !lexico
+ then (
+ (* Comparaison de monomes avec ordre du degre lexicographique
+ = on commence par regarder la 1ere variable*)
+ let res=ref 0 in
+ let i=ref 1 in (* 1 si lexico pur 0 si degre*)
+ while (!res=0) && (!i<=d) do
+ res:= (compare m.(!i) m'.(!i));
+ i:=!i+1;
+ done;
+ !res)
+ else (
+ (* degre lexicographique inverse *)
+ match compare m.(0) m'.(0) with
+ | 0 -> (* meme degre total *)
+ let res=ref 0 in
+ let i=ref d in
+ while (!res=0) && (!i>=1) do
+ res:= - (compare m.(!i) m'.(!i));
+ i:=!i-1;
+ done;
+ !res
+ | x -> x)
+
+(* Division de monome ayant le meme nb de variables *)
+let div_mon d m m' =
+ let m'' = Array.create (d+1) 0 in
+ for i=0 to d do
+ m''.(i)<- (m.(i)-m'.(i));
+ done;
+ m''
+
+let div_pol_coef p c =
+ List.map (fun (a,m) -> (P.divP a c,m)) p
+
+(* m' divise m *)
+let div_mon_test d m m' =
+ let res=ref true in
+ let i=ref 1 in
+ while (!res) && (!i<=d) do
+ res:= (m.(!i) >= m'.(!i));
+ i:=succ !i;
+ done;
+ !res
+
+
+(* Mise en place du degré total du monome *)
+let set_deg d m =
+ m.(0)<-0;
+ for i=1 to d do
+ m.(0)<- m.(i)+m.(0);
+ done;
+ m
+
+
+(* ppcm de monomes *)
+let ppcm_mon d m m' =
+ let m'' = Array.create (d+1) 0 in
+ for i=1 to d do
+ m''.(i)<- (max m.(i) m'.(i));
+ done;
+ set_deg d m''
+
+
+
+(**********************************************************************
+
+ Polynomes
+ liste de couples (coefficient,monome) ordonnee en decroissant
+
+ *)
+
+
+let repr p = p
+
+let equal =
+ Util.list_for_all2eq
+ (fun (c1,m1) (c2,m2) -> P.equal c1 c2 && m1=m2)
+
+let hash p =
+ let c = List.map fst p in
+ let m = List.map snd p in
+ List.fold_left (fun h p -> h * 17 + P.hash p) (Hashtbl.hash m) c
+
+module Hashpol = Hashtbl.Make(
+ struct
+ type t = poly
+ let equal = equal
+ let hash = hash
+ end)
+
+
+(*
+ A pretty printer for polynomials, with Maple-like syntax.
+ *)
+
+open Format
+
+let getvar lv i =
+ try (List.nth lv i)
+ with _ -> (List.fold_left (fun r x -> r^" "^x) "lv= " lv)
+ ^" i="^(string_of_int i)
+
+let string_of_pol zeroP hdP tlP coefterm monterm string_of_coef
+ dimmon string_of_exp lvar p =
+
+ let rec string_of_mon m coefone =
+ let s=ref [] in
+ for i=1 to (dimmon m) do
+ (match (string_of_exp m i) with
+ "0" -> ()
+ | "1" -> s:= (!s) @ [(getvar !lvar (i-1))]
+ | e -> s:= (!s) @ [((getvar !lvar (i-1)) ^ "^" ^ e)]);
+ done;
+ (match !s with
+ [] -> if coefone
+ then "1"
+ else ""
+ | l -> if coefone
+ then (String.concat "*" l)
+ else ( "*" ^
+ (String.concat "*" l)))
+ and string_of_term t start = let a = coefterm t and m = monterm t in
+ match (string_of_coef a) with
+ "0" -> ""
+ | "1" ->(match start with
+ true -> string_of_mon m true
+ |false -> ( "+ "^
+ (string_of_mon m true)))
+ | "-1" ->( "-" ^" "^(string_of_mon m true))
+ | c -> if (String.get c 0)='-'
+ then ( "- "^
+ (String.sub c 1
+ ((String.length c)-1))^
+ (string_of_mon m false))
+ else (match start with
+ true -> ( c^(string_of_mon m false))
+ |false -> ( "+ "^
+ c^(string_of_mon m false)))
+ and stringP p start =
+ if (zeroP p)
+ then (if start
+ then ("0")
+ else "")
+ else ((string_of_term (hdP p) start)^
+ " "^
+ (stringP (tlP p) false))
+ in
+ (stringP p true)
+
+
+
+let print_pol zeroP hdP tlP coefterm monterm string_of_coef
+ dimmon string_of_exp lvar p =
+
+ let rec print_mon m coefone =
+ let s=ref [] in
+ for i=1 to (dimmon m) do
+ (match (string_of_exp m i) with
+ "0" -> ()
+ | "1" -> s:= (!s) @ [(getvar !lvar (i-1))]
+ | e -> s:= (!s) @ [((getvar !lvar (i-1)) ^ "^" ^ e)]);
+ done;
+ (match !s with
+ [] -> if coefone
+ then print_string "1"
+ else ()
+ | l -> if coefone
+ then print_string (String.concat "*" l)
+ else (print_string "*";
+ print_string (String.concat "*" l)))
+ and print_term t start = let a = coefterm t and m = monterm t in
+ match (string_of_coef a) with
+ "0" -> ()
+ | "1" ->(match start with
+ true -> print_mon m true
+ |false -> (print_string "+ ";
+ print_mon m true))
+ | "-1" ->(print_string "-";print_space();print_mon m true)
+ | c -> if (String.get c 0)='-'
+ then (print_string "- ";
+ print_string (String.sub c 1
+ ((String.length c)-1));
+ print_mon m false)
+ else (match start with
+ true -> (print_string c;print_mon m false)
+ |false -> (print_string "+ ";
+ print_string c;print_mon m false))
+ and printP p start =
+ if (zeroP p)
+ then (if start
+ then print_string("0")
+ else ())
+ else (print_term (hdP p) start;
+ if start then open_hovbox 0;
+ print_space();
+ print_cut();
+ printP (tlP p) false)
+ in open_hovbox 3;
+ printP p true;
+ print_flush()
+
+
+let name_var= ref []
+
+let stringP = string_of_pol
+ (fun p -> match p with [] -> true | _ -> false)
+ (fun p -> match p with (t::p) -> t |_ -> failwith "print_pol dans dansideal")
+ (fun p -> match p with (t::p) -> p |_ -> failwith "print_pol dans dansideal")
+ (fun (a,m) -> a)
+ (fun (a,m) -> m)
+ string_of_coef
+ (fun m -> (Array.length m)-1)
+ (fun m i -> (string_of_int (m.(i))))
+ name_var
+
+
+let stringPcut p =
+ if (List.length p)>20
+ then (stringP [List.hd p])^" + "^(string_of_int (List.length p))^" termes"
+ else stringP p
+
+let rec lstringP l =
+ match l with
+ [] -> ""
+ |p::l -> (stringP p)^("\n")^(lstringP l)
+
+let printP = print_pol
+ (fun p -> match p with [] -> true | _ -> false)
+ (fun p -> match p with (t::p) -> t |_ -> failwith "print_pol dans dansideal")
+ (fun p -> match p with (t::p) -> p |_ -> failwith "print_pol dans dansideal")
+ (fun (a,m) -> a)
+ (fun (a,m) -> m)
+ string_of_coef
+ (fun m -> (Array.length m)-1)
+ (fun m i -> (string_of_int (m.(i))))
+ name_var
+
+
+let rec lprintP l =
+ match l with
+ [] -> ()
+ |p::l -> printP p;print_string "\n"; lprintP l
+
+
+(* Operations
+ *)
+
+let zeroP = []
+
+(* Retourne un polynome constant à d variables *)
+let polconst d c =
+ let m = Array.create (d+1) 0 in
+ let m = set_deg d m in
+ [(c,m)]
+
+
+
+(* somme de polynomes= liste de couples (int,monomes) *)
+let plusP d p q =
+ let rec plusP p q =
+ match p with
+ [] -> q
+ |t::p' ->
+ match q with
+ [] -> p
+ |t'::q' ->
+ match compare_mon d (snd t) (snd t') with
+ 1 -> t::(plusP p' q)
+ |(-1) -> t'::(plusP p q')
+ |_ -> let c=P.plusP (fst t) (fst t') in
+ match P.equal c coef0 with
+ true -> (plusP p' q')
+ |false -> (c,(snd t))::(plusP p' q')
+ in plusP p q
+
+
+(* multiplication d'un polynome par (a,monome) *)
+let mult_t_pol d a m p =
+ let rec mult_t_pol p =
+ match p with
+ [] -> []
+ |(b,m')::p -> ((P.multP a b),(mult_mon d m m'))::(mult_t_pol p)
+ in mult_t_pol p
+
+
+(* Retourne un polynome de l dont le monome de tete divise m, [] si rien *)
+let rec selectdiv d m l =
+ match l with
+ [] -> []
+ | q::r ->
+ let m'= snd (List.hd q) in
+ if div_mon_test d m m' then q else selectdiv d m r
+
+
+(* Retourne un polynome générateur 'i' à d variables *)
+let gen d i =
+ let m = Array.create (d+1) 0 in
+ m.(i) <- 1;
+ let m = set_deg d m in
+ [(coef1,m)]
+
+
+
+(* oppose d'un polynome *)
+let oppP p =
+ let rec oppP p =
+ match p with
+ [] -> []
+ |(b,m')::p -> ((P.oppP b),m')::(oppP p)
+ in oppP p
+
+
+(* multiplication d'un polynome par un coefficient par 'a' *)
+let emultP a p =
+ let rec emultP p =
+ match p with
+ [] -> []
+ |(b,m')::p -> ((P.multP a b),m')::(emultP p)
+ in emultP p
+
+
+let multP d p q =
+ let rec aux p =
+ match p with
+ [] -> []
+ |(a,m)::p' -> plusP d (mult_t_pol d a m q) (aux p')
+ in aux p
+
+
+let puisP d p n=
+ let rec puisP n =
+ match n with
+ 0 -> [coef1, Array.create (d+1) 0]
+ | 1 -> p
+ |_ -> multP d p (puisP (n-1))
+ in puisP n
+
+let pgcdpos a b = P.pgcdP a b
+let pgcd1 p q =
+ if P.equal p coef1 || P.equal p coefm1 then p else P.pgcdP p q
+
+let rec contentP p =
+ match p with
+ |[] -> coef1
+ |[a,m] -> a
+ |(a,m)::p1 -> pgcd1 a (contentP p1)
+
+let contentPlist lp =
+ match lp with
+ |[] -> coef1
+ |p::l1 -> List.fold_left (fun r q -> pgcd1 r (contentP q)) (contentP p) l1
+
+(***********************************************************************
+ Division de polynomes
+ *)
+
+let hmon = (Hashtbl.create 1000 : (mon,poly) Hashtbl.t)
+
+let find_hmon m =
+ if !use_hmon
+ then Hashtbl.find hmon m
+ else raise Not_found
+
+let add_hmon m q =
+ if !use_hmon then Hashtbl.add hmon m q
+
+let selectdiv_cache d m l =
+ try find_hmon m
+ with Not_found ->
+ match selectdiv d m l with
+ [] -> []
+ | q -> add_hmon m q; q
+
+let div_pol d p q a b m =
+(* info ".";*)
+ plusP d (emultP a p) (mult_t_pol d b m q)
+
+(* si false on ne divise que le monome de tete *)
+let reduire_les_queues = false
+
+(* reste de la division de p par les polynomes de l
+ rend le reste et le coefficient multiplicateur *)
+
+let reduce2 d p l =
+ let rec reduce p =
+ match p with
+ [] -> (coef1,[])
+ | (a,m)::p' ->
+ let q = selectdiv_cache d m l in
+ (match q with
+ [] ->
+ if reduire_les_queues
+ then
+ let (c,r)=(reduce p') in
+ (c,((P.multP a c,m)::r))
+ else (coef1,p)
+ |(b,m')::q' ->
+ let c=(pgcdpos a b) in
+ let a'= (P.divP b c) in
+ let b'=(P.oppP (P.divP a c)) in
+ let (e,r)=reduce (div_pol d p' q' a' b'
+ (div_mon d m m')) in
+ (P.multP a' e,r)) in
+ reduce p
+
+(* trace des divisions *)
+
+(* liste des polynomes de depart *)
+let poldep = ref []
+let poldepcontent = ref []
+
+
+module HashPolPair = Hashtbl.Make
+ (struct
+ type t = poly * poly
+ let equal (p,q) (p',q') = equal p p' && equal q q'
+ let hash (p,q) =
+ let c = List.map fst p @ List.map fst q in
+ let m = List.map snd p @ List.map snd q in
+ List.fold_left (fun h p -> h * 17 + P.hash p) (Hashtbl.hash m) c
+ end)
+
+(* table des coefficients des polynomes en fonction des polynomes de depart *)
+let coefpoldep = HashPolPair.create 51
+
+(* coefficient de q dans l expression de p = sum_i c_i*q_i *)
+let coefpoldep_find p q =
+ try (HashPolPair.find coefpoldep (p,q))
+ with _ -> []
+
+let coefpoldep_set p q c =
+ HashPolPair.add coefpoldep (p,q) c
+
+let initcoefpoldep d lp =
+ poldep:=lp;
+ poldepcontent:= List.map contentP (!poldep);
+ List.iter
+ (fun p -> coefpoldep_set p p (polconst d coef1))
+ lp
+
+(* garde la trace dans coefpoldep
+ divise sans pseudodivisions *)
+
+let reduce2_trace d p l lcp =
+ (* rend (lq,r), ou r = p + sum(lq) *)
+ let rec reduce p =
+ match p with
+ [] -> ([],[])
+ |t::p' -> let (a,m)=t in
+ let q =
+ (try Hashtbl.find hmon m
+ with Not_found ->
+ let q = selectdiv d m l in
+ match q with
+ t'::q' -> (Hashtbl.add hmon m q;q)
+ |[] -> q) in
+ match q with
+ [] ->
+ if reduire_les_queues
+ then
+ let (lq,r)=(reduce p') in
+ (lq,((a,m)::r))
+ else ([],p)
+ |(b,m')::q' ->
+ let b' = P.oppP (P.divP a b) in
+ let m''= div_mon d m m' in
+ let p1=plusP d p' (mult_t_pol d b' m'' q') in
+ let (lq,r)=reduce p1 in
+ ((b',m'',q)::lq, r)
+ in let (lq,r) = reduce p in
+ (*info "reduce2_trace:\n";
+ List.iter
+ (fun (a,m,s) ->
+ let x = mult_t_pol d a m s in
+ info ((stringP x)^"\n"))
+ lq;
+ info "ok\n";*)
+ (List.map2
+ (fun c0 q ->
+ let c =
+ List.fold_left
+ (fun x (a,m,s) ->
+ if equal s q
+ then
+ plusP d x (mult_t_pol d a m (polconst d coef1))
+ else x)
+ c0
+ lq in
+ c)
+ lcp
+ !poldep,
+ r)
+
+(***********************************************************************
+ Algorithme de Janet (V.P.Gerdt Involutive algorithms...)
+*)
+
+(***********************************
+ deuxieme version, qui elimine des poly inutiles
+*)
+let homogeneous = ref false
+let pol_courant = ref []
+
+
+type pol3 =
+ {pol : poly;
+ anc : poly;
+ nmp : mon}
+
+let fst_mon q = snd (List.hd q.pol)
+let lm_anc q = snd (List.hd q.anc)
+
+let pol_to_pol3 d p =
+ {pol = p; anc = p; nmp = Array.create (d+1) 0}
+
+let is_multiplicative u s i =
+ if i=1
+ then List.for_all (fun q -> (fst_mon q).(1) <= u.(1)) s
+ else
+ List.for_all
+ (fun q ->
+ let v = fst_mon q in
+ let res = ref true in
+ let j = ref 1 in
+ while !j < i && !res do
+ res:= v.(!j) = u.(!j);
+ j:= !j + 1;
+ done;
+ if !res
+ then v.(i) <= u.(i)
+ else true)
+ s
+
+let is_multiplicative_rev d u s i =
+ if i=1
+ then
+ List.for_all
+ (fun q -> (fst_mon q).(d+1-1) <= u.(d+1-1))
+ s
+ else
+ List.for_all
+ (fun q ->
+ let v = fst_mon q in
+ let res = ref true in
+ let j = ref 1 in
+ while !j < i && !res do
+ res:= v.(d+1- !j) = u.(d+1- !j);
+ j:= !j + 1;
+ done;
+ if !res
+ then v.(d+1-i) <= u.(d+1-i)
+ else true)
+ s
+
+let monom_multiplicative d u s =
+ let m = Array.create (d+1) 0 in
+ for i=1 to d do
+ if is_multiplicative u s i
+ then m.(i)<- 1;
+ done;
+ m
+
+(* mu monome des variables multiplicative de u *)
+let janet_div_mon d u mu v =
+ let res = ref true in
+ let i = ref 1 in
+ while !i <= d && !res do
+ if mu.(!i) = 0
+ then res := u.(!i) = v.(!i)
+ else res := u.(!i) <= v.(!i);
+ i:= !i + 1;
+ done;
+ !res
+
+let find_multiplicative p mg =
+ try Hashpol.find mg p.pol
+ with Not_found -> (info "\nPROBLEME DANS LA TABLE DES VAR MULT";
+ info (stringPcut p.pol);
+ failwith "aie")
+
+(* mg hashtable de p -> monome_multiplicatif de g *)
+
+let hashtbl_reductor = ref (Hashtbl.create 51 : (mon,pol3) Hashtbl.t)
+
+(* suppose que la table a été initialisée *)
+let find_reductor d v lt mt =
+ try Hashtbl.find !hashtbl_reductor v
+ with Not_found ->
+ let r =
+ List.find
+ (fun q ->
+ let u = fst_mon q in
+ let mu = find_multiplicative q mt in
+ janet_div_mon d u mu v
+ )
+ lt in
+ Hashtbl.add !hashtbl_reductor v r;
+ r
+
+let normal_form d p g mg onlyhead =
+ let rec aux = function
+ [] -> (coef1,p)
+ | (a,v)::p' ->
+ (try
+ let q = find_reductor d v g mg in
+ let* (b,u)::q' = q.pol in
+ let c = P.pgcdP a b in
+ let a' = P.divP b c in
+ let b' = P.oppP (P.divP a c) in
+ let m = div_mon d v u in
+ (* info ".";*)
+ let (c,r) = aux (plusP d (emultP a' p') (mult_t_pol d b' m q')) in
+ (P.multP c a', r)
+ with _ -> (* TODO: catch only relevant exn *)
+ if onlyhead
+ then (coef1,p)
+ else
+ let (c,r)= (aux p') in
+ (c, plusP d [(P.multP a c,v)] r)) in
+ aux p
+
+let janet_normal_form d p g mg =
+ let (_,r) = normal_form d p g mg false in r
+
+let head_janet_normal_form d p g mg =
+ let (_,r) = normal_form d p g mg true in r
+
+let reduce_rem d r lt lq =
+ Hashtbl.clear hmon;
+ let (_,r) = reduce2 d r (List.map (fun p -> p.pol) (lt @ lq)) in
+ Hashtbl.clear hmon;
+ r
+
+let tail_normal_form d p g mg =
+ let* (a,v)::p' = p in
+ let (c,r)= (normal_form d p' g mg false) in
+ plusP d [(P.multP a c,v)] r
+
+let div_strict d m1 m2 =
+ m1 <> m2 && div_mon_test d m2 m1
+
+let criteria d p g lt =
+ mult_mon d (lm_anc p) (lm_anc g) = fst_mon p
+||
+ (let c = ppcm_mon d (lm_anc p) (lm_anc g) in
+ div_strict d c (fst_mon p))
+||
+ (List.exists
+ (fun t ->
+ let cp = ppcm_mon d (lm_anc p) (fst_mon t) in
+ let cg = ppcm_mon d (lm_anc g) (fst_mon t) in
+ let c = ppcm_mon d (lm_anc p) (lm_anc g) in
+ div_strict d cp c && div_strict d cg c)
+ lt)
+
+let head_normal_form d p lt mt =
+ let h = ref (p.pol) in
+ let res =
+ try (
+ let v = snd(List.hd !h) in
+ let g = ref (find_reductor d v lt mt) in
+ if snd(List.hd !h) <> lm_anc p && criteria d p !g lt
+ then ((* info "=";*) [])
+ else (
+ while !h <> [] && (!g).pol <> [] do
+ let* (a,v)::p' = !h in
+ let* (b,u)::q' = (!g).pol in
+ let c = P.pgcdP a b in
+ let a' = P.divP b c in
+ let b' = P.oppP (P.divP a c) in
+ let m = (div_mon d v u) in
+ h := plusP d (emultP a' p') (mult_t_pol d b' m q');
+ let v = snd(List.hd !h) in
+ g := (
+ try find_reductor d v lt mt
+ with _ -> pol_to_pol3 d []);
+ done;
+ !h)
+ )
+ with _ -> ((* info ".";*) !h)
+ in
+ (*info ("temps de head_normal_form: "
+ ^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));*)
+ res
+
+let deg_hom p =
+ match p with
+ | [] -> -1
+ | (_,m)::_ -> m.(0)
+
+let head_reduce d lq lt mt =
+ let ls = ref lq in
+ let lq = ref [] in
+ while !ls <> [] do
+ let* p::ls1 = !ls in
+ ls := ls1;
+ if !homogeneous && p.pol<>[] && deg_hom p.pol > deg_hom !pol_courant
+ then info "h"
+ else
+ match head_normal_form d p lt mt with
+ (_,v)::_ as h ->
+ if fst_mon p <> v
+ then lq := (!lq) @ [{pol = h; anc = h; nmp = Array.create (d+1) 0}]
+ else lq := (!lq) @ [p]
+ | [] ->
+ (* info "*";*)
+ if fst_mon p = lm_anc p
+ then ls := List.filter (fun q -> not (equal q.anc p.anc)) !ls
+ done;
+ (*info ("temps de head_reduce: "
+ ^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));*)
+ !lq
+
+let choose_irreductible d lf =
+ List.hd lf
+(* bien plus lent
+ (List.sort (fun p q -> compare_mon d (fst_mon p.pol) (fst_mon q.pol)) lf)
+*)
+
+
+let hashtbl_multiplicative d lf =
+ let mg = Hashpol.create 51 in
+ hashtbl_reductor := Hashtbl.create 51;
+ List.iter
+ (fun g ->
+ let (_,u) = List.hd g.pol in
+ Hashpol.add mg g.pol (monom_multiplicative d u lf))
+ lf;
+ (*info ("temps de hashtbl_multiplicative: "
+ ^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));*)
+ mg
+
+let list_diff l x =
+ List.filter (fun y -> y <> x) l
+
+let janet2 d lf p0 =
+ hashtbl_reductor := Hashtbl.create 51;
+ let t1 = Unix.gettimeofday() in
+ info ("------------- Janet2 ------------------\n");
+ pol_courant := p0;
+ let r = ref p0 in
+ let lf = List.map (pol_to_pol3 d) lf in
+ let f = choose_irreductible d lf in
+ let lt = ref [f] in
+ let mt = ref (hashtbl_multiplicative d !lt) in
+ let lq = ref (list_diff lf f) in
+ lq := head_reduce d !lq !lt !mt;
+ r := (* lent head_janet_normal_form d !r !lt !mt ; *)
+ reduce_rem d !r !lt !lq ;
+ info ("reste: "^(stringPcut !r)^"\n");
+ while !lq <> [] && !r <> [] do
+ let p = choose_irreductible d !lq in
+ lq := list_diff !lq p;
+ if p.pol = p.anc
+ then ( (* on enleve de lt les pol divisibles par p et on les met dans lq *)
+ let m = fst_mon p in
+ let lt1 = !lt in
+ List.iter
+ (fun q ->
+ let m'= fst_mon q in
+ if div_strict d m m'
+ then (
+ lq := (!lq) @ [q];
+ lt := list_diff !lt q))
+ lt1;
+ mt := hashtbl_multiplicative d !lt;
+ );
+ let h = tail_normal_form d p.pol !lt !mt in
+ if h = []
+ then info "************ reduction a zero, ne devrait pas arriver *\n"
+ else (
+ lt := (!lt) @ [{pol = h; anc = p.anc; nmp = Array.copy p.nmp}];
+ info ("nouveau polynome: "^(stringPcut h)^"\n");
+ mt := hashtbl_multiplicative d !lt;
+ r := (* lent head_janet_normal_form d !r !lt !mt ; *)
+ reduce_rem d !r !lt !lq ;
+ info ("reste: "^(stringPcut !r)^"\n");
+ if !r <> []
+ then (
+ List.iter
+ (fun q ->
+ let mq = find_multiplicative q !mt in
+ for i=1 to d do
+ if mq.(i) = 1
+ then q.nmp.(i)<- 0
+ else
+ if q.nmp.(i) = 0
+ then (
+ (* info "+";*)
+ lq := (!lq) @
+ [{pol = multP d (gen d i) q.pol;
+ anc = q.anc;
+ nmp = Array.create (d+1) 0}];
+ q.nmp.(i)<-1;
+ );
+ done;
+ )
+ !lt;
+ lq := head_reduce d !lq !lt !mt;
+ info ("["^(string_of_int (List.length !lt))^";"
+ ^(string_of_int (List.length !lq))^"]");
+ ));
+ done;
+ info ("--- Janet2 donne:\n");
+ (* List.iter (fun p -> info ("polynome: "^(stringPcut p.pol)^"\n")) !lt; *)
+ info ("reste: "^(stringPcut !r)^"\n");
+ info ("--- fin Janet2\n");
+ info ("temps: "^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));
+ List. map (fun q -> q.pol) !lt
+
+(**********************************************************************
+ version 3 *)
+
+let head_normal_form3 d p lt mt =
+ let h = ref (p.pol) in
+ let res =
+ try (
+ let v = snd(List.hd !h) in
+ let g = ref (find_reductor d v lt mt) in
+ if snd(List.hd !h) <> lm_anc p && criteria d p !g lt
+ then ((* info "=";*) [])
+ else (
+ while !h <> [] && (!g).pol <> [] do
+ let* (a,v)::p' = !h in
+ let* (b,u)::q' = (!g).pol in
+ let c = P.pgcdP a b in
+ let a' = P.divP b c in
+ let b' = P.oppP (P.divP a c) in
+ let m = div_mon d v u in
+ h := plusP d (emultP a' p') (mult_t_pol d b' m q');
+ let v = snd(List.hd !h) in
+ g := (
+ try find_reductor d v lt mt
+ with _ -> pol_to_pol3 d []);
+ done;
+ !h)
+ )
+ with _ -> ((* info ".";*) !h)
+ in
+ (*info ("temps de head_normal_form: "
+ ^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));*)
+ res
+
+
+let janet3 d lf p0 =
+ hashtbl_reductor := Hashtbl.create 51;
+ let t1 = Unix.gettimeofday() in
+ info ("------------- Janet2 ------------------\n");
+ let r = ref p0 in
+ let lf = List.map (pol_to_pol3 d) lf in
+
+ let* f::lf1 = lf in
+ let lt = ref [f] in
+ let mt = ref (hashtbl_multiplicative d !lt) in
+ let lq = ref lf1 in
+ r := reduce_rem d !r !lt !lq ; (* janet_normal_form d !r !lt !mt ;*)
+ info ("reste: "^(stringPcut !r)^"\n");
+ while !lq <> [] && !r <> [] do
+ let* p::lq1 = !lq in
+ lq := lq1;
+(*
+ if p.pol = p.anc
+ then ( (* on enleve de lt les pol divisibles par p et on les met dans lq *)
+ let m = fst_mon (p.pol) in
+ let lt1 = !lt in
+ List.iter
+ (fun q ->
+ let m'= fst_mon (q.pol) in
+ if div_strict d m m'
+ then (
+ lq := (!lq) @ [q];
+ lt := list_diff !lt q))
+ lt1;
+ mt := hashtbl_multiplicative d !lt;
+ );
+*)
+ let h = head_normal_form3 d p !lt !mt in
+ if h = []
+ then (
+ info "*";
+ if fst_mon p = lm_anc p
+ then
+ lq := List.filter (fun q -> not (equal q.anc p.anc)) !lq;
+ )
+ else (
+ let q =
+ if fst_mon p <> snd(List.hd h)
+ then {pol = h; anc = h; nmp = Array.create (d+1) 0}
+ else {pol = h; anc = p.anc; nmp = Array.copy p.nmp}
+ in
+ lt:= (!lt) @ [q];
+ info ("nouveau polynome: "^(stringPcut q.pol)^"\n");
+ mt := hashtbl_multiplicative d !lt;
+ r := reduce_rem d !r !lt !lq ; (* janet_normal_form d !r !lt !mt ;*)
+ info ("reste: "^(stringPcut !r)^"\n");
+ if !r <> []
+ then (
+ List.iter
+ (fun q ->
+ let mq = find_multiplicative q !mt in
+ for i=1 to d do
+ if mq.(i) = 1
+ then q.nmp.(i)<- 0
+ else
+ if q.nmp.(i) = 0
+ then (
+ (* info "+";*)
+ lq := (!lq) @
+ [{pol = multP d (gen d i) q.pol;
+ anc = q.anc;
+ nmp = Array.create (d+1) 0}];
+ q.nmp.(i)<-1;
+ );
+ done;
+ )
+ !lt;
+ info ("["^(string_of_int (List.length !lt))^";"
+ ^(string_of_int (List.length !lq))^"]");
+ ));
+ done;
+ info ("--- Janet3 donne:\n");
+ (* List.iter (fun p -> info ("polynome: "^(stringPcut p.pol)^"\n")) !lt; *)
+ info ("reste: "^(stringPcut !r)^"\n");
+ info ("--- fin Janet3\n");
+ info ("temps: "^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));
+ List. map (fun q -> q.pol) !lt
+
+
+(***********************************************************************
+ Completion
+ *)
+
+let sugar_flag = ref true
+
+let sugartbl = (Hashpol.create 1000 : int Hashpol.t)
+
+let compute_sugar p =
+ List.fold_left (fun s (a,m) -> max s m.(0)) 0 p
+
+let sugar p =
+ try Hashpol.find sugartbl p
+ with Not_found ->
+ let s = compute_sugar p in
+ Hashpol.add sugartbl p s;
+ s
+
+let spol d p q=
+ let m = snd (List.hd p) in
+ let m'= snd (List.hd q) in
+ let a = fst (List.hd p) in
+ let b = fst (List.hd q) in
+ let p'= List.tl p in
+ let q'= List.tl q in
+ let c=(pgcdpos a b) in
+ let m''=(ppcm_mon d m m') in
+ let m1 = div_mon d m'' m in
+ let m2 = div_mon d m'' m' in
+ let fsp p' q' =
+ plusP d
+ (mult_t_pol d (P.divP b c) m1 p')
+ (mult_t_pol d (P.oppP (P.divP a c)) m2 q') in
+ let sp = fsp p' q' in
+ coefpoldep_set sp p (fsp (polconst d coef1) []);
+ coefpoldep_set sp q (fsp [] (polconst d coef1));
+ Hashpol.add sugartbl sp
+ (max (m1.(0) + (sugar p)) (m2.(0) + (sugar q)));
+ sp
+
+let etrangers d p p'=
+ let m = snd (List.hd p) in
+ let m'= snd (List.hd p') in
+ let res=ref true in
+ let i=ref 1 in
+ while (!res) && (!i<=d) do
+ res:= (m.(!i) = 0) || (m'.(!i)=0);
+ i:=!i+1;
+ done;
+ !res
+
+
+(* teste si le monome dominant de p''
+ divise le ppcm des monomes dominants de p et p' *)
+
+let div_ppcm d p p' p'' =
+ let m = snd (List.hd p) in
+ let m'= snd (List.hd p') in
+ let m''= snd (List.hd p'') in
+ let res=ref true in
+ let i=ref 1 in
+ while (!res) && (!i<=d) do
+ res:= ((max m.(!i) m'.(!i)) >= m''.(!i));
+ i:=!i+1;
+ done;
+ !res
+
+(***********************************************************************
+ Code extrait de la preuve de L.Thery en Coq
+ *)
+type 'poly cpRes =
+ Keep of ('poly list)
+ | DontKeep of ('poly list)
+
+let addRes i = function
+ Keep h'0 -> Keep (i::h'0)
+ | DontKeep h'0 -> DontKeep (i::h'0)
+
+let rec slice d i a = function
+ [] -> if etrangers d i a then DontKeep [] else Keep []
+ | b::q1 ->
+ if div_ppcm d i a b then DontKeep (b::q1)
+ else if div_ppcm d i b a then slice d i a q1
+ else addRes b (slice d i a q1)
+
+let rec addS x l = l @[x]
+
+let addSugar x l =
+ if !sugar_flag
+ then
+ let sx = sugar x in
+ let rec insere l =
+ match l with
+ | [] -> [x]
+ | y::l1 ->
+ if sx <= (sugar y)
+ then x::l
+ else y::(insere l1)
+ in insere l
+ else addS x l
+
+(* ajoute les spolynomes de i avec la liste de polynomes aP,
+ a la liste q *)
+
+let rec genPcPf d i aP q =
+ match aP with
+ [] -> q
+ | a::l1 ->
+ (match slice d i a l1 with
+ Keep l2 -> addSugar (spol d i a) (genPcPf d i l2 q)
+ | DontKeep l2 -> genPcPf d i l2 q)
+
+let rec genOCPf d = function
+ [] -> []
+ | a::l -> genPcPf d a l (genOCPf d l)
+
+let step = ref 0
+
+let infobuch p q =
+ if !step = 0
+ then (info ("[" ^ (string_of_int (List.length p))
+ ^ "," ^ (string_of_int (List.length q))
+ ^ "]"))
+
+(***********************************************************************)
+(* dans lp les nouveaux sont en queue *)
+
+let coef_courant = ref coef1
+
+
+type certificate =
+ { coef : coef; power : int;
+ gb_comb : poly list list; last_comb : poly list }
+
+let test_dans_ideal d p lp lp0 =
+ let (c,r) = reduce2 d !pol_courant lp in
+ info ("reste: "^(stringPcut r)^"\n");
+ coef_courant:= P.multP !coef_courant c;
+ pol_courant:= r;
+ if r=[]
+ then (info "polynome reduit a 0\n";
+ let lcp = List.map (fun q -> []) !poldep in
+ let c = !coef_courant in
+ let (lcq,r) = reduce2_trace d (emultP c p) lp lcp in
+ info ("r: "^(stringP r)^"\n");
+ let res=ref (emultP c p) in
+ List.iter2
+ (fun cq q -> res:=plusP d (!res) (multP d cq q);
+ )
+ lcq !poldep;
+ info ("verif somme: "^(stringP (!res))^"\n");
+ info ("coefficient: "^(stringP (polconst d c))^"\n");
+ let rec aux lp =
+ match lp with
+ |[] -> []
+ |p::lp ->
+ (List.map
+ (fun q -> coefpoldep_find p q)
+ lp)::(aux lp)
+ in
+ let liste_polynomes_de_depart = List.rev lp0 in
+ let polynome_a_tester = p in
+ let liste_des_coefficients_intermediaires =
+ (let lci = List.rev (aux (List.rev lp)) in
+ let lci = ref lci (* (List.map List.rev lci) *) in
+ List.iter (fun x -> lci := List.tl (!lci)) lp0;
+ !lci) in
+ let liste_des_coefficients =
+ List.map
+ (fun cq -> emultP coefm1 cq)
+ (List.rev lcq) in
+ (liste_polynomes_de_depart,
+ polynome_a_tester,
+ {coef=c;
+ power=1;
+ gb_comb = liste_des_coefficients_intermediaires;
+ last_comb = liste_des_coefficients}))
+ else ((*info "polynome non reduit a 0\n";
+ info ("\nreste: "^(stringPcut r)^"\n");*)
+ raise NotInIdeal)
+
+let divide_rem_with_critical_pair = ref false
+
+let choix_spol d p l =
+ if !divide_rem_with_critical_pair
+ then (
+ let (_,m) = List.hd p in
+ try (
+ let q =
+ List.find
+ (fun q ->
+ let (_,m') = List.hd q in
+ div_mon_test d m m')
+ l in
+ q::(list_diff l q))
+ with _ -> l)
+ else l
+
+let pbuchf d pq p lp0=
+ info "calcul de la base de Groebner\n";
+ step:=0;
+ Hashtbl.clear hmon;
+ let rec pbuchf lp lpc =
+ infobuch lp lpc;
+(* step:=(!step+1)mod 10;*)
+ match lpc with
+ [] -> test_dans_ideal d p lp lp0
+ | _ ->
+ let* a::lpc2 = choix_spol d !pol_courant lpc in
+ if !homogeneous && a<>[] && deg_hom a > deg_hom !pol_courant
+ then (info "h";pbuchf lp lpc2)
+ else
+ let sa = sugar a in
+ let (ca,a0)= reduce2 d a lp in
+ match a0 with
+ [] -> info "0";pbuchf lp lpc2
+ | _ ->
+ let a1 = emultP ca a in
+ List.iter
+ (fun q ->
+ coefpoldep_set a1 q (emultP ca (coefpoldep_find a q)))
+ !poldep;
+ let (lca,a0) = reduce2_trace d a1 lp
+ (List.map (fun q -> coefpoldep_find a1 q) !poldep) in
+ List.iter2 (fun c q -> coefpoldep_set a0 q c) lca !poldep;
+ info ("\nnouveau polynome: "^(stringPcut a0)^"\n");
+ let ct = coef1 (* contentP a0 *) in
+ (*info ("contenu: "^(string_of_coef ct)^"\n");*)
+ Hashpol.add sugartbl a0 sa;
+ poldep:=addS a0 lp;
+ poldepcontent:=addS ct (!poldepcontent);
+ try test_dans_ideal d p (addS a0 lp) lp0
+ with NotInIdeal -> pbuchf (addS a0 lp) (genPcPf d a0 lp lpc2)
+ in pbuchf (fst pq) (snd pq)
+
+let is_homogeneous p =
+ match p with
+ | [] -> true
+ | (a,m)::p1 -> let d = m.(0) in
+ List.for_all (fun (b,m') -> m'.(0)=d) p1
+
+(* rend
+ c
+ lp = [pn;...;p1]
+ p
+ lci = [[a(n+1,n);...;a(n+1,1)];
+ [a(n+2,n+1);...;a(n+2,1)];
+ ...
+ [a(n+m,n+m-1);...;a(n+m,1)]]
+ lc = [qn+m; ... q1]
+
+ tels que
+ c*p = sum qi*pi
+ ou pn+k = a(n+k,n+k-1)*pn+k-1 + ... + a(n+k,1)* p1
+ *)
+
+let in_ideal d lp p =
+ homogeneous := List.for_all is_homogeneous (p::lp);
+ if !homogeneous then info "polynomes homogenes\n";
+ (* janet2 d lp p;*)
+ info ("p: "^stringPcut p^"\n");
+ info ("lp:\n"^List.fold_left (fun r p -> r^stringPcut p^"\n") "" lp);
+ initcoefpoldep d lp;
+ coef_courant:=coef1;
+ pol_courant:=p;
+ try test_dans_ideal d p lp lp
+ with NotInIdeal -> pbuchf d (lp, genOCPf d lp) p lp
+
+end
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-30 13:30:55 +0000