diff options
author | barras <barras@85f007b7-540e-0410-9357-904b9bb8a0f7> | 2009-04-16 17:27:00 +0000 |
---|---|---|
committer | barras <barras@85f007b7-540e-0410-9357-904b9bb8a0f7> | 2009-04-16 17:27:00 +0000 |
commit | 665761c4e82504a579d335b529b2e61f0414665c (patch) | |
tree | f24a021023e53a0b092a42e89c4cc65cefaf800f | |
parent | b596e8303c0d22c5df1aed3c7a56b6af863a1b9e (diff) |
nouvelle version de la tactique groebner proposee par Loic:
- algo de Janet (finalement pas utilise)
- le hash-cons des certificats de Benjamin et Laurent pas encore integre
- traitement des puissances jusqu'a 6 (methode totalement naive)
git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/coq/trunk@12088 85f007b7-540e-0410-9357-904b9bb8a0f7
-rw-r--r-- | plugins/groebner/GroebnerR.v | 162 | ||||
-rw-r--r-- | plugins/groebner/groebner.ml4 | 294 | ||||
-rw-r--r-- | plugins/groebner/ideal.ml | 930 | ||||
-rw-r--r-- | plugins/groebner/ideal.mli | 88 | ||||
-rw-r--r-- | plugins/groebner/polynom.ml | 413 | ||||
-rw-r--r-- | plugins/groebner/polynom.mli | 227 | ||||
-rw-r--r-- | plugins/groebner/utile.ml | 2 | ||||
-rw-r--r-- | plugins/setoid_ring/Field_tac.v | 4 | ||||
-rw-r--r-- | tactics/tacinterp.ml | 2 |
9 files changed, 1430 insertions, 692 deletions
diff --git a/plugins/groebner/GroebnerR.v b/plugins/groebner/GroebnerR.v index 71e244126..0ab135500 100644 --- a/plugins/groebner/GroebnerR.v +++ b/plugins/groebner/GroebnerR.v @@ -33,8 +33,8 @@ Require Import Setoid. Require Import BinPos. Require Import BinList. Require Import Znumtheory. -Require Import Ring_polynom Ring_tac InitialRing. Require Import RealField Rdefinitions Rfunctions RIneq DiscrR. +Require Import Ring_polynom Ring_tac InitialRing. Declare ML Module "groebner_plugin". @@ -49,16 +49,6 @@ Lemma psos_r1: forall x y, x = y -> x - y = 0. intros x y H; rewrite H; ring. Qed. -Ltac equalities_to_goal := - match goal with -| H: (@eq R ?x 0) |- _ => - try generalize H; clear H -| H: (@eq R 0 ?x) |- _ => - try generalize (sym_equal H); clear H -| H: (@eq R ?x ?y) |- _ => - try generalize (psos_r1 _ _ H); clear H - end. - Lemma groebnerR_not1: forall x y:R, x<>y -> exists z:R, z*(x-y)-1=0. intros. exists (1/(x-y)). @@ -80,12 +70,26 @@ field. auto. Qed. + +Ltac equalities_to_goal := + lazymatch goal with + | H: (@eq R ?x 0) |- _ => try revert H + | H: (@eq R 0 ?x) |- _ => + try generalize (sym_equal H); clear H + | H: (@eq R ?x ?y) |- _ => + try generalize (psos_r1 _ _ H); clear H + end. + Lemma groebnerR_not2: 1<>0. auto with *. Qed. Lemma groebnerR_diff: forall x y:R, x<>y -> x-y<>0. -Admitted. +intros. +intro; apply H. +replace x with (x-y+y) by ring. +rewrite H0; ring. +Qed. (* Removes x<>0 from hypothesis *) Ltac groebnerR_not_hyp:= @@ -99,11 +103,11 @@ Ltac groebnerR_not_hyp:= |_ => generalize (@groebnerR_diff _ _ H); clear H; intro end end. - + Ltac groebnerR_not_goal := match goal with - | |- ?x<>?y => red; intro; apply groebnerR_not2 - | |- False => apply groebnerR_not2 + | |- ?x<>?y :> R => red; intro; apply groebnerR_not2 + | |- False => apply groebnerR_not2 end. Ltac groebnerR_begin := @@ -234,13 +238,18 @@ Qed. (* fin du code de Benjamin *) -Lemma groebnerR_l3:forall c p:R, ~c=0 -> c*p=0 -> p=0. +Lemma groebnerR_l3:forall c p r, ~c=0 -> c*p^r=0 -> p=0. intros. elim (Rmult_integral _ _ H0);intros. -absurd (c=0);auto. - auto. + absurd (c=0);auto. + + clear H0; induction r; simpl in *. + contradict H1; discrR. + + elim (Rmult_integral _ _ H1); auto. Qed. + Ltac generalise_eq_hyps:= repeat (match goal with @@ -293,7 +302,7 @@ Fixpoint interpret3 t fv {struct t}: R := Ltac parametres_en_tete fv lp := match fv with - | (@nil _) => lp + | (@nil _) => lp | (@cons _ ?x ?fv1) => let res := AddFvTail x lp in parametres_en_tete fv1 res @@ -301,23 +310,43 @@ Ltac parametres_en_tete fv lp := Ltac append1 a l := match l with - | (@nil _) => constr:(cons a l) + | (@nil _) => constr:(cons a l) | (cons ?x ?l) => let l' := append1 a l in constr:(cons x l') end. Ltac rev l := match l with - |(@nil _) => l + |(@nil _) => l | (cons ?x ?l) => let l' := rev l in append1 x l' end. -(* Pompe sur Ring *) -Import Ring_tac. + +Ltac groebner_call_n nparam p rr lp kont := + groebner_compute (PEc nparam :: PEpow p rr :: lp); + match goal with + | |- (?c::PEpow _ ?r::?lq0)::?lci0 = _ -> _ => + intros _; + set (lci:=lci0); + set (lq:=lq0); + kont c rr lq lci + end. + +Ltac groebner_call nparam p lp kont := + let rec try_n n := + lazymatch n with + | 6%N => fail + | _ => +(* idtac "Trying power: " n;*) + groebner_call_n nparam p n lp kont || + let n' := eval compute in (Nsucc n) in try_n n' + end in + try_n 1%N. + Ltac groebnerR_gen lparam lvar n RNG lH _rl := get_Pre RNG (); - let mkFV := get_RingFV RNG in - let mkPol := get_RingMeta RNG in + let mkFV := Ring_tac.get_RingFV RNG in + let mkPol := Ring_tac.get_RingMeta RNG in generalise_eq_hyps; let t := Get_goal in let lpol := lpol_goal t in @@ -342,26 +371,16 @@ Ltac groebnerR_gen lparam lvar n RNG lH _rl := (@nil (PExpr Z)) lpol in let lpol := eval compute in (List.rev lpol) in - let lpol := constr:(PEc nparam :: lpol) in (*idtac lpol;*) - groebner_compute lpol; - let OnResult kont := - match goal with - | |- (?p ::_)::(?c::?lq0)::?lci0 = _ -> _ => - intros _; - set (lci:=lci0); - set (lq:=lq0); - kont p c lq lci - end in let SplitPolyList kont := match lpol with - | _::?p2::?lp2 => kont p2 lp2 + | ?p2::?lp2 => kont p2 lp2 end in - OnResult ltac:(fun p c lq lci => - SplitPolyList ltac:(fun p2 lp2 => - set (p21:=p2) ; - set (lp21:=lp2); - set (q := PEmul c p21); + SplitPolyList ltac:(fun p lp => + set (p21:=p) ; + set (lp21:=lp); + groebner_call nparam p lp ltac:(fun c r lq lci => + set (q := PEmul c (PEpow p21 r)); let Hg := fresh "Hg" in assert (Hg:check lp21 q (lci,lq) = true); [ (vm_compute;reflexivity) || idtac "invalid groebner certificate" @@ -370,7 +389,7 @@ Ltac groebnerR_gen lparam lvar n RNG lH _rl := [ simpl; apply (@check_correct fv lp21 q (lci,lq) Hg); simpl; repeat (split;[assumption|idtac]); exact I | simpl in Hg2; simpl; - apply groebnerR_l3 with (interpret3 c fv);simpl; + apply groebnerR_l3 with (interpret3 c fv) (Nnat.nat_of_N r);simpl; [ discrR || idtac "could not prove discrimination result" | exact Hg2] ] @@ -391,13 +410,67 @@ Ltac groebnerRp lparam := groebnerRpv lparam (@nil R). (*************** Examples *) (* +Section Examples. + +Require Import List Ring_polynom. +Delimit Scope PE_scope with PE. +Infix "+" := PEadd : PE_scope. +Infix "*" := PEmul : PE_scope. +Infix "-" := PEsub : PE_scope. +Infix "^" := PEpow : PE_scope. +Notation "[ n ]" := (@PEc Z n) (at level 0). + Open Scope R_scope. + Goal forall x y, x+y=0 -> x*y=0 -> x^2=0. +groebnerR. +Qed. + +Goal forall x, x^2=0 -> x=0. +groebnerR. +Qed. + +(* +Notation X := (PEX Z 3). +Notation Y := (PEX Z 2). +Notation Z_ := (PEX Z 1). +*) +Goal forall x y z, + x+y+z=0 -> + x*y+x*z+y*z=0-> + x*y*z=0 -> x^3=0. Time groebnerR. Qed. + +(* +Notation X := (PEX Z 4). +Notation Y := (PEX Z 3). +Notation Z_ := (PEX Z 2). +Notation U := (PEX Z 1). +*) +Goal forall x y z u, + x+y+z+u=0 -> + x*y+x*z+x*u+y*z+y*u+z*u=0-> + x*y*z+x*y*u+x*z*u+y*z*u=0-> + x*y*z*u=0 -> x^4=0. +Time groebnerR. +Qed. + +(* +Notation x_ := (PEX Z 5). +Notation y_ := (PEX Z 4). +Notation z_ := (PEX Z 3). +Notation u_ := (PEX Z 2). +Notation v_ := (PEX Z 1). +Notation "x :: y" := (List.cons x y) +(at level 60, right associativity, format "'[hv' x :: '/' y ']'"). +Notation "x :: y" := (List.app x y) +(at level 60, right associativity, format "x :: y"). +*) + Goal forall x y z u v, x+y+z+u+v=0 -> x*y+x*z+x*u+x*v+y*z+y*u+y*v+z*u+z*v+u*v=0-> @@ -406,10 +479,9 @@ Goal forall x y z u v, x*y*z*u*v=0 -> x^5=0. Time groebnerR. Qed. -*) - - +End Examples. +*) diff --git a/plugins/groebner/groebner.ml4 b/plugins/groebner/groebner.ml4 index 00cd8ae5d..da41a89b6 100644 --- a/plugins/groebner/groebner.ml4 +++ b/plugins/groebner/groebner.ml4 @@ -36,7 +36,10 @@ open Entries open Num open Unix open Utile -open Ideal + +(*********************************************************************** + 1. Opérations sur les coefficients. +*) let num_0 = Int 0 and num_1 = Int 1 @@ -48,6 +51,94 @@ let numdom r = num_of_big_int(Ratio.numerator_ratio r'), num_of_big_int(Ratio.denominator_ratio r') +module BigInt = struct + open Big_int + + type t = big_int + let of_int = big_int_of_int + let coef0 = of_int 0 + let coef1 = of_int 1 + let of_num = Num.big_int_of_num + let to_num = Num.num_of_big_int + let equal = eq_big_int + let lt = lt_big_int + let le = le_big_int + let abs = abs_big_int + let plus =add_big_int + let mult = mult_big_int + let sub = sub_big_int + let opp = minus_big_int + let div = div_big_int + let modulo = mod_big_int + let to_string = string_of_big_int + let to_int x = int_of_big_int x + let hash x = + try (int_of_big_int x) + with _-> 1 + let puis = power_big_int_positive_int + + (* a et b positifs, résultat positif *) + let rec pgcd a b = + if equal b coef0 + then a + else if lt a b then pgcd b a else pgcd b (modulo a b) + + + (* signe du pgcd = signe(a)*signe(b) si non nuls. *) + let pgcd2 a b = + if equal a coef0 then b + else if equal b coef0 then a + else let c = pgcd (abs a) (abs b) in + if ((lt coef0 a)&&(lt b coef0)) + ||((lt coef0 b)&&(lt a coef0)) + then opp c else c +end + +(* +module Ent = struct + type t = Entiers.entiers + let of_int = Entiers.ent_of_int + let of_num x = Entiers.ent_of_string(Num.string_of_num x) + let to_num x = Num.num_of_string (Entiers.string_of_ent x) + let equal = Entiers.eq_ent + let lt = Entiers.lt_ent + let le = Entiers.le_ent + let abs = Entiers.abs_ent + let plus =Entiers.add_ent + let mult = Entiers.mult_ent + let sub = Entiers.moins_ent + let opp = Entiers.opp_ent + let div = Entiers.div_ent + let modulo = Entiers.mod_ent + let coef0 = Entiers.ent0 + let coef1 = Entiers.ent1 + let to_string = Entiers.string_of_ent + let to_int x = Entiers.int_of_ent x + let hash x =Entiers.hash_ent x + let signe = Entiers.signe_ent + + let rec puis p n = match n with + 0 -> coef1 + |_ -> (mult p (puis p (n-1))) + + (* a et b positifs, résultat positif *) + let rec pgcd a b = + if equal b coef0 + then a + else if lt a b then pgcd b a else pgcd b (modulo a b) + + + (* signe du pgcd = signe(a)*signe(b) si non nuls. *) + let pgcd2 a b = + if equal a coef0 then b + else if equal b coef0 then a + else let c = pgcd (abs a) (abs b) in + if ((lt coef0 a)&&(lt b coef0)) + ||((lt coef0 b)&&(lt a coef0)) + then opp c else c +end +*) + (* ------------------------------------------------------------------------- *) (* term?? *) (* ------------------------------------------------------------------------- *) @@ -59,10 +150,24 @@ type term = | Const of Num.num | Var of vname | Opp of term - | Add of (term * term) - | Sub of (term * term) - | Mul of (term * term) - | Pow of (term * int) + | Add of term * term + | Sub of term * term + | Mul of term * term + | Pow of term * int + +let const n = + if eq_num n num_0 then Zero else Const n +let pow(p,i) = if i=1 then p else Pow(p,i) +let add = function + (Zero,q) -> q + | (p,Zero) -> p + | (p,q) -> Add(p,q) +let mul = function + (Zero,_) -> Zero + | (_,Zero) -> Zero + | (p,Const n) when eq_num n num_1 -> p + | (Const n,q) when eq_num n num_1 -> q + | (p,q) -> Mul(p,q) let unconstr = mkRel 1 @@ -201,6 +306,17 @@ let string_of_term p = | Pow (t1,n) -> (aux t1)^"^"^(string_of_int n) in aux p + +(*********************************************************************** + Coefficients: polynomes recursifs + *) + +module Coef = BigInt +(*module Coef = Ent*) +module Poly = Polynom.Make(Coef) +module PIdeal = Ideal.Make(Poly) +open PIdeal + (* terme vers polynome creux *) (* les variables d'indice <=np sont dans les coeff *) @@ -212,11 +328,11 @@ let term_pol_sparse np t= | Const r -> if r = num_0 then zeroP - else polconst d (Polynom.Pint (Polynom.coef_of_num r)) + else polconst d (Poly.Pint (Coef.of_num r)) | Var v -> let v = int_of_string v in if v <= np - then polconst d (Polynom.x v) + then polconst d (Poly.x v) else gen d v | Opp t1 -> oppP (aux t1) | Add (t1,t2) -> plusP d (aux t1) (aux t2) @@ -230,16 +346,17 @@ let term_pol_sparse np t= (* polynome recusrsif vers terme *) + let polrec_to_term p = let rec aux p = match p with - |Polynom.Pint n -> Const (Polynom.num_of_coef n) - |Polynom.Prec (v,coefs) -> + |Poly.Pint n -> const (Coef.to_num n) + |Poly.Prec (v,coefs) -> let res = ref Zero in Array.iteri (fun i c -> - res:=Add(!res, Mul(aux c, - Pow (Var (string_of_int v), + res:=add(!res, mul(aux c, + pow (Var (string_of_int v), i)))) coefs; !res @@ -248,52 +365,50 @@ let polrec_to_term p = (* on approche la forme de Horner utilisee par la tactique ring. *) let pol_sparse_to_term n2 p = + let p = PIdeal.repr p in let rec aux p = match p with - |[] -> Const (num_of_string "0") - |(a,m)::p1 -> - let n = (Array.length m)-1 in - let (i0,e0) = - (List.fold_left - (fun (r,d) (a,m) -> - let i0= ref 0 in - for k=1 to n do - if m.(k)>0 - then i0:=k - done; - if !i0 = 0 - then (r,d) - else if !i0 > r - then (!i0, m.(!i0)) - else if !i0 = r && m.(!i0)<d - then (!i0, m.(!i0)) - else (r,d)) - (0,0) - p) in - if i0=0 - then - let mp = ref (polrec_to_term a) in - if p1=[] - then !mp - else Add(!mp,aux p1) - else ( - let p1=ref [] in - let p2=ref [] in - List.iter - (fun (a,m) -> - if m.(i0)>=e0 - then (m.(i0)<-m.(i0)-e0; - p1:=(a,m)::(!p1)) - else p2:=(a,m)::(!p2)) - p; - let vm = - if e0=1 - then Var (string_of_int (i0)) - else Pow (Var (string_of_int (i0)),e0) in - Add(Mul(vm, aux (List.rev (!p1))), aux (List.rev (!p2)))) - in let pp = aux p in - - pp + [] -> const (num_of_string "0") + | (a,m)::p1 -> + let n = (Array.length m)-1 in + let (i0,e0) = + List.fold_left (fun (r,d) (a,m) -> + let i0= ref 0 in + for k=1 to n do + if m.(k)>0 + then i0:=k + done; + if !i0 = 0 + then (r,d) + else if !i0 > r + then (!i0, m.(!i0)) + else if !i0 = r && m.(!i0)<d + then (!i0, m.(!i0)) + else (r,d)) + (0,0) + p in + if i0=0 + then + let mp = ref (polrec_to_term a) in + if p1=[] + then !mp + else add(!mp,aux p1) + else ( + let p1=ref [] in + let p2=ref [] in + List.iter + (fun (a,m) -> + if m.(i0)>=e0 + then (m.(i0)<-m.(i0)-e0; + p1:=(a,m)::(!p1)) + else p2:=(a,m)::(!p2)) + p; + let vm = + if e0=1 + then Var (string_of_int (i0)) + else pow (Var (string_of_int (i0)),e0) in + add(mul(vm, aux (List.rev (!p1))), aux (List.rev (!p2)))) + in aux p let rec remove_list_tail l i = @@ -360,12 +475,10 @@ let remove_zeros zero lci = let theoremedeszeros lpol p = let t1 = Unix.gettimeofday() in let m = !nvars in - let (c,lp0,p,lci,lc) = in_ideal m lpol p in + let (lp0,p,cert) = in_ideal m lpol p in let lpc = List.rev !poldepcontent in info ("temps: "^Format.sprintf "@[%10.3f@]s\n" (Unix.gettimeofday ()-.t1)); - if lc<>[] - then (c,lp0,p,lci,lc,1,lpc) - else (c,[],zeroP,[],[],0,[]) + (cert,lp0,p,lpc) let theoremedeszeros_termes lp = @@ -383,33 +496,46 @@ let theoremedeszeros_termes lp = | [] -> assert false | p::lp1 -> let lpol = List.rev lp1 in - let (c,lp0,p,lci,lq,d,lct)= theoremedeszeros lpol p in - - let lci1 = List.rev lci in - match remove_zeros (fun x -> x=zeroP) (lq::lci1) with + let (cert,lp0,p,_lct) = theoremedeszeros lpol p in + let lc = cert.last_comb::List.rev cert.gb_comb in + match remove_zeros (fun x -> x=zeroP) lc with | [] -> assert false | (lq::lci) -> (* lci commence par les nouveaux polynomes *) let m= !nvars in + let c = pol_sparse_to_term m (polconst m cert.coef) in + let r = Pow(Zero,cert.power) in let lci = List.rev lci in - let c = pol_sparse_to_term m [c, Array.create (m+1) 0] in - let lp0 = List.map (pol_sparse_to_term m) lp0 in - let p = (Pow ((pol_sparse_to_term m p),d)) in - let lci = List.map (fun lc -> List.map (pol_sparse_to_term m) lc) lci in + let lci = List.map (List.map (pol_sparse_to_term m)) lci in let lq = List.map (pol_sparse_to_term m) lq in info ("nombre de parametres: "^string_of_int nparam^"\n"); info "terme calcule\n"; - (c,lp0,p,lci,lq) + (c,r,lci,lq) ) |_ -> assert false - +(* version avec hash-consing du certificat: +let groebner lpol = + Hashtbl.clear Dansideal.hmon; + Hashtbl.clear Dansideal.coefpoldep; + Hashtbl.clear Dansideal.sugartbl; + Hashtbl.clear Polynomesrec.hcontentP; + init_constants (); + let lp= parse_request lpol in + let (_lp0,_p,c,r,_lci,_lq as rthz) = theoremedeszeros_termes lp in + let certif = certificat_vers_polynome_creux rthz in + let certif = hash_certif certif in + let certif = certif_term certif in + let c = mkt_term c in + info "constr calcule\n"; + (c, certif) +*) let groebner lpol = let lp= parse_request lpol in - let (c,lp0,p,lci,lq) = theoremedeszeros_termes lp in - let res = [p::lp0;c::lq]@lci in + let (c,r,lci,lq) = theoremedeszeros_termes lp in + let res = [c::r::lq]@lci in let res = List.map (fun lx -> List.map mkt_term lx) res in let res = List.fold_right @@ -426,24 +552,20 @@ let groebner lpol = info "terme calcule\n"; res -(* -open Util -open Term -open Tactics -open Coqlib -open Utile -open Groebner -*) - -let groebner_compute t gl = - let lpol = groebner t in +let return_term t = let a = - mkApp(gen_constant "CC" ["Init";"Logic"] "refl_equal",[|tllp ();lpol|]) in - generalize [a] gl + mkApp(gen_constant "CC" ["Init";"Logic"] "refl_equal",[|tllp ();t|]) in + generalize [a] + +let groebner_compute t = + let lpol = + try groebner t + with Ideal.NotInIdeal -> + error "groebner cannot solve this problem" in + return_term lpol TACTIC EXTEND groebner_compute -| [ "groebner_compute" constr(lt) ] -> - [ groebner_compute lt ] +| [ "groebner_compute" constr(lt) ] -> [ groebner_compute lt ] END diff --git a/plugins/groebner/ideal.ml b/plugins/groebner/ideal.ml index a87a9972b..4570f69d8 100644 --- a/plugins/groebner/ideal.ml +++ b/plugins/groebner/ideal.ml @@ -21,50 +21,101 @@ open Utile -(*********************************************************************** - Coefficients: polynomes recursifs - on n'utilise pas les big_int car leur egalite n'est pas generique: on ne peut pas utiliser de tables de hash simplement. - *) +exception NotInIdeal -type coef = Polynom.poly -let coef_of_int = Polynom.cf -let eq_coef = Polynom.eqP -let plus_coef =Polynom.plusP -let mult_coef = Polynom.multP -let sub_coef = Polynom.moinsP -let opp_coef = Polynom.oppP -let div_coef = Polynom.divP (* division exacte *) -let coef0 = Polynom.cf0 -let coef1 = Polynom.cf1 -let string_of_coef c = "["^(Polynom.string_of_P c)^"]" +module type S = sig +(* Monomials *) +type mon = int array -let rec pgcd a b = Polynom.pgcdP a b +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 -(* -let htblpgcd = Hashtbl.create 1000 - -let pgcd a b = - try - (let c = Hashtbl.find htblpgcd (a,b) in - info "*"; - c) - with _ -> - (let c = Polynom.pgcdP a b in - info "-"; - Hashtbl.add htblpgcd (a,b) c; - c) +(*********************************************************************** + 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 *) -let hmon = Hashtbl.create 1000 - type mon = int array +type deg = int +type poly = (coef * mon) list (* d représente le nb de variables du monome *) @@ -80,30 +131,30 @@ let mult_mon d m m' = let deg m = m.(0) -(* Comparaison de monomes avec ordre du degre lexicographique = on commence par regarder la 1ere variable*) -(*let compare_mon d m m' = - 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 -*) - -(* degre lexicographique inverse *) let compare_mon d m m' = - 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 - + 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' = @@ -114,7 +165,7 @@ let div_mon d m m' = m'' let div_pol_coef p c = - List.map (fun (a,m) -> (div_coef a c,m)) p + List.map (fun (a,m) -> (P.divP a c,m)) p (* m' divise m *) let div_mon_test d m m' = @@ -153,11 +204,25 @@ let ppcm_mon d m m' = *) -type poly = (coef * mon) list -let eq_poly = +let repr p = p + +let equal = Util.list_for_all2eq - (fun (c1,m1) (c2,m2) -> eq_coef c1 c2 && m1=m2) + (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. @@ -337,8 +402,8 @@ let plusP d p q = match compare_mon d (snd t) (snd t') with 1 -> t::(plusP p' q) |(-1) -> t'::(plusP p q') - |_ -> let c=plus_coef (fst t) (fst t') in - match eq_coef c coef0 with + |_ -> 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 @@ -349,7 +414,7 @@ let mult_t_pol d a m p = let rec mult_t_pol p = match p with [] -> [] - |(b,m')::p -> ((mult_coef a b),(mult_mon d m m'))::(mult_t_pol p) + |(b,m')::p -> ((P.multP a b),(mult_mon d m m'))::(mult_t_pol p) in mult_t_pol p @@ -357,10 +422,9 @@ let mult_t_pol d a m p = let rec selectdiv d m l = match l with [] -> [] - |q::r -> let m'= snd (List.hd q) in - match (div_mon_test d m m') with - true -> q - |false -> selectdiv d m r + | 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 *) @@ -368,7 +432,7 @@ let gen d i = let m = Array.create (d+1) 0 in m.(i) <- 1; let m = set_deg d m in - [((coef_of_int 1),m)] + [(coef1,m)] @@ -377,7 +441,7 @@ let oppP p = let rec oppP p = match p with [] -> [] - |(b,m')::p -> ((opp_coef b),m')::(oppP p) + |(b,m')::p -> ((P.oppP b),m')::(oppP p) in oppP p @@ -386,7 +450,7 @@ let emultP a p = let rec emultP p = match p with [] -> [] - |(b,m')::p -> ((mult_coef a b),m')::(emultP p) + |(b,m')::p -> ((P.multP a b),m')::(emultP p) in emultP p @@ -406,23 +470,41 @@ let puisP d p n= |_ -> 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 -> pgcd a (contentP p1) + |(a,m)::p1 -> pgcd1 a (contentP p1) let contentPlist lp = match lp with - |[] -> coef1 - |p::l1 -> List.fold_left (fun r q -> pgcd r (contentP q)) (contentP p) l1 + |[] -> coef1 + |p::l1 -> List.fold_left (fun r q -> pgcd1 r (contentP q)) (contentP p) l1 (*********************************************************************** Division de polynomes *) -let pgcdpos a b = pgcd a b +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 ".";*) @@ -438,28 +520,23 @@ let reduce2 d p l = let rec reduce p = match p with [] -> (coef1,[]) - |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 (c,r)=(reduce p') in - (c,((mult_coef a c,m)::r)) - else (coef1,p) + | (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'= (div_coef b c) in - let b'=(opp_coef (div_coef a c)) 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 - (mult_coef a' e,r) - in let (c,r) = reduce p in - (c,r) + (P.multP a' e,r)) in + reduce p (* trace des divisions *) @@ -467,22 +544,33 @@ let reduce2 d p l = 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 = Hashtbl.create 51 +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 (Hashtbl.find coefpoldep (p,q)) + try (HashPolPair.find coefpoldep (p,q)) with _ -> [] let coefpoldep_set p q c = - Hashtbl.add coefpoldep (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 (coef_of_int 1))) + (fun p -> coefpoldep_set p p (polconst d coef1)) lp (* garde la trace dans coefpoldep @@ -509,7 +597,7 @@ let reduce2_trace d p l lcp = (lq,((a,m)::r)) else ([],p) |(b,m')::q' -> - let b'=(opp_coef (div_coef a b)) in + 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 @@ -527,9 +615,9 @@ let reduce2_trace d p l lcp = let c = List.fold_left (fun x (a,m,s) -> - if eq_poly s q + if equal s q then - plusP d x (mult_t_pol d a m (polconst d (coef_of_int 1))) + plusP d x (mult_t_pol d a m (polconst d coef1)) else x) c0 lq in @@ -539,9 +627,457 @@ let reduce2_trace d p l lcp = 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)::_ = 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 homogeneous = ref false + +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 @@ -552,17 +1088,17 @@ let spol d p q= 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 - (div_coef b c) - (div_mon d m'' m) p') - (mult_t_pol d - (opp_coef (div_coef a c)) - (div_mon d m'' m') q') in + (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 (coef_of_int 1)) []); - coefpoldep_set sp q (fsp [] (polconst d (coef_of_int 1))); + 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'= @@ -599,48 +1135,47 @@ type 'poly cpRes = Keep of ('poly list) | DontKeep of ('poly list) -let list_rec f0 f1 = - let rec f2 = function - [] -> f0 - | a0::l0 -> f1 a0 l0 (f2 l0) - in f2 - let addRes i = function Keep h'0 -> Keep (i::h'0) | DontKeep h'0 -> DontKeep (i::h'0) -let slice d i a q = - list_rec - (match etrangers d i a with - true -> DontKeep [] - | false -> Keep []) - (fun b q1 rec_ren -> - match div_ppcm d i a b with - true -> DontKeep (b::q1) - | false -> - (match div_ppcm d i b a with - true -> rec_ren - | false -> addRes b rec_ren)) q +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 genPcPf d i aP q = - (let rec genPc aP0 = - match aP0 with - [] -> (fun r -> r) +let rec genPcPf d i aP q = + match aP with + [] -> q | a::l1 -> - (fun l -> - (match slice d i a l1 with - Keep l2 -> addS (spol d i a) (genPc l2 l) - | DontKeep l2 -> genPc l2 l)) - in genPc aP) q + (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 genOCPf d h' = - list_rec [] (fun a l rec_ren -> - genPcPf d a l rec_ren) h' +let rec genOCPf d = function + [] -> [] + | a::l -> genPcPf d a l (genOCPf d l) let step = ref 0 @@ -654,12 +1189,16 @@ let infobuch p q = (* dans lp les nouveaux sont en queue *) let coef_courant = ref coef1 -let pol_courant = ref [] + + +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:= mult_coef !coef_courant c; + coef_courant:= P.multP !coef_courant c; pol_courant:= r; if r=[] then (info "polynome reduit a 0\n"; @@ -673,7 +1212,7 @@ let test_dans_ideal d p lp lp0 = ) lcq !poldep; info ("verif somme: "^(stringP (!res))^"\n"); - + info ("coefficient: "^(stringP (polconst d c))^"\n"); let rec aux lp = match lp with |[] -> [] @@ -682,7 +1221,6 @@ let test_dans_ideal d p lp lp0 = (fun q -> coefpoldep_find p q) lp)::(aux lp) in - let coefficient_multiplicateur = c in let liste_polynomes_de_depart = List.rev lp0 in let polynome_a_tester = p in let liste_des_coefficients_intermediaires = @@ -692,68 +1230,71 @@ let test_dans_ideal d p lp lp0 = !lci) in let liste_des_coefficients = List.map - (fun cq -> emultP (coef_of_int (-1)) cq) + (fun cq -> emultP coefm1 cq) (List.rev lcq) in - (coefficient_multiplicateur, - liste_polynomes_de_depart, + (liste_polynomes_de_depart, polynome_a_tester, - liste_des_coefficients_intermediaires, - liste_des_coefficients)) + {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");*) - (c,[],[],[],[])) - -let choix_spol p l = l - -let deg_hom p = - match p with - | [] -> -1 - | (a,m)::_ -> m.(0) + 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)::_ = p in + try ( + let q = + List.find + (fun q -> + let (_,m')::_ = 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 x = - let (lp, lpc) = x in + let rec pbuchf lp lpc = infobuch lp lpc; (* step:=(!step+1)mod 10;*) match lpc with [] -> test_dans_ideal d p lp lp0 | _ -> - (match choix_spol !pol_courant lpc with - |[] -> assert false - |(a::lpc2) -> + 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)) + then (info "h";pbuchf lp lpc2) else - 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");*) - poldep:=addS a0 lp; - poldepcontent:=addS ct (!poldepcontent); - - let (c1,lp1,p1,lci1,lc1) - = test_dans_ideal d p (addS a0 lp) lp0 in - if lc1<>[] - then (c1,lp1,p1,lci1,lc1) - else - pbuchf - (((addS a0 lp), - (genPcPf d a0 lp lpc2)))) - in pbuchf pq + 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 @@ -777,16 +1318,15 @@ let is_homogeneous p = *) let in_ideal d lp p = - info ("p: "^(stringPcut p)^"\n"); - info ("lp:\n"^(List.fold_left (fun r p -> r^(stringPcut p)^"\n") "" lp)); + 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; - homogeneous := List.for_all is_homogeneous (p::lp); - if !homogeneous then info "polynomes homogenes\n"; - let (c1,lp1,p1,lci1,lc1) = test_dans_ideal d p lp lp in - if lc1<>[] - then (c1,lp1,p1,lci1,lc1) - else pbuchf d (lp, (genOCPf d lp)) p lp - + try test_dans_ideal d p lp lp + with NotInIdeal -> pbuchf d (lp, genOCPf d lp) p lp +end diff --git a/plugins/groebner/ideal.mli b/plugins/groebner/ideal.mli index 48b81102e..9c8f43d7d 100644 --- a/plugins/groebner/ideal.mli +++ b/plugins/groebner/ideal.mli @@ -6,77 +6,75 @@ (* * GNU Lesser General Public License Version 2.1 *) (************************************************************************) -type coef = Polynom.poly -val coef_of_int : int -> coef -val eq_coef : coef -> coef -> bool -val plus_coef : coef -> coef -> coef -val mult_coef : coef -> coef -> coef -val sub_coef : coef -> coef -> coef -val opp_coef : coef -> coef -val div_coef : coef -> coef -> coef -val coef0 : coef -val coef1 : coef -val string_of_coef : coef -> string -val pgcd : coef -> coef -> coef +exception NotInIdeal +val lexico : bool ref +val use_hmon : bool ref +module type S = sig + +(* Monomials *) type mon = int array -type poly = (coef * mon) list -val eq_poly : poly -> poly -> bool 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_pol_coef : poly -> coef -> poly val div_mon_test : int -> mon -> mon -> bool -val set_deg : int -> mon -> mon 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 stringP : poly -> string -val stringPcut : poly -> string val lstringP : poly list -> string val printP : poly -> unit val lprintP : poly list -> unit -val polconst : int -> coef -> poly -val zeroP : poly -val plusP : int -> poly -> poly -> poly -val mult_t_pol : int -> coef -> mon -> poly -> poly -val selectdiv : int -> mon -> poly list -> poly -val gen : int -> int -> poly +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 : int -> poly -> poly -> poly -val puisP : int -> poly -> int -> 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 : int -> poly -> poly -> coef -> coef -> mon -> poly -val reduce2 : int -> poly -> poly list -> coef * poly +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 : int -> poly list -> unit -val reduce2_trace : - int -> poly -> poly list -> poly list -> poly list * poly -val homogeneous : bool ref -val spol : int -> poly -> poly -> poly -val etrangers : int -> ('a * mon) list -> ('b * mon) list -> bool -val div_ppcm : int -> poly -> poly -> poly -> bool +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 -type 'a cpRes = Keep of 'a list | DontKeep of 'a list -val list_rec : 'a -> ('b -> 'b list -> 'a -> 'a) -> 'b list -> 'a -val addRes : 'a -> 'a cpRes -> 'a cpRes -val slice : int -> poly -> poly -> poly list -> poly cpRes -val addS : 'a -> 'a list -> 'a list -val genPcPf : int -> poly -> poly list -> poly list -> poly list -val genOCPf : int -> poly list -> poly list +val genPcPf : deg -> poly -> poly list -> poly list -> poly list +val genOCPf : deg -> poly list -> poly list val is_homogeneous : poly -> bool -val in_ideal : - int -> poly list -> poly -> - coef * poly list * poly * poly list list * poly list +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 + +module Make (P:Polynom.S) : S with type coef = P.t diff --git a/plugins/groebner/polynom.ml b/plugins/groebner/polynom.ml index 5c859f755..6d2ed26e8 100644 --- a/plugins/groebner/polynom.ml +++ b/plugins/groebner/polynom.ml @@ -1,5 +1,4 @@ (************************************************************************) - (* v * The Coq Proof Assistant / The Coq Development Team *) (* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) (* \VV/ **************************************************************) @@ -13,89 +12,130 @@ open Utile open Util +module type Coef = sig + type t + val equal : t -> t -> bool + val lt : t -> t -> bool + val le : t -> t -> bool + val abs : t -> t + val plus : t -> t -> t + val mult : t -> t -> t + val sub : t -> t -> t + val opp : t -> t + val div : t -> t -> t + val modulo : t -> t -> t + val puis : t -> int -> t + val pgcd : t -> t -> t + + val hash : t -> int + val of_num : Num.num -> t + val to_string : t -> string +end + +module type S = sig + type coef + type variable = int + type t = Pint of coef | Prec of variable * t array + + val of_num : Num.num -> t + val x : variable -> t + val monome : variable -> int -> t + val is_constantP : t -> bool + val is_zero : t -> bool + + val max_var_pol : t -> variable + val max_var_pol2 : t -> variable + val max_var : t array -> variable + val equal : t -> t -> bool + val norm : t -> t + val deg : variable -> t -> int + val deg_total : t -> int + val copyP : t -> t + val coef : variable -> int -> t -> t + + val plusP : t -> t -> t + val content : t -> coef + val div_int : t -> coef -> t + val vire_contenu : t -> t + val vars : t -> variable list + val int_of_Pint : t -> coef + val multx : int -> variable -> t -> t + val multP : t -> t -> t + val deriv : variable -> t -> t + val oppP : t -> t + val moinsP : t -> t -> t + val puisP : t -> int -> t + val ( @@ ) : t -> t -> t + val ( -- ) : t -> t -> t + val ( ^^ ) : t -> int -> t + val coefDom : variable -> t -> t + val coefConst : variable -> t -> t + val remP : variable -> t -> t + val coef_int_tete : t -> coef + val normc : t -> t + val coef_constant : t -> coef + val univ : bool ref + val string_of_var : int -> string + val nsP : int ref + val to_string : t -> string + val printP : t -> unit + val print_tpoly : t array -> unit + val print_lpoly : t list -> unit + val quo_rem_pol : t -> t -> variable -> t * t + val div_pol : t -> t -> variable -> t + val divP : t -> t -> t + val div_pol_rat : t -> t -> bool + val pseudo_div : t -> t -> variable -> t * t * int * t + val pgcdP : t -> t -> t + val pgcd_pol : t -> t -> variable -> t + val content_pol : t -> variable -> t + val pgcd_coef_pol : t -> t -> variable -> t + val pgcd_pol_rec : t -> t -> variable -> t + val gcd_sub_res : t -> t -> variable -> t + val gcd_sub_res_rec : t -> t -> t -> t -> int -> variable -> t + val lazard_power : t -> t -> int -> variable -> t + val sans_carre : t -> variable -> t + val facteurs : t -> variable -> t list + val facteurs_impairs : t -> variable -> t list + val hcontentP : (string, t) Hashtbl.t + val prcontentP : unit -> unit + val contentP : t * variable -> t + val hash : t -> int + module Hashpol : Hashtbl.S with type key=t + val memoP : string -> 'a Hashpol.t -> (t -> 'a) -> t -> 'a + val hfactorise : t list list Hashpol.t + val prfactorise : unit -> unit + val factorise : t -> t list list + val facteurs2 : t -> t list + val pol_de_factorisation : t list list -> t + val set_of_array_facteurs : t list array -> t list + val factorise_tableauP2 : + t array -> t list array -> t array * (t * int list) array + val factorise_tableauP : t array -> t array * (t * int list) array + val is_positif : t -> bool + val is_negatif : t -> bool + val pseudo_euclide : + t list -> t -> t -> variable -> + t * t * int * t * t * (t * int) list * (t * int) list + val implique_non_nul : t list -> t -> bool + val ajoute_non_nul : t -> t list -> t list +end + (*********************************************************************** - 1. Opérations sur les coefficients. + 2. Le type des polynômes, operations. *) +module Make (C:Coef) = struct -open Big_int - -type coef = big_int -let coef_of_int = big_int_of_int -let coef_of_num = Num.big_int_of_num -let num_of_coef = Num.num_of_big_int -let eq_coef = eq_big_int -let lt_coef = lt_big_int -let le_coef = le_big_int -let abs_coef = abs_big_int -let plus_coef =add_big_int -let mult_coef = mult_big_int -let sub_coef = sub_big_int -let opp_coef = minus_big_int -let div_coef = div_big_int -let mod_coef = mod_big_int +type coef = C.t +let coef_of_int i = C.of_num (Num.Int i) let coef0 = coef_of_int 0 let coef1 = coef_of_int 1 -let string_of_coef = string_of_big_int -let int_of_coef x = int_of_big_int x -let hash_coef x = - try (int_of_big_int x) - with _-> 1 -let puis_coef = power_big_int_positive_int - -(* -type coef = Entiers.entiers -let coef_of_int = Entiers.ent_of_int -let coef_of_num x = (* error if x is a ratio *) - Entiers.ent_of_string - (Big_int.string_of_big_int (Num.big_int_of_num x)) -let num_of_coef x = Num.num_of_string (Entiers.string_of_ent x) -let eq_coef = Entiers.eq_ent -let lt_coef = Entiers.lt_ent -let le_coef = Entiers.le_ent -let abs_coef = Entiers.abs_ent -let plus_coef =Entiers.add_ent -let mult_coef = Entiers.mult_ent -let sub_coef = Entiers.moins_ent -let opp_coef = Entiers.opp_ent -let div_coef = Entiers.div_ent -let mod_coef = Entiers.mod_ent -let coef0 = Entiers.ent0 -let coef1 = Entiers.ent1 -let string_of_coef = Entiers.string_of_ent -let int_of_coef x = Entiers.int_of_ent x -let hash_coef x =Entiers.hash_ent x -let signe_coef = Entiers.signe_ent - -let rec puis_coef p n = match n with - 0 -> coef1 - |_ -> (mult_coef p (puis_coef p (n-1))) -*) - -(* a et b positifs, résultat positif *) -let rec pgcd a b = - if eq_coef b coef0 - then a - else if lt_coef a b then pgcd b a else pgcd b (mod_coef a b) - - -(* signe du pgcd = signe(a)*signe(b) si non nuls. *) -let pgcd2 a b = - if eq_coef a coef0 then b - else if eq_coef b coef0 then a - else let c = pgcd (abs_coef a) (abs_coef b) in - if ((lt_coef coef0 a)&&(lt_coef b coef0)) - ||((lt_coef coef0 b)&&(lt_coef a coef0)) - then opp_coef c else c - -(*********************************************************************** - 2. Le type des polynômes, operations. -*) type variable = int -type poly = - Pint of coef (* polynome constant *) - | Prec of variable * (poly array) (* coefficients par degre croissant *) +type t = + Pint of coef (* polynome constant *) + | Prec of variable * (t array) (* coefficients par degre croissant *) (* sauf mention du contraire, les opérations ne concernent que des polynomes normalisés: @@ -106,13 +146,38 @@ type poly = *) (* Polynômes constant *) -let cf x = Pint (coef_of_int x) -let cf0 = (cf 0) -let cf1 = (cf 1) +let of_num x = Pint (C.of_num x) +let cf0 = of_num (Num.Int 0) +let cf1 = of_num (Num.Int 1) (* la n-ième variable *) let x n = Prec (n,[|cf0;cf1|]) +(* crée rapidement v^n *) +let monome v n = + match n with + 0->Pint coef1; + |_->let tmp = Array.create (n+1) (Pint coef0) in + tmp.(n)<-(Pint coef1); + Prec (v, tmp) + + +(* teste si un polynome est constant *) +let is_constantP = function + Pint _ -> true + | Prec _ -> false + + +(* conversion d'un poly cst en entier*) +let int_of_Pint = function + Pint x -> x + | _ -> failwith "non" + + +(* teste si un poly est identiquement nul *) +let is_zero p = + match p with Pint n -> if C.equal n coef0 then true else false |_-> false + (* variable max *) let max_var_pol p = match p with @@ -134,13 +199,13 @@ let rec max_var l = Array.fold_right (fun p m -> max (max_var_pol2 p) m) l 0 (* Egalité de deux polynômes On ne peut pas utiliser = car elle ne marche pas sur les Big_int. *) -let rec eqP p q = +let rec equal p q = match (p,q) with - (Pint a,Pint b) -> eq_coef a b + (Pint a,Pint b) -> C.equal a b |(Prec(x,p1),Prec(y,q1)) -> if x<>y then false else if (Array.length p1)<>(Array.length q1) then false - else (try (Array.iteri (fun i a -> if not (eqP a q1.(i)) + else (try (Array.iteri (fun i a -> if not (equal a q1.(i)) then failwith "raté") p1; true) @@ -157,7 +222,7 @@ let rec norm p = match p with |Prec (x,a)-> let d = (Array.length a -1) in let n = ref d in - while !n>0 && (eqP a.(!n) (Pint coef0)) do + while !n>0 && (equal a.(!n) (Pint coef0)) do n:=!n-1; done; if !n<0 then Pint coef0 @@ -201,7 +266,7 @@ let coef v i p = let rec plusP p q = let res = (match (p,q) with - (Pint a,Pint b) -> Pint (plus_coef a b) + (Pint a,Pint b) -> Pint (C.plus a b) |(Pint a, Prec (y,q1)) -> let q2=Array.map copyP q1 in q2.(0)<- plusP p q1.(0); Prec (y,q2) @@ -228,21 +293,21 @@ let rec plusP p q = (* contenu entier positif *) let rec content p = match p with - Pint a -> abs_coef a + Pint a -> C.abs a | Prec (x ,p1) -> - Array.fold_left pgcd coef0 (Array.map content p1) + Array.fold_left C.pgcd coef0 (Array.map content p1) (* divise tous les coefficients de p par l'entier a*) let rec div_int p a= match p with - Pint b -> Pint (div_coef b a) + Pint b -> Pint (C.div b a) | Prec(x,p1) -> Prec(x,Array.map (fun x -> div_int x a) p1) (* divise p par son contenu entier positif. *) let vire_contenu p = let c = content p in - if eq_coef c coef0 then p else div_int p c + if C.equal c coef0 then p else div_int p c (* liste triee des variables impliquees dans un poly *) @@ -251,12 +316,6 @@ let rec vars=function | Prec (x,l)->(List.flatten ([x]::(List.map vars (Array.to_list l)))) -(* conversion d'un poly cst en entier*) -let int_of_Pint=function - (Pint x)->x - |_->failwith "non" - - (* multiplie p par v^n, v >= max_var p *) let rec multx n v p = match p with @@ -274,14 +333,14 @@ let rec multx n v p = (* produit de 2 polynomes *) let rec multP p q = match (p,q) with - (Pint a,Pint b) -> Pint (mult_coef a b) + (Pint a,Pint b) -> Pint (C.mult a b) |(Pint a, Prec (y,q1)) -> - if eq_coef a coef0 then Pint coef0 + if C.equal a coef0 then Pint coef0 else let q2 = Array.map (fun z-> multP p z) q1 in Prec (y,q2) |(Prec (x,p1), Pint b) -> - if eq_coef b coef0 then Pint coef0 + if C.equal b coef0 then Pint coef0 else let p2 = Array.map (fun z-> multP z q) p1 in Prec (x,p2) |(Prec (x,p1), Prec(y,q1)) -> @@ -300,7 +359,7 @@ let rec multP p q = let rec deriv v p = match p with Pint a -> Pint coef0 - |Prec(x,p1) when x=v -> + | Prec(x,p1) when x=v -> let d = Array.length p1 -1 in if d=1 then p1.(1) else @@ -309,13 +368,13 @@ let rec deriv v p = p2.(i)<- multP (Pint (coef_of_int (i+1))) p1.(i+1); done; Prec (x,p2)) - |Prec(x,p1)-> Pint coef0 + | Prec(x,p1)-> Pint coef0 (* opposé de p *) let rec oppP p = match p with - Pint a -> Pint (opp_coef a) + Pint a -> Pint (C.opp a) |Prec(x,p1) -> Prec(x,Array.map oppP p1) @@ -361,26 +420,7 @@ let rec coef_int_tete p = let normc p = let p = vire_contenu p in let a = coef_int_tete p in - if le_coef coef0 a then p else oppP p - - -(* teste si un polynome est constant *) -let is_constantP p = match p with - Pint _ -> true - |Prec (_,_) -> false - - -(* teste si un poly est identiquement nul *) -let is_zero p = - match p with Pint n -> if eq_coef n coef0 then true else false |_-> false - -(* crée rapidement v^n *) -let monome v n = - match n with - 0->Pint coef1; - |_->let tmp = Array.create (n+1) (Pint coef0) in - tmp.(n)<-(Pint coef1); - Prec (v, tmp) + if C.le coef0 a then p else oppP p (*coef constant d'un polynome normalise*) @@ -415,9 +455,9 @@ let rec string_of_Pcut p = else match p with |Pint a-> nsP:=(!nsP)-1; - if le_coef coef0 a - then string_of_coef a - else "("^(string_of_coef a)^")" + if C.le coef0 a + then C.to_string a + else "("^(C.to_string a)^")" |Prec (x,t)-> let v=string_of_var x and s=ref "" @@ -462,16 +502,16 @@ let rec string_of_Pcut p = (s:="0")); !s -let string_of_P p = +let to_string p = nsP:=20; string_of_Pcut p -let printP p = Format.printf "@[%s@]" (string_of_P p) +let printP p = Format.printf "@[%s@]" (to_string p) let print_tpoly lp = let s = ref "\n{ " in - Array.iter (fun p -> s:=(!s)^(string_of_P p)^"\n") lp; + Array.iter (fun p -> s:=(!s)^(to_string p)^"\n") lp; prt0 ((!s)^"}") @@ -488,8 +528,8 @@ let rec quo_rem_pol p q x = if x=0 then (match (p,q) with |(Pint a, Pint b) -> - if eq_coef (mod_coef a b) coef0 - then (Pint (div_coef a b), cf 0) + if C.equal (C.modulo a b) coef0 + then (Pint (C.div a b), cf0) else failwith "div_pol1" |_ -> assert false) else @@ -499,7 +539,7 @@ let rec quo_rem_pol p q x = let r = ref p in let s = ref cf0 in let continue =ref true in - while (!continue) && (not (eqP !r cf0)) do + while (!continue) && (not (equal !r cf0)) do let n = deg x !r in if n<m then continue:=false @@ -517,12 +557,12 @@ let rec quo_rem_pol p q x = (* echoue si q ne divise pas p, rend le quotient sinon *) and div_pol p q x = let (s,r) = quo_rem_pol p q x in - if eqP r cf0 + if equal r cf0 then s else failwith ("div_pol:\n" - ^"p:"^(string_of_P p)^"\n" - ^"q:"^(string_of_P q)^"\n" - ^"r:"^(string_of_P r)^"\n" + ^"p:"^(to_string p)^"\n" + ^"q:"^(to_string q)^"\n" + ^"r:"^(to_string r)^"\n" ^"x:"^(string_of_int x)^"\n" ) @@ -611,11 +651,11 @@ and pgcd_coef_pol c p x = and pgcd_pol_rec p q x = match (p,q) with - (Pint a,Pint b) -> Pint (pgcd (abs_coef a) (abs_coef b)) + (Pint a,Pint b) -> Pint (C.pgcd (C.abs a) (C.abs b)) |_ -> - if eqP p cf0 + if equal p cf0 then q - else if eqP q cf0 + else if equal q cf0 then p else if (deg x q) = 0 then pgcd_coef_pol q p x @@ -634,46 +674,6 @@ and pgcd_pol_rec p q x = res ) -(* version avec memo*) -(* -and htblpgcd = Hashtbl.create 1000 - -and pgcd_pol_rec p q x = - try - (let c = Hashtbl.find htblpgcd (p,q,x) in - (*info "*";*) - c) - with _ -> - (let c = - (match (p,q) with - (Pint a,Pint b) -> Pint (pgcd (abs_coef a) (abs_coef b)) - |_ -> - if eqP p cf0 - then q - else if eqP q cf0 - then p - else if (deg x q) = 0 - then pgcd_coef_pol q p x - else if (deg x p) = 0 - then pgcd_coef_pol p q x - else ( - let a = content_pol p x in - let b = content_pol q x in - let c = pgcd_pol_rec a b (x-1) in - pr (string_of_int x); - let p1 = div_pol p c x in - let q1 = div_pol q c x in - let r = gcd_sub_res p1 q1 x in - let cr = content_pol r x in - let res = c @@ (div_pol r cr x) in - res - ) - ) - in - (*info "-";*) - Hashtbl.add htblpgcd (p,q,x) c; - c) -*) (* Sous-résultants: ai*Ai = Qi*Ai+1 + bi*Ai+2 @@ -692,7 +692,7 @@ and pgcd_pol_rec p q x = *) and gcd_sub_res p q x = - if eqP q cf0 + if equal q cf0 then p else let d = deg x p in @@ -706,7 +706,7 @@ and gcd_sub_res p q x = gcd_sub_res_rec q r (c'^^delta) c' d' x and gcd_sub_res_rec p q s c d x = - if eqP q cf0 + if equal q cf0 then p else ( let d' = deg x q in @@ -768,40 +768,34 @@ let facteurs_impairs p x = (* décomposition sans carrés en toutes les variables *) -let hcontentP = Hashtbl.create 51 +let hcontentP = (Hashtbl.create 51 : (string,t) Hashtbl.t) let prcontentP () = Hashtbl.iter (fun s c -> - prn (s^" -> "^(string_of_P c))) + prn (s^" -> "^(to_string c))) hcontentP let contentP = - memos "c" hcontentP (fun (p,x) -> ((string_of_P p)^":"^(string_of_int x))) + memos "c" hcontentP (fun (p,x) -> ((to_string p)^":"^(string_of_int x))) (fun (p,x) -> content_pol p x) (* Tables de hash et polynômes, mémo *) -(* dans le mli on écrit: - module Hashpol:Hashtbl.S with type t = poly -*) - (* fonction de hachage des polynômes *) -let hashP p = - let rec hP p = - match p with - Pint a -> (hash_coef a) - |Prec (v,p) -> - (Array.fold_right (fun q h -> h+(hP q)) p 0) - in (* Hashtbl.hash *) (hP p) +let rec hash = function + Pint a -> (C.hash a) + | Prec (v,p) -> + Array.fold_right (fun q h -> h + hash q) p 0 module Hashpol = Hashtbl.Make( struct + type poly = t type t = poly - let equal = eqP - let hash = hashP + let equal = equal + let hash = hash end) @@ -817,7 +811,7 @@ let hfactorise = Hashpol.create 51 let prfactorise () = Hashpol.iter (fun p c -> - prn ((string_of_P p)^" -> "); + prn ((to_string p)^" -> "); print_lpoly (List.flatten c)) hfactorise @@ -888,7 +882,7 @@ let factorise_tableauP2 f l1 = (* puis on décompose les polynômes de f avec ces facteurs *) pr "-"; let res = factorise_tableau (fun a b -> div_pol a b x) - (fun p -> eqP p cf0) + (fun p -> equal p cf0) cf0 f l1 in pr ">"; @@ -908,14 +902,14 @@ let rec is_positif p = let res = match p with - Pint a -> le_coef coef0 a + Pint a -> C.le coef0 a |Prec(x,p1) -> (array_for_all is_positif p1) - && (try (Array.iteri (fun i c -> if (i mod 2)<>0 && not (eqP c cf0) + && (try (Array.iteri (fun i c -> if (i mod 2)<>0 && not (equal c cf0) then failwith "pas pair") p1; true) - with _ -> false) + with Failure _ -> false) in res @@ -936,11 +930,11 @@ let pseudo_euclide coef_non_nuls p q x = *) (* vérification de la pseudo-division: let verif = ((c^^d)@@p)--(plusP (s@@q) r) in - if not (eqP verif cf0) - then (prn ("p:"^(string_of_P p)); - prn ("q:"^(string_of_P q)); - prn ("c:"^(string_of_P c)); - prn ("r:"^(string_of_P r)); + if not (equal verif cf0) + then (prn ("p:"^(to_string p)); + prn ("q:"^(to_string q)); + prn ("c:"^(to_string c)); + prn ("r:"^(to_string r)); prn ("d:"^(string_of_int d)); failwith "erreur dans la pseudo-division"); *) @@ -951,7 +945,7 @@ let pseudo_euclide coef_non_nuls p q x = let d = if d mod 2 = 1 then d+1 else d in (* on encore c^d * p = s*q + r, mais d pair *) - if eqP r cf0 + if equal r cf0 then ((*pr "reste nul"; *) (r,c,d,s,cf1,[],[])) else ( let r1 = vire_contenu r in @@ -967,7 +961,7 @@ let pseudo_euclide coef_non_nuls p q x = (try (while true do let rd = div_pol !r f x in (* verification de la division - if not (eqP cf0 ((!r)--(f@@rd))) + if not (equal cf0 ((!r)--(f@@rd))) then failwith "erreur dans la division"; *) k:=(!k)+1; @@ -1003,7 +997,7 @@ let pseudo_euclide coef_non_nuls p q x = !lf; (* pr ((* "------reste de même signe: " - ^(string_of_P c) + ^(to_string c) ^" variable: "^(string_of_int x) ^" c:"^(string_of_int (deg_total c)) ^" d:"^(string_of_int d) @@ -1026,7 +1020,7 @@ let pseudo_euclide coef_non_nuls p q x = divise un des polynômes de lp (dans Q[x1...xn]) *) let implique_non_nul lp p = - if eqP p cf0 then false + if equal p cf0 then false else( pr "["; let lf = facteurs2 p in @@ -1050,7 +1044,8 @@ let ajoute_non_nul p lp = then( let p = normc p in let lf = facteurs2 p in - let r = set_of_list_eq eqP (lp@lf@[p]) in + let r = set_of_list_eq equal (lp@lf@[p]) in r) else lp +end diff --git a/plugins/groebner/polynom.mli b/plugins/groebner/polynom.mli index 99c1701b8..79680262a 100644 --- a/plugins/groebner/polynom.mli +++ b/plugins/groebner/polynom.mli @@ -1,111 +1,120 @@ -type coef (*= Entiers.entiers*) -val coef_of_num : Num.num -> coef -val num_of_coef : coef -> Num.num -val eq_coef : coef -> coef -> bool -val lt_coef : coef -> coef -> bool -val le_coef : coef -> coef -> bool -val abs_coef : coef -> coef -val plus_coef : coef -> coef -> coef -val mult_coef : coef -> coef -> coef -val sub_coef : coef -> coef -> coef -val opp_coef : coef -> coef -val div_coef : coef -> coef -> coef -val mod_coef : coef -> coef -> coef -val coef0 : coef -val coef1 : coef -val string_of_coef : coef -> string -val int_of_coef : coef -> int -val hash_coef : coef -> int -val puis_coef : coef -> int -> coef -val pgcd : coef -> coef -> coef -val pgcd2 : coef -> coef -> coef +(************************************************************************) +(* 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 *) +(************************************************************************) -type variable = int -type poly = Pint of coef | Prec of variable * poly array +(* Building recursive polynom operations from a type of coefficients *) -val cf : int -> poly -val cf0 : poly -val cf1 : poly -val x : variable -> poly -val max_var_pol : poly -> variable -val max_var_pol2 : poly -> variable -val max_var : poly array -> variable -val eqP : poly -> poly -> bool -val norm : poly -> poly -val deg : variable -> poly -> int -val deg_total : poly -> int -val copyP : poly -> poly -val coef : variable -> int -> poly -> poly -val plusP : poly -> poly -> poly -val content : poly -> coef -val div_int : poly -> coef -> poly -val vire_contenu : poly -> poly -val vars : poly -> variable list -val int_of_Pint : poly -> coef -val multx : int -> variable -> poly -> poly -val multP : poly -> poly -> poly -val deriv : variable -> poly -> poly -val oppP : poly -> poly -val moinsP : poly -> poly -> poly -val puisP : poly -> int -> poly -val ( @@ ) : poly -> poly -> poly -val ( -- ) : poly -> poly -> poly -val ( ^^ ) : poly -> int -> poly -val coefDom : variable -> poly -> poly -val coefConst : variable -> poly -> poly -val remP : variable -> poly -> poly -val coef_int_tete : poly -> coef -val normc : poly -> poly -val is_constantP : poly -> bool -val is_zero : poly -> bool -val monome : variable -> int -> poly -val coef_constant : poly -> coef -val univ : bool ref -val string_of_var : int -> string -val nsP : int ref -val string_of_Pcut : poly -> string -val string_of_P : poly -> string -val printP : poly -> unit -val print_tpoly : poly array -> unit -val print_lpoly : poly list -> unit -val quo_rem_pol : poly -> poly -> variable -> poly * poly -val div_pol : poly -> poly -> variable -> poly -val divP : poly -> poly -> poly -val div_pol_rat : poly -> poly -> bool -val pseudo_div : poly -> poly -> variable -> poly * poly * int * poly -val pgcdP : poly -> poly -> poly -val pgcd_pol : poly -> poly -> variable -> poly -val content_pol : poly -> variable -> poly -val pgcd_coef_pol : poly -> poly -> variable -> poly -val pgcd_pol_rec : poly -> poly -> variable -> poly -val gcd_sub_res : poly -> poly -> variable -> poly -val gcd_sub_res_rec : poly -> poly -> poly -> poly -> int -> variable -> poly -val lazard_power : poly -> poly -> int -> variable -> poly -val sans_carre : poly -> variable -> poly -val facteurs : poly -> variable -> poly list -val facteurs_impairs : poly -> variable -> poly list -val hcontentP : (string, poly) Hashtbl.t -val prcontentP : unit -> unit -val contentP : poly * variable -> poly -val hashP : poly -> int -module Hashpol : Hashtbl.S with type key=poly -val memoP : - string -> 'a Hashpol.t -> (poly -> 'a) -> poly -> 'a -val hfactorise : poly list list Hashpol.t -val prfactorise : unit -> unit -val factorise : poly -> poly list list -val facteurs2 : poly -> poly list -val pol_de_factorisation : poly list list -> poly -val set_of_array_facteurs : poly list array -> poly list -val factorise_tableauP2 : - poly array -> - poly list array -> poly array * (poly * int list) array -val factorise_tableauP : - poly array -> poly array * (poly * int list) array -val is_positif : poly -> bool -val is_negatif : poly -> bool -val pseudo_euclide : - poly list -> poly -> poly -> variable -> - poly * poly * int * poly * poly * (poly * int) list * (poly * int) list -val implique_non_nul : poly list -> poly -> bool -val ajoute_non_nul : poly -> poly list -> poly list +module type Coef = sig + type t + val equal : t -> t -> bool + val lt : t -> t -> bool + val le : t -> t -> bool + val abs : t -> t + val plus : t -> t -> t + val mult : t -> t -> t + val sub : t -> t -> t + val opp : t -> t + val div : t -> t -> t + val modulo : t -> t -> t + val puis : t -> int -> t + val pgcd : t -> t -> t + + val hash : t -> int + val of_num : Num.num -> t + val to_string : t -> string +end + +module type S = sig + type coef + type variable = int + type t = Pint of coef | Prec of variable * t array + + val of_num : Num.num -> t + val x : variable -> t + val monome : variable -> int -> t + val is_constantP : t -> bool + val is_zero : t -> bool + + val max_var_pol : t -> variable + val max_var_pol2 : t -> variable + val max_var : t array -> variable + val equal : t -> t -> bool + val norm : t -> t + val deg : variable -> t -> int + val deg_total : t -> int + val copyP : t -> t + val coef : variable -> int -> t -> t + + val plusP : t -> t -> t + val content : t -> coef + val div_int : t -> coef -> t + val vire_contenu : t -> t + val vars : t -> variable list + val int_of_Pint : t -> coef + val multx : int -> variable -> t -> t + val multP : t -> t -> t + val deriv : variable -> t -> t + val oppP : t -> t + val moinsP : t -> t -> t + val puisP : t -> int -> t + val ( @@ ) : t -> t -> t + val ( -- ) : t -> t -> t + val ( ^^ ) : t -> int -> t + val coefDom : variable -> t -> t + val coefConst : variable -> t -> t + val remP : variable -> t -> t + val coef_int_tete : t -> coef + val normc : t -> t + val coef_constant : t -> coef + val univ : bool ref + val string_of_var : int -> string + val nsP : int ref + val to_string : t -> string + val printP : t -> unit + val print_tpoly : t array -> unit + val print_lpoly : t list -> unit + val quo_rem_pol : t -> t -> variable -> t * t + val div_pol : t -> t -> variable -> t + val divP : t -> t -> t + val div_pol_rat : t -> t -> bool + val pseudo_div : t -> t -> variable -> t * t * int * t + val pgcdP : t -> t -> t + val pgcd_pol : t -> t -> variable -> t + val content_pol : t -> variable -> t + val pgcd_coef_pol : t -> t -> variable -> t + val pgcd_pol_rec : t -> t -> variable -> t + val gcd_sub_res : t -> t -> variable -> t + val gcd_sub_res_rec : t -> t -> t -> t -> int -> variable -> t + val lazard_power : t -> t -> int -> variable -> t + val sans_carre : t -> variable -> t + val facteurs : t -> variable -> t list + val facteurs_impairs : t -> variable -> t list + val hcontentP : (string, t) Hashtbl.t + val prcontentP : unit -> unit + val contentP : t * variable -> t + val hash : t -> int + module Hashpol : Hashtbl.S with type key=t + val memoP : string -> 'a Hashpol.t -> (t -> 'a) -> t -> 'a + val hfactorise : t list list Hashpol.t + val prfactorise : unit -> unit + val factorise : t -> t list list + val facteurs2 : t -> t list + val pol_de_factorisation : t list list -> t + val set_of_array_facteurs : t list array -> t list + val factorise_tableauP2 : + t array -> t list array -> t array * (t * int list) array + val factorise_tableauP : t array -> t array * (t * int list) array + val is_positif : t -> bool + val is_negatif : t -> bool + val pseudo_euclide : + t list -> t -> t -> variable -> + t * t * int * t * t * (t * int) list * (t * int) list + val implique_non_nul : t list -> t -> bool + val ajoute_non_nul : t -> t list -> t list +end + +module Make (C:Coef) : S with type coef = C.t diff --git a/plugins/groebner/utile.ml b/plugins/groebner/utile.ml index 52a69dc10..fc7de1e33 100644 --- a/plugins/groebner/utile.ml +++ b/plugins/groebner/utile.ml @@ -16,7 +16,7 @@ let prt s = if !Flags.debug then (print_string (s^"\n");flush(stdout)) else () let info s = - output_string stderr s; flush stderr + Flags.if_verbose prerr_string s (********************************************************************** Listes diff --git a/plugins/setoid_ring/Field_tac.v b/plugins/setoid_ring/Field_tac.v index b31ebda55..0082eb9af 100644 --- a/plugins/setoid_ring/Field_tac.v +++ b/plugins/setoid_ring/Field_tac.v @@ -488,6 +488,10 @@ Ltac reduce_field_expr ope kont FLD fv expr := (* Hack to let a Ltac return a term in the context of a primitive tactic *) Ltac return_term x := generalize (refl_equal x). +Ltac get_term := + match goal with + | |- ?x = _ -> _ => x + end. (* Turn an operation on field expressions (FExpr) into a reduction on terms (in the field carrier). Because of field_lookup, diff --git a/tactics/tacinterp.ml b/tactics/tacinterp.ml index 6fe7b9288..a7b3a15e3 100644 --- a/tactics/tacinterp.ml +++ b/tactics/tacinterp.ml @@ -2077,8 +2077,6 @@ and interp_genarg ist gl x = in_gen (wit_tactic n) (TacArg(valueIn(VFun(ist.trace,ist.lfun,[], out_gen (globwit_tactic n) x)))) -(* in_gen (wit_tactic n) - (TacArg(valueIn(val_interp ist gl (out_gen (globwit_tactic n) x))))*) | None -> lookup_interp_genarg s ist gl x |