From ea60a60b66de43b0c51580587582ad82f21ebfb4 Mon Sep 17 00:00:00 2001 From: pottier Date: Thu, 3 Jun 2010 10:50:48 +0000 Subject: nsatz ajoute git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/coq/trunk@13058 85f007b7-540e-0410-9357-904b9bb8a0f7 --- plugins/nsatz/NsatzR.v | 406 +++++++++++++++ plugins/nsatz/NsatzZ.v | 73 +++ plugins/nsatz/ideal.ml | 1066 ++++++++++++++++++++++++++++++++++++++ plugins/nsatz/nsatz.ml4 | 608 ++++++++++++++++++++++ plugins/nsatz/nsatz_plugin.mllib | 5 + plugins/nsatz/polynom.ml | 679 ++++++++++++++++++++++++ plugins/nsatz/polynom.mli | 97 ++++ plugins/nsatz/utile.ml | 130 +++++ plugins/nsatz/utile.mli | 22 + plugins/nsatz/vo.itarget | 1 + 10 files changed, 3087 insertions(+) create mode 100644 plugins/nsatz/NsatzR.v create mode 100644 plugins/nsatz/NsatzZ.v create mode 100644 plugins/nsatz/ideal.ml create mode 100644 plugins/nsatz/nsatz.ml4 create mode 100644 plugins/nsatz/nsatz_plugin.mllib create mode 100644 plugins/nsatz/polynom.ml create mode 100644 plugins/nsatz/polynom.mli create mode 100644 plugins/nsatz/utile.ml create mode 100644 plugins/nsatz/utile.mli create mode 100644 plugins/nsatz/vo.itarget (limited to 'plugins') diff --git a/plugins/nsatz/NsatzR.v b/plugins/nsatz/NsatzR.v new file mode 100644 index 000000000..0664a3f78 --- /dev/null +++ b/plugins/nsatz/NsatzR.v @@ -0,0 +1,406 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* x = y. +intros x y H; replace x with ((x - y) + y); + [rewrite H | idtac]; ring. +Qed. + +Lemma psos_r1: forall x y, x = y -> x - y = 0. +intros x y H; rewrite H; ring. +Qed. + +Lemma nsatzR_not1: forall x y:R, x<>y -> exists z:R, z*(x-y)-1=0. +intros. +exists (1/(x-y)). +field. +unfold not. +unfold not in H. +intros. +apply H. +replace x with ((x-y)+y). +rewrite H0. +ring. +ring. +Qed. + +Lemma nsatzR_not1_0: forall x:R, x<>0 -> exists z:R, z*x-1=0. +intros. +exists (1/(x)). +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 nsatzR_not2: 1<>0. +auto with *. +Qed. + +Lemma nsatzR_diff: forall x y:R, x<>y -> x-y<>0. +intros. +intro; apply H. +replace x with (x-y+y) by ring. +rewrite H0; ring. +Qed. + +(* Removes x<>0 from hypothesis *) +Ltac nsatzR_not_hyp:= + match goal with + | H: ?x<>?y |- _ => + match y with + |0 => + let H1:=fresh "Hnsatz" in + let y:=fresh "x" in + destruct (@nsatzR_not1_0 _ H) as (y,H1); clear H + |_ => generalize (@nsatzR_diff _ _ H); clear H; intro + end + end. + +Ltac nsatzR_not_goal := + match goal with + | |- ?x<>?y :> R => red; intro; apply nsatzR_not2 + | |- False => apply nsatzR_not2 + end. + +Ltac nsatzR_begin := + intros; + repeat nsatzR_not_hyp; + try nsatzR_not_goal; + try apply psos_r1b; + repeat equalities_to_goal. + +(* code de Benjamin *) + +Definition PolZ := Pol Z. +Definition PEZ := PExpr Z. + +Definition P0Z : PolZ := @P0 Z 0%Z. + +Definition PolZadd : PolZ -> PolZ -> PolZ := + @Padd Z 0%Z Zplus Zeq_bool. + +Definition PolZmul : PolZ -> PolZ -> PolZ := + @Pmul Z 0%Z 1%Z Zplus Zmult Zeq_bool. + +Definition PolZeq := @Peq Z Zeq_bool. + +Definition norm := + @norm_aux Z 0%Z 1%Z Zplus Zmult Zminus Zopp Zeq_bool. + +Fixpoint mult_l (la : list PEZ) (lp: list PolZ) : PolZ := + match la, lp with + | a::la, p::lp => PolZadd (PolZmul (norm a) p) (mult_l la lp) + | _, _ => P0Z + end. + +Fixpoint compute_list (lla: list (list PEZ)) (lp:list PolZ) := + match lla with + | List.nil => lp + | la::lla => compute_list lla ((mult_l la lp)::lp) + end. + +Definition check (lpe:list PEZ) (qe:PEZ) (certif: list (list PEZ) * list PEZ) := + let (lla, lq) := certif in + let lp := List.map norm lpe in + PolZeq (norm qe) (mult_l lq (compute_list lla lp)). + + +(* Correction *) +Definition PhiR : list R -> PolZ -> R := + (Pphi 0 Rplus Rmult (gen_phiZ 0 1 Rplus Rmult Ropp)). + +Definition PEevalR : list R -> PEZ -> R := + PEeval 0 Rplus Rmult Rminus Ropp (gen_phiZ 0 1 Rplus Rmult Ropp) + Nnat.nat_of_N pow. + +Lemma P0Z_correct : forall l, PhiR l P0Z = 0. +Proof. trivial. Qed. + + +Lemma PolZadd_correct : forall P' P l, + PhiR l (PolZadd P P') = (PhiR l P + PhiR l P'). +Proof. + refine (Padd_ok Rset Rext (Rth_ARth Rset Rext (F_R Rfield)) + (gen_phiZ_morph Rset Rext (F_R Rfield))). +Qed. + +Lemma PolZmul_correct : forall P P' l, + PhiR l (PolZmul P P') = (PhiR l P * PhiR l P'). +Proof. + refine (Pmul_ok Rset Rext (Rth_ARth Rset Rext (F_R Rfield)) + (gen_phiZ_morph Rset Rext (F_R Rfield))). +Qed. + +Lemma norm_correct : + forall (l : list R) (pe : PEZ), PEevalR l pe = PhiR l (norm pe). +Proof. + intros;apply (norm_aux_spec Rset Rext (Rth_ARth Rset Rext (F_R Rfield)) + (gen_phiZ_morph Rset Rext (F_R Rfield)) R_power_theory) with (lmp:= List.nil). + compute;trivial. +Qed. + +Lemma PolZeq_correct : forall P P' l, + PolZeq P P' = true -> + PhiR l P = PhiR l P'. +Proof. + intros;apply + (Peq_ok Rset Rext (gen_phiZ_morph Rset Rext (F_R Rfield)));trivial. +Qed. + +Fixpoint Cond0 (A:Type) (Interp:A->R) (l:list A) : Prop := + match l with + | List.nil => True + | a::l => Interp a = 0 /\ Cond0 A Interp l + end. + +Lemma mult_l_correct : forall l la lp, + Cond0 PolZ (PhiR l) lp -> + PhiR l (mult_l la lp) = 0. +Proof. + induction la;simpl;intros;trivial. + destruct lp;trivial. + simpl in H;destruct H. + rewrite PolZadd_correct, PolZmul_correct, H, IHla;[ring | trivial]. +Qed. + +Lemma compute_list_correct : forall l lla lp, + Cond0 PolZ (PhiR l) lp -> + Cond0 PolZ (PhiR l) (compute_list lla lp). +Proof. + induction lla;simpl;intros;trivial. + apply IHlla;simpl;split;trivial. + apply mult_l_correct;trivial. +Qed. + +Lemma check_correct : + forall l lpe qe certif, + check lpe qe certif = true -> + Cond0 PEZ (PEevalR l) lpe -> + PEevalR l qe = 0. +Proof. + unfold check;intros l lpe qe (lla, lq) H2 H1. + apply PolZeq_correct with (l:=l) in H2. + rewrite norm_correct, H2. + apply mult_l_correct. + apply compute_list_correct. + clear H2 lq lla qe;induction lpe;simpl;trivial. + simpl in H1;destruct H1. + rewrite <- norm_correct;auto. +Qed. + +(* fin du code de Benjamin *) + +Lemma nsatzR_l3:forall c p r, ~c=0 -> c*p^r=0 -> p=0. +intros. +elim (Rmult_integral _ _ H0);intros. + 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 + |h : (?p = ?q)|- _ => revert h + end). + +Ltac lpol_goal t := + match t with + | ?a = 0 -> ?b => + let r:= lpol_goal b in + constr:(a::r) + | ?a = 0 => constr:(a::nil) + end. + +Fixpoint IPR p {struct p}: R := + match p with + xH => 1 + | xO xH => 1 + 1 + | xO p1 => 2 * (IPR p1) + | xI xH => 1 + (1 + 1) + | xI p1 => 1 + 2 * (IPR p1) + end. + +Definition IZR1 z := + match z with Z0 => 0 + | Zpos p => IPR p + | Zneg p => -(IPR p) + end. + +Fixpoint interpret3 t fv {struct t}: R := + match t with + | (PEadd t1 t2) => + let v1 := interpret3 t1 fv in + let v2 := interpret3 t2 fv in (v1 + v2) + | (PEmul t1 t2) => + let v1 := interpret3 t1 fv in + let v2 := interpret3 t2 fv in (v1 * v2) + | (PEsub t1 t2) => + let v1 := interpret3 t1 fv in + let v2 := interpret3 t2 fv in (v1 - v2) + | (PEopp t1) => + let v1 := interpret3 t1 fv in (-v1) + | (PEpow t1 t2) => + let v1 := interpret3 t1 fv in v1 ^(Nnat.nat_of_N t2) + | (PEc t1) => (IZR1 t1) + | (PEX n) => List.nth (pred (nat_of_P n)) fv 0 + end. + +(* lp est incluse dans fv. La met en tete. *) + +Ltac parametres_en_tete fv lp := + match fv with + | (@nil _) => lp + | (@cons _ ?x ?fv1) => + let res := AddFvTail x lp in + parametres_en_tete fv1 res + end. + +Ltac append1 a l := + match l with + | (@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 + | (cons ?x ?l) => let l' := rev l in append1 x l' + end. + + +Ltac nsatz_call_n info nparam p rr lp kont := + nsatz_compute (PEc info :: 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 nsatz_call radicalmax info nparam p lp kont := + let rec try_n n := + lazymatch n with + | 0%N => fail + | _ => +(* idtac "Trying power: " n;*) + (let r := eval compute in (Nminus radicalmax (Npred n)) in + nsatz_call_n info nparam p r lp kont) || + let n' := eval compute in (Npred n) in try_n n' + end in + try_n radicalmax. + + +Ltac nsatzR_gen radicalmax info lparam lvar n RNG lH _rl := + get_Pre RNG (); + 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 + intros; + let fv := + match lvar with + | nil => + let fv1 := FV_hypo_tac mkFV ltac:(get_Eq RNG) lH in + let fv1 := list_fold_right mkFV fv1 lpol in + rev fv1 + (* heuristique: les dernieres variables auront le poid le plus fort *) + | _ => lvar + end in + check_fv fv; + (*idtac "variables:";idtac fv;*) + let nparam := eval compute in (Z_of_nat (List.length lparam)) in + let fv := parametres_en_tete fv lparam in + (* idtac "variables:"; idtac fv; + idtac "nparam:"; idtac nparam;*) + let lpol := list_fold_right + ltac:(fun p l => let p' := mkPol p fv in constr:(p'::l)) + (@nil (PExpr Z)) + lpol in + let lpol := eval compute in (List.rev lpol) in + (*idtac lpol;*) + let SplitPolyList kont := + match lpol with + | ?p2::?lp2 => kont p2 lp2 + end in + SplitPolyList ltac:(fun p lp => + set (p21:=p) ; + set (lp21:=lp); + nsatz_call radicalmax info 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 nsatz certificate" + | let Hg2 := fresh "Hg" in + assert (Hg2: interpret3 q fv = 0); + [ simpl; apply (@check_correct fv lp21 q (lci,lq) Hg); simpl; + repeat (split;[assumption|idtac]); exact I + | simpl in Hg2; simpl; + apply nsatzR_l3 with (interpret3 c fv) (Nnat.nat_of_N r);simpl; + [ discrR || idtac "could not prove discrimination result" + | exact Hg2] + ] + ])). + +Ltac nsatzRpv radicalmax info lparam lvar:= + nsatzR_begin; + intros; + let G := Get_goal in + ring_lookup + (PackRing ltac:(nsatzR_gen radicalmax info lparam lvar ring_subst_niter)) + [] G. + +Ltac nsatzR := nsatzRpv 6%N 1%Z (@nil R) (@nil R). +Ltac nsatzRradical radicalmax := nsatzRpv radicalmax 1%Z (@nil R) (@nil R). +Ltac nsatzRparameters lparam := nsatzRpv 6%N 1%Z lparam (@nil R). + +Tactic Notation "nsatz" := nsatzR. +Tactic Notation "nsatz" "with" "lexico" := + nsatzRpv 6%N 2%Z (@nil R) (@nil R). +Tactic Notation "nsatz" "with" "lexico" "sugar":= + nsatzRpv 6%N 3%Z (@nil R) (@nil R). +Tactic Notation "nsatz" "without" "sugar":= + nsatzRpv 6%N 0%Z (@nil R) (@nil R). + + diff --git a/plugins/nsatz/NsatzZ.v b/plugins/nsatz/NsatzZ.v new file mode 100644 index 000000000..a65efac27 --- /dev/null +++ b/plugins/nsatz/NsatzZ.v @@ -0,0 +1,73 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* IZR x = IZR y. +Proof IZR_eq. (* or f_equal ... *) + +Lemma nsatzZconclR: forall x y:Z, IZR x = IZR y -> x = y. +Proof eq_IZR. + +Lemma nsatzZhypnotR: forall x y:Z, x<>y -> IZR x <> IZR y. +Proof IZR_neq. + +Lemma nsatzZconclnotR: forall x y:Z, IZR x <> IZR y -> x <> y. +Proof. +intros x y H. contradict H. f_equal. assumption. +Qed. + +Ltac nsatzZtoR1 := + repeat + (match goal with + | H:(@eq Z ?x ?y) |- _ => + generalize (@nsatzZhypR _ _ H); clear H; intro H + | |- (@eq Z ?x ?y) => apply nsatzZconclR + | H:not (@eq Z ?x ?y) |- _ => + generalize (@nsatzZhypnotR _ _ H); clear H; intro H + | |- not (@eq Z ?x ?y) => apply nsatzZconclnotR + end). + +Lemma nsatzZR1: forall x y:Z, IZR(x+y) = (IZR x + IZR y)%R. +Proof plus_IZR. + +Lemma nsatzZR2: forall x y:Z, IZR(x*y) = (IZR x * IZR y)%R. +Proof mult_IZR. + +Lemma nsatzZR3: forall x y:Z, IZR(x-y) = (IZR x - IZR y)%R. +Proof. +intros; symmetry. apply Z_R_minus. +Qed. + +Lemma nsatzZR4: forall (x:Z) p, IZR(x ^ Zpos p) = (IZR x ^ nat_of_P p)%R. +Proof. +intros. rewrite pow_IZR. +do 2 f_equal. +apply Zpos_eq_Z_of_nat_o_nat_of_P. +Qed. + +Ltac nsatzZtoR2:= + repeat + (rewrite nsatzZR1 in * || + rewrite nsatzZR2 in * || + rewrite nsatzZR3 in * || + rewrite nsatzZR4 in *). + +Ltac nsatzZ_begin := + intros; + nsatzZtoR1; + nsatzZtoR2; + simpl in *. + (*cbv beta iota zeta delta [nat_of_P Pmult_nat plus mult] in *.*) + +Ltac nsatzZ := + nsatzZ_begin; (*idtac "nsatzZ_begin;";*) + nsatzR. diff --git a/plugins/nsatz/ideal.ml b/plugins/nsatz/ideal.ml new file mode 100644 index 000000000..b4a0e7ade --- /dev/null +++ b/plugins/nsatz/ideal.ml @@ -0,0 +1,1066 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* mon -> mon +val deg : mon -> int +val compare_mon : mon -> mon -> int +val div_mon : mon -> mon -> mon +val div_mon_test : mon -> mon -> bool +val ppcm_mon : mon -> mon -> mon + +(* Polynomials *) + +type deg = int +type coef +type poly +type polynom + +val repr : poly -> (coef * mon) list +val polconst : coef -> poly +val zeroP : poly +val gen : 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 : poly -> poly -> poly +val mult_t_pol : coef -> mon -> poly -> poly +val selectdiv : mon -> poly list -> poly +val oppP : poly -> poly +val emultP : coef -> poly -> poly +val multP : poly -> poly -> poly +val puisP : poly -> int -> poly +val contentP : poly -> coef +val contentPlist : poly list -> coef +val pgcdpos : coef -> coef -> coef +val div_pol : poly -> poly -> coef -> coef -> mon -> poly +val reduce2 : 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 : poly list -> unit +val reduce2_trace : poly -> poly list -> poly list -> poly list * poly +val spol : poly -> poly -> poly +val etrangers : poly -> poly -> bool +val div_ppcm : poly -> poly -> poly -> bool + +val genPcPf : poly -> poly list -> poly list -> poly list +val genOCPf : 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 : 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 + +(* division of tail monomials *) + +let reduire_les_queues = false + +(* divide first with new polynomials *) + +let nouveaux_pol_en_tete = 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)^"]" + +(*********************************************************************** + Monomials + array of integers, first is the degree +*) + +type mon = int array +type deg = int +type poly = (coef * mon) list +type polynom = + {pol : poly ref; + num : int; + sugar : int} + +let nvar m = Array.length m - 1 + +let deg m = m.(0) + +let mult_mon m m' = + let d = nvar m in + let m'' = Array.create (d+1) 0 in + for i=0 to d do + m''.(i)<- (m.(i)+m'.(i)); + done; + m'' + + +let compare_mon m m' = + let d = nvar m in + 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) + +let div_mon m m' = + let d = nvar m in + 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' divides m *) +let div_mon_test m m' = + let d = nvar m in + let res=ref true in + let i=ref 0 in (*il faut que le degre total soit bien mis sinon, i=ref 1*) + while (!res) && (!i<=d) do + res:= (m.(!i) >= m'.(!i)); + i:=succ !i; + done; + !res + +let set_deg m = + let d = nvar m in + m.(0)<-0; + for i=1 to d do + m.(0)<- m.(i)+m.(0); + done; + m + +(* lcm *) +let ppcm_mon m m' = + let d = nvar m in + let m'' = Array.create (d+1) 0 in + for i=1 to d do + m''.(i)<- (max m.(i) m'.(i)); + done; + set_deg m'' + + + +(********************************************************************** + Polynomials + list of (coefficient, monomial) decreasing order +*) + +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 = map fst p in + let m = map snd p in + 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 (nth lv i) + with _ -> (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 p = + 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 + p + +let nsP2 = ref max_int + +let stringPcut p = + (*Polynomesrec.nsP1:=20;*) + nsP2:=10; + let res = + if (length p)> !nsP2 + then (stringP [hd p])^" + "^(string_of_int (length p))^" termes" + else stringP p in + (*Polynomesrec.nsP1:= max_int;*) + nsP2:= max_int; + res + +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 = [] + +(* returns a constant polynom ial with d variables *) +let polconst d c = + let m = Array.create (d+1) 0 in + let m = set_deg m in + [(c,m)] + +let plusP p q = + let rec plusP p q = + match p with + [] -> q + |t::p' -> + match q with + [] -> p + |t'::q' -> + match compare_mon (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 by (a,monomial) *) +let mult_t_pol a m p = + let rec mult_t_pol p = + match p with + [] -> [] + |(b,m')::p -> ((P.multP a b),(mult_mon m m'))::(mult_t_pol p) + in mult_t_pol p + +let coef_of_int x = P.of_num (Num.Int x) + +(* variable i *) +let gen d i = + let m = Array.create (d+1) 0 in + m.(i) <- 1; + let m = set_deg m in + [((coef_of_int 1),m)] + +let oppP p = + let rec oppP p = + match p with + [] -> [] + |(b,m')::p -> ((P.oppP b),m')::(oppP p) + in oppP p + +(* multiplication by a coefficient *) +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 p q = + let rec aux p = + match p with + [] -> [] + |(a,m)::p' -> plusP (mult_t_pol a m q) (aux p') + in aux p + +let puisP p n= + match p with + [] -> [] + |_ -> + let d = nvar (snd (hd p)) in + let rec puisP n = + match n with + 0 -> [coef1, Array.create (d+1) 0] + | 1 -> p + |_ -> multP p (puisP (n-1)) + in puisP n + +let rec contentP p = + match p with + |[] -> coef1 + |[a,m] -> a + |(a,m)::p1 -> + if P.equal a coef1 || P.equal a coefm1 + then a + else P.pgcdP a (contentP p1) + +let contentPlist lp = + match lp with + |[] -> coef1 + |p::l1 -> + fold_left + (fun r q -> + if P.equal r coef1 || P.equal r coefm1 + then r + else P.pgcdP r (contentP q)) + (contentP p) l1 + +(*********************************************************************** + Division of polynomials + *) + +let pgcdpos a b = P.pgcdP a b + +let polynom0 = {pol = ref []; num = 0; sugar = 0} + +let ppol p = !(p.pol) + +let lm p = snd (hd (ppol p)) + +let nallpol = ref 0 + +let allpol = ref (Array.create 1000 polynom0) + +let new_allpol p s = + nallpol := !nallpol + 1; + if !nallpol >= Array.length !allpol + then + allpol := Array.append !allpol (Array.create !nallpol polynom0); + let p = {pol = ref p; num = !nallpol; sugar = s} in + !allpol.(!nallpol)<- p; + p + +(* returns a polynomial of l whose head monomial divides m, else [] *) + +let rec selectdiv m l = + match l with + [] -> polynom0 + |q::r -> let m'= snd (hd (ppol q)) in + match (div_mon_test m m') with + true -> q + |false -> selectdiv m r + +let div_pol p q a b m = +(* info ".";*) + plusP (emultP a p) (mult_t_pol b m q) + +let hmon = Hashtbl.create 1000 + +let use_hmon = ref false + +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 + else () + +let div_coef a b = P.divP a b + + +(* remainder r of the division of p by polynomials of l, returns (c,r) where c is the coefficient for pseudo-division : c p = sum_i q_i p_i + r *) + +let reduce2 p l = + let l = if nouveaux_pol_en_tete then rev l else l in + let rec reduce p = + match p with + [] -> (coef1,[]) + |t::p' -> + let (a,m)=t in + let q = (try find_hmon m + with Not_found -> + let q = selectdiv m l in + match (ppol q) with + t'::q' -> (add_hmon m q; + q) + |[] -> q) in + match (ppol 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'=(P.oppP (div_coef a c)) in + let (e,r)=reduce (div_pol p' q' a' b' + (div_mon m m')) in + (P.multP a' e,r) + in let (c,r) = reduce p in + (c,r) + +(* trace of divisions *) + +(* list of initial polynomials *) +let poldep = ref [] +let poldepcontent = ref [] + +(* coefficients of polynomials when written with initial polynomials *) +let coefpoldep = Hashtbl.create 51 + +(* coef of q in p = sum_i c_i*q_i *) +let coefpoldep_find p q = + try (Hashtbl.find coefpoldep (p.num,q.num)) + with _ -> [] + +let coefpoldep_remove p q = + Hashtbl.remove coefpoldep (p.num,q.num) + +let coefpoldep_set p q c = + Hashtbl.add coefpoldep (p.num,q.num) c + +let initcoefpoldep d lp = + poldep:=lp; + poldepcontent:= map (fun p -> contentP (ppol p)) lp; + iter + (fun p -> coefpoldep_set p p (polconst d (coef_of_int 1))) + lp + +(* keeps trace in coefpoldep + divides without pseudodivisions *) + +let reduce2_trace p l lcp = + let l = if nouveaux_pol_en_tete then rev l else l in + (* 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 find_hmon m + with Not_found -> + let q = selectdiv m l in + match (ppol q) with + t'::q' -> (add_hmon m q; + q) + |[] -> q) in + match (ppol 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 (div_coef a b)) in + let m''= div_mon m m' in + let p1=plusP p' (mult_t_pol 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"; + iter + (fun (a,m,s) -> + let x = mult_t_pol a m s in + info ((stringP x)^"\n")) + lq; + info "ok\n";*) + (map2 + (fun c0 q -> + let c = + fold_left + (fun x (a,m,s) -> + if equal (ppol s) (ppol q) + then + plusP x (mult_t_pol a m (polconst (nvar m) (coef_of_int 1))) + else x) + c0 + lq in + c) + lcp + !poldep, + r) + +let homogeneous = ref false +let pol_courant = ref polynom0 + +(*********************************************************************** + Completion + *) + +let sugar_flag = ref true + +let compute_sugar p = + fold_left (fun s (a,m) -> max s m.(0)) 0 p + +let mk_polynom p = + new_allpol p (compute_sugar p) + +let spol ps qs= + let p = ppol ps in + let q = ppol qs in + let m = snd (hd p) in + let m'= snd (hd q) in + let a = fst (hd p) in + let b = fst (hd q) in + let p'= tl p in + let q'= tl q in + let c = (pgcdpos a b) in + let m''=(ppcm_mon m m') in + let m1 = div_mon m'' m in + let m2 = div_mon m'' m' in + let fsp p' q' = + plusP + (mult_t_pol + (div_coef b c) + m1 p') + (mult_t_pol + (P.oppP (div_coef a c)) + m2 q') in + let sp = fsp p' q' in + let sps = + new_allpol + sp + (max (m1.(0) + ps.sugar) (m2.(0) + qs.sugar)) in + coefpoldep_set sps ps (fsp (polconst (nvar m) (coef_of_int 1)) []); + coefpoldep_set sps qs (fsp [] (polconst (nvar m) (coef_of_int 1))); + sps + + +let etrangers p p'= + let m = snd (hd p) in + let m'= snd (hd p') in + let d = nvar m 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 if head monomial of p'' divides lcm of lhead monomials of p and p' *) + +let div_ppcm p p' p'' = + let m = snd (hd p) in + let m'= snd (hd p') in + let m''= snd (hd p'') in + let d = nvar m 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 from extraction of Laurent Théry Coq program *) + +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 i a q = + list_rec + (match etrangers (ppol i) (ppol a) with + true -> DontKeep [] + | false -> Keep []) + (fun b q1 rec_ren -> + match div_ppcm (ppol i) (ppol a) (ppol b) with + true -> DontKeep (b::q1) + | false -> + (match div_ppcm (ppol i) (ppol b) (ppol a) with + true -> rec_ren + | false -> addRes b rec_ren)) q + +(* sugar strategy *) + +let rec addS x l = l @ [x] (* oblige de mettre en queue sinon le certificat deconne *) + +let addSsugar x l = + if !sugar_flag + then + let sx = x.sugar in + let rec insere l = + match l with + | [] -> [x] + | y::l1 -> + if sx <= y.sugar + 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 i aP q = + (let rec genPc aP0 = + match aP0 with + [] -> (fun r -> r) + | a::l1 -> + (fun l -> + (match slice i a l1 with + Keep l2 -> addSsugar (spol i a) (genPc l2 l) + | DontKeep l2 -> genPc l2 l)) + in genPc aP) q + +let genOCPf h' = + list_rec [] (fun a l rec_ren -> + genPcPf a l rec_ren) h' + +(*********************************************************************** + critical pairs/s-polynomials + *) + +let ordcpair ((i1,j1),m1) ((i2,j2),m2) = +(* let s1 = (max + (!allpol.(i1).sugar + m1.(0) + - (snd (hd (ppol !allpol.(i1)))).(0)) + (!allpol.(j1).sugar + m1.(0) + - (snd (hd (ppol !allpol.(j1)))).(0))) in + let s2 = (max + (!allpol.(i2).sugar + m2.(0) + - (snd (hd (ppol !allpol.(i2)))).(0)) + (!allpol.(j2).sugar + m2.(0) + - (snd (hd (ppol !allpol.(j2)))).(0))) in + match compare s1 s2 with + | 1 -> 1 + |(-1) -> -1 + |0 -> compare_mon m1 m2*) + + compare_mon m1 m2 + +let sortcpairs lcp = + sort ordcpair lcp + +let mergecpairs l1 l2 = + merge ordcpair l1 l2 + +let ord i j = + if i r @ (cpair p q)) [] lq) + +let cpairs lp = + let rec aux l = + match l with + []|[_] -> [] + |p::l1 -> mergecpairs (cpairs1 p l1) (aux l1) + in aux lp + + +let critere2 ((i,j),m) lp lcp = + exists + (fun h -> + h.num <> i && h.num <> j + && (div_mon_test m (lm h)) + && (let c1 = ord i h.num in + not (exists (fun (c,_) -> c1 = c) lcp)) + && (let c1 = ord j h.num in + not (exists (fun (c,_) -> c1 = c) lcp))) + lp + +let critere3 ((i,j),m) lp lcp = + exists + (fun h -> + h.num <> i && h.num <> j + && (div_mon_test m (lm h)) + && (h.num < j + || not (m = ppcm_mon + (lm (!allpol.(i))) + (lm h))) + && (h.num < i + || not (m = ppcm_mon + (lm (!allpol.(j))) + (lm h)))) + lp + +let add_cpairs p lp lcp = + mergecpairs (cpairs1 p lp) lcp + +let prod_rec f = function + (x0, x) -> f x0 x + +let step = ref 0 + +let infobuch p q = + if !step = 0 + then (info ("[" ^ (string_of_int (length p)) + ^ "," ^ (string_of_int (length q)) + ^ "]")) + +(* in lp new polynomials are at the end *) + +let coef_courant = ref coef1 + +type certificate = + { coef : coef; power : int; + gb_comb : poly list list; last_comb : poly list } + +let test_dans_ideal p lp lp0 = + let (c,r) = reduce2 (ppol !pol_courant) lp in + info ("remainder: "^(stringPcut r)^"\n"); + coef_courant:= P.multP !coef_courant c; + pol_courant:= mk_polynom r; + if r=[] + then (info "polynomial reduced to 0\n"; + let lcp = map (fun q -> []) !poldep in + let c = !coef_courant in + let (lcq,r) = reduce2_trace (emultP c p) lp lcp in + info "r ok\n"; + info ("r: "^(stringP r)^"\n"); + let res=ref (emultP c p) in + iter2 + (fun cq q -> res:=plusP (!res) (multP cq (ppol q)); + ) + lcq !poldep; + info ("verif sum: "^(stringP (!res))^"\n"); + info ("coefficient: "^(stringP (polconst 1 c))^"\n"); + let rec aux lp = + match lp with + |[] -> [] + |p::lp -> + (map + (fun q -> coefpoldep_find p q) + lp)::(aux lp) + in + let coefficient_multiplicateur = c in + let liste_polynomes_de_depart = rev lp0 in + let polynome_a_tester = p in + let liste_des_coefficients_intermediaires = + (let lci = rev (aux (rev lp)) in + let lci = ref lci (* (map rev lci) *) in + iter (fun x -> lci := tl (!lci)) lp0; + !lci) in + let liste_des_coefficients = + map + (fun cq -> emultP (coef_of_int (-1)) cq) + (rev lcq) in + (liste_polynomes_de_depart, + polynome_a_tester, + {coef = coefficient_multiplicateur; + power = 1; + gb_comb = liste_des_coefficients_intermediaires; + last_comb = liste_des_coefficients}) + ) + else ((*info "polynomial not reduced to 0\n"; + info ("\nremainder: "^(stringPcut r)^"\n");*) + raise NotInIdeal) + +let divide_rem_with_critical_pair = ref false + +let list_diff l x = + filter (fun y -> y <> x) l + +let choix_spol p l = + l + +let deg_hom p = + match p with + | [] -> -1 + | (a,m)::_ -> m.(0) + +let pbuchf pq p lp0= + info "computation of the Groebner basis\n"; + step:=0; + Hashtbl.clear hmon; + let rec pbuchf x = + prod_rec (fun lp lpc -> + infobuch lp lpc; +(* step:=(!step+1)mod 10;*) + match lpc with + [] -> + + (* info ("List of polynomials:\n"^(fold_left (fun r p -> r^(stringP p)^"\n") "" lp)); + info "--------------------\n";*) + test_dans_ideal (ppol p) lp lp0 + | _ -> + let (((i,j),m)::lpc2)= choix_spol !pol_courant lpc in +(* info "choosen pair\n";*) + if critere3 ((i,j),m) lp lpc2 + then (info "c"; pbuchf (lp, lpc2)) + else + let a = spol !allpol.(i) !allpol.(j) in + if !homogeneous && (ppol a)<>[] && deg_hom (ppol a) + > deg_hom (ppol !pol_courant) + then (info "h"; pbuchf (lp, lpc2)) + else +(* let sa = a.sugar in*) + let (ca,a0)= reduce2 (ppol a) lp in + match a0 with + [] -> info "0";pbuchf (lp, lpc2) + | _ -> +(* info "pair reduced\n";*) + a.pol := emultP ca (ppol a); + let (lca,a0) = reduce2_trace (ppol a) lp + (map (fun q -> emultP ca (coefpoldep_find a q)) + !poldep) in +(* info "paire re-reduced";*) + a.pol := a0; +(* let a0 = new_allpol a0 sa in*) + iter2 (fun c q -> + coefpoldep_remove a q; + coefpoldep_set a q c) lca !poldep; + let a0 = a in + info ("\nnew polynomials: "^(stringPcut (ppol a0))^"\n"); + let ct = coef1 (* contentP a0 *) in + (*info ("content: "^(string_of_coef ct)^"\n");*) + poldep:=addS a0 lp; + poldepcontent:=addS ct (!poldepcontent); + + try test_dans_ideal (ppol p) (addS a0 lp) lp0 + with NotInIdeal -> + let newlpc = add_cpairs a0 lp lpc2 in + pbuchf (((addS a0 lp), newlpc))) + x + in pbuchf pq + +let is_homogeneous p = + match p with + | [] -> true + | (a,m)::p1 -> let d = m.(0) in + for_all (fun (b,m') -> m'.(0)=d) p1 + +(* returns + 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] + + such that + c*p = sum qi*pi + where pn+k = a(n+k,n+k-1)*pn+k-1 + ... + a(n+k,1)* p1 + *) + +let in_ideal d lp p = + Hashtbl.clear hmon; + Hashtbl.clear coefpoldep; + nallpol := 0; + allpol := Array.create 1000 polynom0; + homogeneous := for_all is_homogeneous (p::lp); + if !homogeneous then info "homogeneous polynomials\n"; + info ("p: "^(stringPcut p)^"\n"); + info ("lp:\n"^(fold_left (fun r p -> r^(stringPcut p)^"\n") "" lp)); + (*info ("p: "^(stringP p)^"\n"); + info ("lp:\n"^(fold_left (fun r p -> r^(stringP p)^"\n") "" lp));*) + + let lp = map mk_polynom lp in + let p = mk_polynom p in + initcoefpoldep d lp; + coef_courant:=coef1; + pol_courant:=p; + + let (lp1,p1,cert) = + try test_dans_ideal (ppol p) lp lp + with NotInIdeal -> pbuchf (lp, (cpairs lp)) p lp in + info "computed\n"; + + (map ppol lp1, p1, cert) + +(* *) +end + + + diff --git a/plugins/nsatz/nsatz.ml4 b/plugins/nsatz/nsatz.ml4 new file mode 100644 index 000000000..892d60373 --- /dev/null +++ b/plugins/nsatz/nsatz.ml4 @@ -0,0 +1,608 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* 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 +*) + +(* ------------------------------------------------------------------------- *) +(* ------------------------------------------------------------------------- *) + +type vname = string + +type term = + | Zero + | 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 + +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 + +let tpexpr = + lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PExpr") +let ttconst = lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PEc") +let ttvar = lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PEX") +let ttadd = lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PEadd") +let ttsub = lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PEsub") +let ttmul = lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PEmul") +let ttopp = lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PEopp") +let ttpow = lazy (gen_constant "CC" ["setoid_ring";"Ring_polynom"] "PEpow") + +let tlist = lazy (gen_constant "CC" ["Lists";"List"] "list") +let lnil = lazy (gen_constant "CC" ["Lists";"List"] "nil") +let lcons = lazy (gen_constant "CC" ["Lists";"List"] "cons") + +let tz = lazy (gen_constant "CC" ["ZArith";"BinInt"] "Z") +let z0 = lazy (gen_constant "CC" ["ZArith";"BinInt"] "Z0") +let zpos = lazy (gen_constant "CC" ["ZArith";"BinInt"] "Zpos") +let zneg = lazy(gen_constant "CC" ["ZArith";"BinInt"] "Zneg") + +let pxI = lazy(gen_constant "CC" ["NArith";"BinPos"] "xI") +let pxO = lazy(gen_constant "CC" ["NArith";"BinPos"] "xO") +let pxH = lazy(gen_constant "CC" ["NArith";"BinPos"] "xH") + +let nN0 = lazy (gen_constant "CC" ["NArith";"BinNat"] "N0") +let nNpos = lazy(gen_constant "CC" ["NArith";"BinNat"] "Npos") + +let mkt_app name l = mkApp (Lazy.force name, Array.of_list l) + +let tlp () = mkt_app tlist [mkt_app tpexpr [Lazy.force tz]] +let tllp () = mkt_app tlist [tlp()] + +let rec mkt_pos n = + if n =/ num_1 then Lazy.force pxH + else if mod_num n num_2 =/ num_0 then + mkt_app pxO [mkt_pos (quo_num n num_2)] + else + mkt_app pxI [mkt_pos (quo_num n num_2)] + +let mkt_n n = + if n=num_0 + then Lazy.force nN0 + else mkt_app nNpos [mkt_pos n] + +let mkt_z z = + if z =/ num_0 then Lazy.force z0 + else if z >/ num_0 then + mkt_app zpos [mkt_pos z] + else + mkt_app zneg [mkt_pos ((Int 0) -/ z)] + +let rec mkt_term t = match t with +| Zero -> mkt_term (Const num_0) +| Const r -> let (n,d) = numdom r in + mkt_app ttconst [Lazy.force tz; mkt_z n] +| Var v -> mkt_app ttvar [Lazy.force tz; mkt_pos (num_of_string v)] +| Opp t1 -> mkt_app ttopp [Lazy.force tz; mkt_term t1] +| Add (t1,t2) -> mkt_app ttadd [Lazy.force tz; mkt_term t1; mkt_term t2] +| Sub (t1,t2) -> mkt_app ttsub [Lazy.force tz; mkt_term t1; mkt_term t2] +| Mul (t1,t2) -> mkt_app ttmul [Lazy.force tz; mkt_term t1; mkt_term t2] +| Pow (t1,n) -> if (n = 0) then + mkt_app ttconst [Lazy.force tz; mkt_z num_1] +else + mkt_app ttpow [Lazy.force tz; mkt_term t1; mkt_n (num_of_int n)] + +let rec parse_pos p = + match kind_of_term p with +| App (a,[|p2|]) -> + if a = Lazy.force pxO then num_2 */ (parse_pos p2) + else num_1 +/ (num_2 */ (parse_pos p2)) +| _ -> num_1 + +let parse_z z = + match kind_of_term z with +| App (a,[|p2|]) -> + if a = Lazy.force zpos then parse_pos p2 else (num_0 -/ (parse_pos p2)) +| _ -> num_0 + +let parse_n z = + match kind_of_term z with +| App (a,[|p2|]) -> + parse_pos p2 +| _ -> num_0 + +let rec parse_term p = + match kind_of_term p with +| App (a,[|_;p2|]) -> + if a = Lazy.force ttvar then Var (string_of_num (parse_pos p2)) + else if a = Lazy.force ttconst then Const (parse_z p2) + else if a = Lazy.force ttopp then Opp (parse_term p2) + else Zero +| App (a,[|_;p2;p3|]) -> + if a = Lazy.force ttadd then Add (parse_term p2, parse_term p3) + else if a = Lazy.force ttsub then Sub (parse_term p2, parse_term p3) + else if a = Lazy.force ttmul then Mul (parse_term p2, parse_term p3) + else if a = Lazy.force ttpow then + Pow (parse_term p2, int_of_num (parse_n p3)) + else Zero +| _ -> Zero + +let rec parse_request lp = + match kind_of_term lp with + | App (_,[|_|]) -> [] + | App (_,[|_;p;lp1|]) -> + (parse_term p)::(parse_request lp1) + |_-> assert false + +let nvars = ref 0 + +let set_nvars_term t = + let rec aux t = + match t with + | Zero -> () + | Const r -> () + | Var v -> let n = int_of_string v in + nvars:= max (!nvars) n + | Opp t1 -> aux t1 + | Add (t1,t2) -> aux t1; aux t2 + | Sub (t1,t2) -> aux t1; aux t2 + | Mul (t1,t2) -> aux t1; aux t2 + | Pow (t1,n) -> aux t1 + in aux t + +let string_of_term p = + let rec aux p = + match p with + | Zero -> "0" + | Const r -> string_of_num r + | Var v -> "x"^v + | Opp t1 -> "(-"^(aux t1)^")" + | Add (t1,t2) -> "("^(aux t1)^"+"^(aux t2)^")" + | Sub (t1,t2) -> "("^(aux t1)^"-"^(aux t2)^")" + | Mul (t1,t2) -> "("^(aux t1)^"*"^(aux t2)^")" + | Pow (t1,n) -> (aux t1)^"^"^(string_of_int n) + in aux p + + +(*********************************************************************** + Coefficients: recursive polynomials + *) + +module Coef = BigInt +(*module Coef = Ent*) +module Poly = Polynom.Make(Coef) +module PIdeal = Ideal.Make(Poly) +open PIdeal + +(* term to sparse polynomial + varaibles <=np are in the coefficients +*) + +let term_pol_sparse np t= + let d = !nvars in + let rec aux t = + match t with + | Zero -> zeroP + | Const r -> + if r = num_0 + then zeroP + else polconst d (Poly.Pint (Coef.of_num r)) + | Var v -> + let v = int_of_string v in + if v <= np + then polconst d (Poly.x v) + else gen d v + | Opp t1 -> oppP (aux t1) + | Add (t1,t2) -> plusP (aux t1) (aux t2) + | Sub (t1,t2) -> plusP (aux t1) (oppP (aux t2)) + | Mul (t1,t2) -> multP (aux t1) (aux t2) + | Pow (t1,n) -> puisP (aux t1) n + in (*info ("conversion de: "^(string_of_term t)^"\n");*) + let res= aux t in + (*info ("donne: "^(stringP res)^"\n");*) + res + +(* sparse polynomial to term *) + +let polrec_to_term p = + let rec aux p = + match p with + |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), + i)))) + coefs; + !res + in aux p + +(* approximation of the Horner form used in the tactic ring *) + +let pol_sparse_to_term n2 p = + info "pol_sparse_to_term ->\n"; + 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) + 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 info "-> pol_sparse_to_term\n"; + aux p + + +let rec remove_list_tail l i = + let rec aux l i = + if l=[] + then [] + else if i<0 + then l + else if i=0 + then List.tl l + else + match l with + |(a::l1) -> + a::(aux l1 (i-1)) + |_ -> assert false + in + List.rev (aux (List.rev l) i) + +(* + lq = [cn+m+1 n+m ...cn+m+1 1] + lci=[[cn+1 n,...,cn1 1] + ... + [cn+m n+m-1,...,cn+m 1]] + + removes intermediate polynomials not useful to compute the last one. + *) + +let remove_zeros zero lci = + let n = List.length (List.hd lci) in + let m=List.length lci in + let u = Array.create m false in + let rec utiles k = + if k>=m + then () + else ( + u.(k)<-true; + let lc = List.nth lci k in + for i=0 to List.length lc - 1 do + if not (zero (List.nth lc i)) + then utiles (i+k+1); + done) + in utiles 0; + let lr = ref [] in + for i=0 to m-1 do + if u.(i) + then lr:=(List.nth lci i)::(!lr) + done; + let lr=List.rev !lr in + let lr = List.map + (fun lc -> + let lcr=ref lc in + for i=0 to m-1 do + if not u.(i) + then lcr:=remove_list_tail !lcr (m-i+(n-m)) + done; + !lcr) + lr in + info ("unuseful spolynomials: " + ^string_of_int (m-List.length lr)^"\n"); + info ("useful spolynomials: " + ^string_of_int (List.length lr)^"\n"); + lr + +let theoremedeszeros lpol p = + let t1 = Unix.gettimeofday() in + let m = !nvars in + let (lp0,p,cert) = in_ideal m lpol p in + let lpc = List.rev !poldepcontent in + info ("time: "^Format.sprintf "@[%10.3f@]s\n" (Unix.gettimeofday ()-.t1)); + (cert,lp0,p,lpc) + +open Ideal + +let theoremedeszeros_termes lp = + nvars:=0;(* mise a jour par term_pol_sparse *) + List.iter set_nvars_term lp; + match lp with + | Const (Int sugarparam)::Const (Int nparam)::lp -> + ((match sugarparam with + |0 -> info "calcul sans sugar\n"; + lexico:=false; + sugar_flag := false; + divide_rem_with_critical_pair := false + |1 -> info "calcul avec sugar\n"; + lexico:=false; + sugar_flag := true; + divide_rem_with_critical_pair := false + |2 -> info "ordre lexico calcul sans sugar\n"; + lexico:=true; + sugar_flag := false; + divide_rem_with_critical_pair := false + |3 -> info "ordre lexico calcul avec sugar\n"; + lexico:=true; + sugar_flag := true; + divide_rem_with_critical_pair := false + |4 -> info "calcul sans sugar, division par les paires\n"; + lexico:=false; + sugar_flag := false; + divide_rem_with_critical_pair := true + |5 -> info "calcul avec sugar, division par les paires\n"; + lexico:=false; + sugar_flag := true; + divide_rem_with_critical_pair := true + |6 -> info "ordre lexico calcul sans sugar, division par les paires\n"; + lexico:=true; + sugar_flag := false; + divide_rem_with_critical_pair := true + |7 -> info "ordre lexico calcul avec sugar, division par les paires\n"; + lexico:=true; + sugar_flag := true; + divide_rem_with_critical_pair := true + | _ -> error "nsatz: bad parameter" + ); + let m= !nvars in + let lvar=ref [] in + for i=m downto 1 do lvar:=["x"^(string_of_int i)^""]@(!lvar); done; + lvar:=["a";"b";"c";"d";"e";"f";"g";"h";"i";"j";"k";"l";"m";"n";"o";"p";"q";"r";"s";"t";"u";"v";"w";"x";"y";"z"] @ (!lvar); (* pour macaulay *) + name_var:=!lvar; + let lp = List.map (term_pol_sparse nparam) lp in + match lp with + | [] -> assert false + | p::lp1 -> + let lpol = List.rev lp1 in + 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 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,r,lci,lq) + ) + |_ -> assert false + + +(* version avec hash-consing du certificat: +let nsatz 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 nsatz lpol = + let lp= parse_request lpol 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 + (fun lt r -> + let ltterm = + List.fold_right + (fun t r -> + mkt_app lcons [mkt_app tpexpr [Lazy.force tz];t;r]) + lt + (mkt_app lnil [mkt_app tpexpr [Lazy.force tz]]) in + mkt_app lcons [tlp ();ltterm;r]) + res + (mkt_app lnil [tlp ()]) in + info "terme calcule\n"; + res + +let return_term t = + let a = + mkApp(gen_constant "CC" ["Init";"Logic"] "refl_equal",[|tllp ();t|]) in + generalize [a] + +let nsatz_compute t = + let lpol = + try nsatz t + with Ideal.NotInIdeal -> + error "nsatz cannot solve this problem" in + return_term lpol + +TACTIC EXTEND nsatz_compute +| [ "nsatz_compute" constr(lt) ] -> [ nsatz_compute lt ] +END + + diff --git a/plugins/nsatz/nsatz_plugin.mllib b/plugins/nsatz/nsatz_plugin.mllib new file mode 100644 index 000000000..a25e649d0 --- /dev/null +++ b/plugins/nsatz/nsatz_plugin.mllib @@ -0,0 +1,5 @@ +Utile +Polynom +Ideal +Nsatz +Nsatz_plugin_mod diff --git a/plugins/nsatz/polynom.ml b/plugins/nsatz/polynom.ml new file mode 100644 index 000000000..14e279b56 --- /dev/null +++ b/plugins/nsatz/polynom.ml @@ -0,0 +1,679 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* 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 hash : t -> int + module Hashpol : Hashtbl.S with type key=t +end + +(*********************************************************************** + 2. Type of polynomials, operations. +*) +module Make (C:Coef) = struct + +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 + +type variable = int + +type t = + Pint of coef (* constant polynomial *) + | Prec of variable * (t array) (* coefficients, increasing degree *) + +(* by default, operations work with normalized polynomials: +- variables are positive integers +- coefficients of a polynomial in x only use variables < x +- no zero coefficient at beginning +- no Prec(x,a) where a is constant in x +*) + +(* constant polynomials *) +let of_num x = Pint (C.of_num x) +let cf0 = of_num (Num.Int 0) +let cf1 = of_num (Num.Int 1) + +(* nth variable *) +let x n = Prec (n,[|cf0;cf1|]) + +(* create 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) + +let is_constantP = function + Pint _ -> true + | Prec _ -> false + +let int_of_Pint = function + Pint x -> x + | _ -> failwith "non" + +let is_zero p = + match p with Pint n -> if C.equal n coef0 then true else false |_-> false + +let max_var_pol p = + match p with + Pint _ -> 0 + |Prec(x,_) -> x + +(* p not normalized *) +let rec max_var_pol2 p = + match p with + Pint _ -> 0 + |Prec(v,c)-> Array.fold_right (fun q m -> max (max_var_pol2 q) m) c v + +let rec max_var l = Array.fold_right (fun p m -> max (max_var_pol2 p) m) l 0 + +(* equality between polynomials *) + +let rec equal p q = + match (p,q) with + (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 (equal a q1.(i)) + then failwith "raté") + p1; + true) + with _ -> false) + | (_,_) -> false + +(* normalize polynomial: remove head zeros, coefficients are normalized + if constant, returns the coefficient +*) + +let rec norm p = match p with + Pint _ -> p + |Prec (x,a)-> + let d = (Array.length a -1) in + let n = ref d in + while !n>0 && (equal a.(!n) (Pint coef0)) do + n:=!n-1; + done; + if !n<0 then Pint coef0 + else if !n=0 then a.(0) + else if !n=d then p + else (let b=Array.create (!n+1) (Pint coef0) in + for i=0 to !n do b.(i)<-a.(i);done; + Prec(x,b)) + + +(* degree in v, v >= max var of p *) +let rec deg v p = + match p with + Prec(x,p1) when x=v -> Array.length p1 -1 + |_ -> 0 + + +(* total degree *) +let rec deg_total p = + match p with + Prec (x,p1) -> let d = ref 0 in + Array.iteri (fun i q -> d:= (max !d (i+(deg_total q)))) p1; + !d + |_ -> 0 + +let rec copyP p = + match p with + Pint i -> Pint i + |Prec(x,q) -> Prec(x,Array.map copyP q) + +(* coefficient of degree i in v, v >= max var of p *) +let coef v i p = + match p with + Prec (x,p1) when x=v -> if i<(Array.length p1) then p1.(i) else Pint coef0 + |_ -> if i=0 then p else Pint coef0 + +(* addition *) + +let rec plusP p q = + let res = + (match (p,q) with + (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) + |(Prec (x,p1),Pint b) -> let p2=Array.map copyP p1 in + p2.(0)<- plusP p1.(0) q; + Prec (x,p2) + |(Prec (x,p1),Prec (y,q1)) -> + if xy then (let p2=Array.map copyP p1 in + p2.(0)<- plusP p1.(0) q; + Prec (x,p2)) + else + (let n=max (deg x p) (deg x q) in + let r=Array.create (n+1) (Pint coef0) in + for i=0 to n do + r.(i)<- plusP (coef x i p) (coef x i q); + done; + Prec(x,r))) + in norm res + + +(* content, positive integer *) +let rec content p = + match p with + Pint a -> C.abs a + | Prec (x ,p1) -> + Array.fold_left C.pgcd coef0 (Array.map content p1) + +let rec div_int p a= + match p with + Pint b -> Pint (C.div b a) + | Prec(x,p1) -> Prec(x,Array.map (fun x -> div_int x a) p1) + +let vire_contenu p = + let c = content p in + if C.equal c coef0 then p else div_int p c + +(* sorted list of variables of a polynomial *) + +let rec vars=function + Pint _->[] + | Prec (x,l)->(List.flatten ([x]::(List.map vars (Array.to_list l)))) + + +(* multiply p by v^n, v >= max_var p *) +let rec multx n v p = + match p with + Prec (x,p1) when x=v -> let p2= Array.create ((Array.length p1)+n) (Pint coef0) in + for i=0 to (Array.length p1)-1 do + p2.(i+n)<-p1.(i); + done; + Prec (x,p2) + |_ -> if p = (Pint coef0) then (Pint coef0) + else (let p2=Array.create (n+1) (Pint coef0) in + p2.(n)<-p; + Prec (v,p2)) + + +(* product *) +let rec multP p q = + match (p,q) with + (Pint a,Pint b) -> Pint (C.mult a b) + |(Pint a, Prec (y,q1)) -> + 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 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)) -> + if x multP p z) q1 in + Prec (y,q2)) + else if x>y + then (let p2 = Array.map (fun z-> multP z q) p1 in + Prec (x,p2)) + else Array.fold_left plusP (Pint coef0) + (Array.mapi (fun i z-> (multx i x (multP z q))) p1) + + + +(* derive p with variable v, v >= max_var p *) +let rec deriv v p = + match p with + Pint a -> Pint coef0 + | Prec(x,p1) when x=v -> + let d = Array.length p1 -1 in + if d=1 then p1.(1) + else + (let p2 = Array.create d (Pint coef0) in + for i=0 to d-1 do + p2.(i)<- multP (Pint (coef_of_int (i+1))) p1.(i+1); + done; + Prec (x,p2)) + | Prec(x,p1)-> Pint coef0 + + +(* opposite *) +let rec oppP p = + match p with + Pint a -> Pint (C.opp a) + |Prec(x,p1) -> Prec(x,Array.map oppP p1) + +let moinsP p q=plusP p (oppP q) + +let rec puisP p n = match n with + 0 -> cf1 + |_ -> (multP p (puisP p (n-1))) + + +(* infix notations *) +(*let (++) a b = plusP a b +*) +let (@@) a b = multP a b + +let (--) a b = moinsP a b + +let (^^) a b = puisP a b + + +(* leading coefficient in v, v>= max_var p *) + +let coefDom v p= coef v (deg v p) p + +let coefConst v p = coef v 0 p + +(* tail of a polynomial *) +let remP v p = + moinsP p (multP (coefDom v p) (puisP (x v) (deg v p))) + + +(* first interger coefficient of p *) +let rec coef_int_tete p = + let v = max_var_pol p in + if v>0 + then coef_int_tete (coefDom v p) + else (match p with | Pint a -> a |_ -> assert false) + + +(* divide by the content and make the head int coef positive *) +let normc p = + let p = vire_contenu p in + let a = coef_int_tete p in + if C.le coef0 a then p else oppP p + + +(* constant coef of normalized polynomial *) +let rec coef_constant p = + match p with + Pint a->a + |Prec(_,q)->coef_constant q.(0) + + +(*********************************************************************** + 3. Printing polynomials. +*) + +(* if univ = false, we use x,y,z,a,b,c,d... as variables, else x1,x2,... +*) +let univ=ref true + +let string_of_var x= + if !univ then + "u"^(string_of_int x) + else + if x<=3 then String.make 1 (Char.chr(x+(Char.code 'w'))) + else String.make 1 (Char.chr(x-4+(Char.code 'a'))) + +let nsP = ref 0 + +let rec string_of_Pcut p = + if (!nsP)<=0 + then "..." + else + match p with + |Pint a-> nsP:=(!nsP)-1; + 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 "" + and sp=ref "" in + let st0 = string_of_Pcut t.(0) in + if st0<>"0" + then s:=st0; + let fin = ref false in + for i=(Array.length t)-1 downto 1 do + if (!nsP)<0 + then (sp:="..."; + if not (!fin) then s:=(!s)^"+"^(!sp); + fin:=true) + else ( + let si=string_of_Pcut t.(i) in + sp:=""; + if i=1 + then ( + if si<>"0" + then (nsP:=(!nsP)-1; + if si="1" + then sp:=v + else + (if (String.contains si '+') + then sp:="("^si^")*"^v + else sp:=si^"*"^v))) + else ( + if si<>"0" + then (nsP:=(!nsP)-1; + if si="1" + then sp:=v^"^"^(string_of_int i) + else (if (String.contains si '+') + then sp:="("^si^")*"^v^"^"^(string_of_int i) + else sp:=si^"*"^v^"^"^(string_of_int i)))); + if !sp<>"" && not (!fin) + then (nsP:=(!nsP)-1; + if !s="" + then s:=!sp + else s:=(!s)^"+"^(!sp))); + done; + if !s="" then (nsP:=(!nsP)-1; + (s:="0")); + !s + +let to_string p = + nsP:=20; + string_of_Pcut 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)^(to_string p)^"\n") lp; + prt0 ((!s)^"}") + +let print_lpoly lp = print_tpoly (Array.of_list lp) + +(*********************************************************************** + 4. Exact division of polynomials. +*) + +(* return (s,r) s.t. p = s*q+r *) +let rec quo_rem_pol p q x = + if x=0 + then (match (p,q) with + |(Pint a, Pint b) -> + if C.equal (C.modulo a b) coef0 + then (Pint (C.div a b), cf0) + else failwith "div_pol1" + |_ -> assert false) + else + let m = deg x q in + let b = coefDom x q in + let q1 = remP x q in (* q = b*x^m+q1 *) + let r = ref p in + let s = ref cf0 in + let continue =ref true in + while (!continue) && (not (equal !r cf0)) do + let n = deg x !r in + if n false + +(*********************************************************************** + 5. Pseudo-division and gcd with subresultants. +*) + +(* pseudo division : + q = c*x^m+q1 + retruns (r,c,d,s) s.t. c^d*p = s*q + r. +*) + +let pseudo_div p q x = + match q with + Pint _ -> (cf0, q,1, p) + | Prec (v,q1) when x<>v -> (cf0, q,1, p) + | Prec (v,q1) -> + ( + (* pr "pseudo_division: c^d*p = s*q + r";*) + let delta = ref 0 in + let r = ref p in + let c = coefDom x q in + let q1 = remP x q in + let d' = deg x q in + let s = ref cf0 in + while (deg x !r)>=(deg x q) do + let d = deg x !r in + let a = coefDom x !r in + let r1=remP x !r in + let u = a @@ ((monome x (d-d'))) in + r:=(c @@ r1) -- (u @@ q1); + s:=plusP (c @@ (!s)) u; + delta := (!delta) + 1; + done; + (* + pr ("deg d: "^(string_of_int (!delta))^", deg c: "^(string_of_int (deg_total c))); + pr ("deg r:"^(string_of_int (deg_total !r))); + *) + (!r,c,!delta, !s) + ) + +(* gcd with subresultants *) + +let rec pgcdP p q = + let x = max (max_var_pol p) (max_var_pol q) in + pgcd_pol p q x + +and pgcd_pol p q x = + pgcd_pol_rec p q x + +and content_pol p x = + match p with + Prec(v,p1) when v=x -> + Array.fold_left (fun a b -> pgcd_pol_rec a b (x-1)) cf0 p1 + | _ -> p + +and pgcd_coef_pol c p x = + match p with + Prec(v,p1) when x=v -> + Array.fold_left (fun a b -> pgcd_pol_rec a b (x-1)) c p1 + |_ -> pgcd_pol_rec c p (x-1) + +and pgcd_pol_rec p q x = + match (p,q) with + (Pint a,Pint b) -> Pint (C.pgcd (C.abs a) (C.abs b)) + |_ -> + if equal p cf0 + then q + else if equal 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 + ) + +(* Sub-résultants: + + ai*Ai = Qi*Ai+1 + bi*Ai+2 + + deg Ai+2 < deg Ai+1 + + Ai = ci*X^ni + ... + di = ni - ni+1 + + ai = (- ci+1)^(di + 1) + b1 = 1 + bi = ci*si^di si i>1 + + s1 = 1 + si+1 = ((ci+1)^di*si)/si^di + +*) +and gcd_sub_res p q x = + if equal q cf0 + then p + else + let d = deg x p in + let d' = deg x q in + if d (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 = equal + let hash = hash + end) + +end diff --git a/plugins/nsatz/polynom.mli b/plugins/nsatz/polynom.mli new file mode 100644 index 000000000..623d901e8 --- /dev/null +++ b/plugins/nsatz/polynom.mli @@ -0,0 +1,97 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* 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 hash : t -> int + module Hashpol : Hashtbl.S with type key=t +end + +module Make (C:Coef) : S with type coef = C.t diff --git a/plugins/nsatz/utile.ml b/plugins/nsatz/utile.ml new file mode 100644 index 000000000..c16bd4256 --- /dev/null +++ b/plugins/nsatz/utile.ml @@ -0,0 +1,130 @@ +(* Printing *) + +let pr x = + if !Flags.debug then (Format.printf "@[%s@]" x; flush(stdout);)else () + +let prn x = + if !Flags.debug then (Format.printf "@[%s\n@]" x; flush(stdout);) else () + +let prt0 s = () (* print_string s;flush(stdout)*) + +let prt s = + if !Flags.debug then (print_string (s^"\n");flush(stdout)) else () + +let info s = + Flags.if_verbose prerr_string s + +(* Lists *) + +let rec list_mem_eq eq x l = + match l with + [] -> false + |y::l1 -> if (eq x y) then true else (list_mem_eq eq x l1) + +let set_of_list_eq eq l = + let res = ref [] in + List.iter (fun x -> if not (list_mem_eq eq x (!res)) then res:=x::(!res)) l; + List.rev !res + + +(* Memoization + f is compatible with nf: f(nf(x)) = f(x) +*) + +let memos s memoire nf f x = + try (let v = Hashtbl.find memoire (nf x) in pr s;v) + with _ -> (pr "#"; + let v = f x in + Hashtbl.add memoire (nf x) v; + v) + + +(********************************************************************** + Eléments minimaux pour un ordre partiel de division. + E est un ensemble, avec une multiplication + et une division partielle div (la fonction div peut échouer), + constant est un prédicat qui définit un sous-ensemble C de E. +*) +(* + Etant donnée une partie A de E, on calcule une partie B de E disjointe de C + telle que: + - les éléments de A sont des produits d'éléments de B et d'un de C. + - B est minimale pour cette propriété. +*) + +let facteurs_liste div constant lp = + let lp = List.filter (fun x -> not (constant x)) lp in + let rec factor lmin lp = (* lmin: ne se divisent pas entre eux *) + match lp with + [] -> lmin + |p::lp1 -> + (let l1 = ref [] in + let p_dans_lmin = ref false in + List.iter (fun q -> try (let r = div p q in + if not (constant r) + then l1:=r::(!l1) + else p_dans_lmin:=true) + with _ -> ()) + lmin; + if !p_dans_lmin + then factor lmin lp1 + else if (!l1)=[] + (* aucun q de lmin ne divise p *) + then (let l1=ref lp1 in + let lmin1=ref [] in + List.iter (fun q -> try (let r = div q p in + if not (constant r) + then l1:=r::(!l1)) + with _ -> lmin1:=q::(!lmin1)) + lmin; + factor (List.rev (p::(!lmin1))) !l1) + (* au moins un q de lmin divise p non trivialement *) + else factor lmin ((!l1)@lp1)) + in + factor [] lp + + +(* On suppose que tout élément de A est produit d'éléments de B et d'un de C: + A et B sont deux tableaux, rend un tableau de couples + (élément de C, listes d'indices l) + tels que A.(i) = l.(i)_1*Produit(B.(j), j dans l.(i)_2) + zero est un prédicat sur E tel que (zero x) => (constant x): + si (zero x) est vrai on ne decompose pas x + c est un élément quelconque de E. +*) +let factorise_tableau div zero c f l1 = + let res = Array.create (Array.length f) (c,[]) in + Array.iteri (fun i p -> + let r = ref p in + let li = ref [] in + if not (zero p) + then + Array.iteri (fun j q -> + try (while true do + let rr = div !r q in + li:=j::(!li); + r:=rr; + done) + with _ -> ()) + l1; + res.(i)<-(!r,!li)) + f; + (l1,res) + + +(* exemples: + +let l = [1;2;6;24;720] +and div1 = (fun a b -> if a mod b =0 then a/b else failwith "div") +and constant = (fun x -> x<2) +and zero = (fun x -> x=0) + + +let f = facteurs_liste div1 constant l + + +factorise_tableau div1 zero 0 (Array.of_list l) (Array.of_list f) + +*) + + diff --git a/plugins/nsatz/utile.mli b/plugins/nsatz/utile.mli new file mode 100644 index 000000000..83b2ac39b --- /dev/null +++ b/plugins/nsatz/utile.mli @@ -0,0 +1,22 @@ + +(* Printing *) +val pr : string -> unit +val prn : string -> unit +val prt0 : 'a -> unit +val prt : string -> unit +val info : string -> unit + +(* Listes *) +val list_mem_eq : ('a -> 'b -> bool) -> 'a -> 'b list -> bool +val set_of_list_eq : ('a -> 'a -> bool) -> 'a list -> 'a list + +(* Memoization *) +val memos : + string -> ('a, 'b) Hashtbl.t -> ('c -> 'a) -> ('c -> 'b) -> 'c -> 'b + + +val facteurs_liste : ('a -> 'a -> 'a) -> ('a -> bool) -> 'a list -> 'a list +val factorise_tableau : + ('a -> 'b -> 'a) -> + ('a -> bool) -> + 'a -> 'a array -> 'b array -> 'b array * ('a * int list) array diff --git a/plugins/nsatz/vo.itarget b/plugins/nsatz/vo.itarget new file mode 100644 index 000000000..80d69e4af --- /dev/null +++ b/plugins/nsatz/vo.itarget @@ -0,0 +1 @@ +NsatzR.vo -- cgit v1.2.3