aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--CHANGES4
-rw-r--r--CREDITS2
-rw-r--r--Makefile.common14
-rw-r--r--doc/refman/Reference-Manual.tex1
-rw-r--r--plugins/groebner/GroebnerR.v410
-rw-r--r--plugins/groebner/GroebnerZ.v73
-rw-r--r--plugins/groebner/groebner.ml4571
-rw-r--r--plugins/groebner/groebner_plugin.mllib5
-rw-r--r--plugins/groebner/ideal.ml1356
-rw-r--r--plugins/groebner/ideal.mli80
-rw-r--r--plugins/groebner/polynom.ml1051
-rw-r--r--plugins/groebner/polynom.mli120
-rw-r--r--plugins/groebner/utile.ml161
-rw-r--r--plugins/groebner/utile.mli24
-rw-r--r--plugins/groebner/vo.itarget2
-rw-r--r--plugins/pluginsbyte.itarget2
-rw-r--r--plugins/pluginsdyn.itarget2
-rw-r--r--plugins/pluginsopt.itarget2
-rw-r--r--plugins/pluginsvo.itarget3
-rw-r--r--test-suite/success/Nsatz.v216
20 files changed, 230 insertions, 3869 deletions
diff --git a/CHANGES b/CHANGES
index 6e1e218f5..3212a4d42 100644
--- a/CHANGES
+++ b/CHANGES
@@ -38,8 +38,8 @@ Automation tactics
hints (wish #2104).
- An inductive type as argument of the "using" option of "auto/eauto/firstorder"
is interpreted as using the collection of its constructors.
-- New decision tactic "gb" to solve systems of polynomial equations
- by (exact) computation of Groebner bases.
+- New decision tactic "nsatz" to prove polynomial equations
+ by computation of Groebner bases.
Other tactics
diff --git a/CREDITS b/CREDITS
index 68009d390..53bd9e93c 100644
--- a/CREDITS
+++ b/CREDITS
@@ -47,7 +47,7 @@ plugins/funind
and Yves Bertot (INRIA-Marelle, 2005-2006)
plugins/omega
developed by Pierre Crégut (France Telecom R&D, 1996)
-plugins/groebner
+plugins/nsatz
developed by Loïc Pottier (INRIA-Marelle, 2009)
plugins/ring
developed by Samuel Boutin (INRIA-Coq, 1996) and Patrick
diff --git a/Makefile.common b/Makefile.common
index 046a85c62..7ccd5c901 100644
--- a/Makefile.common
+++ b/Makefile.common
@@ -73,7 +73,7 @@ SRCDIRS:=\
omega romega micromega quote ring dp \
setoid_ring xml extraction fourier \
cc funind firstorder field subtac \
- rtauto groebner syntax decl_mode)
+ rtauto nsatz syntax decl_mode)
# Order is relevent here because kernel and checker contain files
# with the same name
@@ -118,7 +118,7 @@ REFMANCOQTEXFILES:=$(addprefix doc/refman/, \
RefMan-oth.v.tex RefMan-ltac.v.tex \
RefMan-decl.v.tex \
Cases.v.tex Coercion.v.tex Extraction.v.tex \
- Program.v.tex Omega.v.tex Micromega.v.tex Polynom.v.tex \
+ Program.v.tex Omega.v.tex Micromega.v.tex Polynom.v.tex Nsatz.v.tex \
Setoid.v.tex Helm.tex Classes.v.tex )
REFMANTEXFILES:=$(addprefix doc/refman/, \
@@ -166,7 +166,7 @@ MICROMEGACMA:=plugins/micromega/micromega_plugin.cma
QUOTECMA:=plugins/quote/quote_plugin.cma
RINGCMA:=plugins/ring/ring_plugin.cma
NEWRINGCMA:=plugins/setoid_ring/newring_plugin.cma
-GBCMA:=plugins/groebner/groebner_plugin.cma
+NSATZCMA:=plugins/nsatz/nsatz_plugin.cma
DPCMA:=plugins/dp/dp_plugin.cma
FIELDCMA:=plugins/field/field_plugin.cma
XMLCMA:=plugins/xml/xml_plugin.cma
@@ -191,7 +191,7 @@ PLUGINSCMA:=$(OMEGACMA) $(ROMEGACMA) $(MICROMEGACMA) \
$(FOURIERCMA) $(EXTRACTIONCMA) $(XMLCMA) \
$(DECLMODECMA) \
$(CCCMA) $(FOCMA) $(SUBTACCMA) $(RTAUTOCMA) \
- $(FUNINDCMA) $(GBCMA) $(NATSYNTAXCMA) $(OTHERSYNTAXCMA)
+ $(FUNINDCMA) $(NSATZCMA) $(NATSYNTAXCMA) $(OTHERSYNTAXCMA)
ifneq ($(HASNATDYNLINK),false)
STATICPLUGINS:=
@@ -296,20 +296,18 @@ QUOTEVO:=$(call cat_vo_itarget, plugins/quote)
RINGVO:=$(call cat_vo_itarget, plugins/ring)
FIELDVO:=$(call cat_vo_itarget, plugins/field)
NEWRINGVO:=$(call cat_vo_itarget, plugins/setoid_ring)
-GBVO:=$(call cat_vo_itarget, plugins/groebner)
+NSATZVO:=$(call cat_vo_itarget, plugins/nsatz)
FOURIERVO:=$(call cat_vo_itarget, plugins/fourier)
FUNINDVO:=$(call cat_vo_itarget, plugins/funind)
DPVO:=$(call cat_vo_itarget, plugins/dp)
RTAUTOVO:=$(call cat_vo_itarget, plugins/rtauto)
-EXTRACTIONVO:=$(call cat_vo_itarget, plugins/extraction)
XMLVO:=
CCVO:=
-
PLUGINSVO:= $(OMEGAVO) $(ROMEGAVO) $(MICROMEGAVO) $(RINGVO) $(FIELDVO) \
$(XMLVO) $(FOURIERVO) $(CCVO) $(FUNINDVO) \
$(RTAUTOVO) $(NEWRINGVO) $(DPVO) $(QUOTEVO) \
- $(GBVO) $(EXTRACTIONVO)
+ $(NSATZVO)
ALLVO:= $(THEORIESVO) $(PLUGINSVO)
VFILES:= $(ALLVO:.vo=.v)
diff --git a/doc/refman/Reference-Manual.tex b/doc/refman/Reference-Manual.tex
index 0c90e61ba..d9641c2c8 100644
--- a/doc/refman/Reference-Manual.tex
+++ b/doc/refman/Reference-Manual.tex
@@ -109,6 +109,7 @@ Options A and B of the licence are {\em not} elected.}
\include{Extraction.v}%
\include{Program.v}%
\include{Polynom.v}% = Ring
+\include{Nsatz.v}%
\include{Setoid.v}% Tactique pour les setoides
%BEGIN LATEX
\RefManCutCommand{ENDADDENDUM=\thepage}
diff --git a/plugins/groebner/GroebnerR.v b/plugins/groebner/GroebnerR.v
deleted file mode 100644
index fc01c5886..000000000
--- a/plugins/groebner/GroebnerR.v
+++ /dev/null
@@ -1,410 +0,0 @@
-(************************************************************************)
-(* 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 *)
-(************************************************************************)
-
-(*
- Tactic groebnerR: proofs of polynomials equalities with variables in R.
- Use Hilbert Nullstellensatz and Buchberger algorithm (adapted version of
- L.Thery Coq proven implementation).
- Thanks to B.Gregoire and L.Thery for help on ring tactic.
- Examples at the end of the file.
-
- 3 versions:
-
- - groebnerR.
-
- - groebnerRp (a::b::c::nil) : give the list of variables are considered as
- parameters. Computation will be performed with rational fractions in these
- variables.
-
- - groebnerRpv (a::b::c::nil) (x::y::z::nil): give also the order of the
- variables used in Buchberger algorithm. Here x>y>z.
-
- Loic Pottier
- 19-12-08
-*)
-
-Require Import List.
-Require Import Setoid.
-Require Import BinPos.
-Require Import BinList.
-Require Import Znumtheory.
-Require Import RealField Rdefinitions Rfunctions RIneq DiscrR.
-Require Import Ring_polynom Ring_tac InitialRing.
-
-Declare ML Module "groebner_plugin".
-
-Local Open Scope R_scope.
-
-Lemma psos_r1b: forall x y, x - y = 0 -> 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 groebnerR_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 groebnerR_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 groebnerR_not2: 1<>0.
-auto with *.
-Qed.
-
-Lemma groebnerR_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 groebnerR_not_hyp:=
- match goal with
- | H: ?x<>?y |- _ =>
- match y with
- |0 =>
- let H1:=fresh "Hgroebner" in
- let y:=fresh "x" in
- destruct (@groebnerR_not1_0 _ H) as (y,H1); clear H
- |_ => generalize (@groebnerR_diff _ _ H); clear H; intro
- end
- end.
-
-Ltac groebnerR_not_goal :=
- match goal with
- | |- ?x<>?y :> R => red; intro; apply groebnerR_not2
- | |- False => apply groebnerR_not2
- end.
-
-Ltac groebnerR_begin :=
- intros;
- repeat groebnerR_not_hyp;
- try groebnerR_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 groebnerR_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 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 := 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);
- 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"
- | 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 groebnerR_l3 with (interpret3 c fv) (Nnat.nat_of_N r);simpl;
- [ discrR || idtac "could not prove discrimination result"
- | exact Hg2]
- ]
- ])).
-
-Ltac groebnerRpv lparam lvar:=
- groebnerR_begin;
- intros;
- let G := Get_goal in
- ring_lookup
- (PackRing ltac:(groebnerR_gen lparam lvar ring_subst_niter))
- [] G.
-
-Ltac groebnerR := groebnerRpv (@nil R) (@nil R).
-
-Ltac groebnerRp lparam := groebnerRpv lparam (@nil R).
-
-(**** Examples : see test-suite/success/Groebner.v ****)
diff --git a/plugins/groebner/GroebnerZ.v b/plugins/groebner/GroebnerZ.v
deleted file mode 100644
index 7c40bbb70..000000000
--- a/plugins/groebner/GroebnerZ.v
+++ /dev/null
@@ -1,73 +0,0 @@
-(************************************************************************)
-(* 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 *)
-(************************************************************************)
-
-Require Import Reals ZArith.
-Require Export GroebnerR.
-
-Open Scope Z_scope.
-
-Lemma groebnerZhypR: forall x y:Z, x=y -> IZR x = IZR y.
-Proof IZR_eq. (* or f_equal ... *)
-
-Lemma groebnerZconclR: forall x y:Z, IZR x = IZR y -> x = y.
-Proof eq_IZR.
-
-Lemma groebnerZhypnotR: forall x y:Z, x<>y -> IZR x <> IZR y.
-Proof IZR_neq.
-
-Lemma groebnerZconclnotR: forall x y:Z, IZR x <> IZR y -> x <> y.
-Proof.
-intros x y H. contradict H. f_equal. assumption.
-Qed.
-
-Ltac groebnerZversR1 :=
- repeat
- (match goal with
- | H:(@eq Z ?x ?y) |- _ =>
- generalize (@groebnerZhypR _ _ H); clear H; intro H
- | |- (@eq Z ?x ?y) => apply groebnerZconclR
- | H:not (@eq Z ?x ?y) |- _ =>
- generalize (@groebnerZhypnotR _ _ H); clear H; intro H
- | |- not (@eq Z ?x ?y) => apply groebnerZconclnotR
- end).
-
-Lemma groebnerZR1: forall x y:Z, IZR(x+y) = (IZR x + IZR y)%R.
-Proof plus_IZR.
-
-Lemma groebnerZR2: forall x y:Z, IZR(x*y) = (IZR x * IZR y)%R.
-Proof mult_IZR.
-
-Lemma groebnerZR3: forall x y:Z, IZR(x-y) = (IZR x - IZR y)%R.
-Proof.
-intros; symmetry. apply Z_R_minus.
-Qed.
-
-Lemma groebnerZR4: 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 groebnerZversR2:=
- repeat
- (rewrite groebnerZR1 in * ||
- rewrite groebnerZR2 in * ||
- rewrite groebnerZR3 in * ||
- rewrite groebnerZR4 in *).
-
-Ltac groebnerZ_begin :=
- intros;
- groebnerZversR1;
- groebnerZversR2;
- simpl in *.
- (*cbv beta iota zeta delta [nat_of_P Pmult_nat plus mult] in *.*)
-
-Ltac groebnerZ :=
- groebnerZ_begin; (*idtac "groebnerZ_begin;";*)
- groebnerR.
diff --git a/plugins/groebner/groebner.ml4 b/plugins/groebner/groebner.ml4
deleted file mode 100644
index cc1b08a63..000000000
--- a/plugins/groebner/groebner.ml4
+++ /dev/null
@@ -1,571 +0,0 @@
-(************************************************************************)
-(* v * The Coq Proof Assistant / The Coq Development Team *)
-(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *)
-(* \VV/ **************************************************************)
-(* // * This file is distributed under the terms of the *)
-(* * GNU Lesser General Public License Version 2.1 *)
-(************************************************************************)
-
-(*i camlp4deps: "parsing/grammar.cma" i*)
-
-open Pp
-open Util
-open Names
-open Term
-open Closure
-open Environ
-open Libnames
-open Tactics
-open Rawterm
-open Tacticals
-open Tacexpr
-open Pcoq
-open Tactic
-open Constr
-open Proof_type
-open Coqlib
-open Tacmach
-open Mod_subst
-open Tacinterp
-open Libobject
-open Printer
-open Declare
-open Decl_kinds
-open Entries
-
-open Num
-open Unix
-open Utile
-
-(***********************************************************************
- 1. Opérations sur les coefficients.
-*)
-
-let num_0 = Int 0
-and num_1 = Int 1
-and num_2 = Int 2
-and num_10 = Int 10
-
-let numdom r =
- let r' = Ratio.normalize_ratio (ratio_of_num r) in
- 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?? *)
-(* ------------------------------------------------------------------------- *)
-
-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: 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 *)
-
-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 d (aux t1) (aux t2)
- | Sub (t1,t2) -> plusP d (aux t1) (oppP (aux t2))
- | Mul (t1,t2) -> multP d (aux t1) (aux t2)
- | Pow (t1,n) -> puisP d (aux t1) n
- in (*info ("conversion de: "^(string_of_term t)^"\n");*)
- let res= aux t in
- (*info ("donne: "^(stringP res)^"\n");*)
- res
-
-(* polynome recusrsif vers terme *)
-
-
-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
-
-(* 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 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]]
-
- enleve les polynomes intermediaires inutiles pour calculer le dernier
- *)
-
-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 ("spolynomes inutiles: "
- ^string_of_int (m-List.length lr)^"\n");
- info ("spolynomes utiles: "
- ^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 ("temps: "^Format.sprintf "@[%10.3f@]s\n" (Unix.gettimeofday ()-.t1));
- (cert,lp0,p,lpc)
-
-
-let theoremedeszeros_termes lp =
- nvars:=0;(* mise a jour par term_pol_sparse *)
- List.iter set_nvars_term lp;
- match lp with
- | Const (Int nparam)::lp ->
- (let m= !nvars in
- let lvar=ref [] in
- for i=m downto 1 do lvar:=["x"^string_of_int i^""]@(!lvar); done;
- 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 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,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 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 ]
-END
-
-
diff --git a/plugins/groebner/groebner_plugin.mllib b/plugins/groebner/groebner_plugin.mllib
deleted file mode 100644
index e227b5e09..000000000
--- a/plugins/groebner/groebner_plugin.mllib
+++ /dev/null
@@ -1,5 +0,0 @@
-Utile
-Polynom
-Ideal
-Groebner
-Groebner_plugin_mod
diff --git a/plugins/groebner/ideal.ml b/plugins/groebner/ideal.ml
deleted file mode 100644
index b41d6d8e3..000000000
--- a/plugins/groebner/ideal.ml
+++ /dev/null
@@ -1,1356 +0,0 @@
-(************************************************************************)
-(* 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 *)
-(************************************************************************)
-
-(*
- Nullstellensatz par calcul de base de Grobner
-
- On utilise une representation creuse des polynomes:
- un monome est un tableau d'exposants (un par variable),
- avec son degre en tete.
- un polynome est une liste de (coefficient,monome).
-
- L'algorithme de Buchberger a proprement parler est tire du code caml
- extrait du code Coq ecrit par L.Thery.
-
- *)
-
-open Utile
-
-(* NB: Loic is using a coding-style "let x::q = ... in ..." that
- produces lots of warnings about non-exhaustive pattern matchs.
- Even worse, it is not clear whether a [Match_failure] could
- happen (and be catched by a "with _ ->") or not. Loic told me
- it shouldn't happen.
-
- In a first time, I used a camlp4 extension (search history
- for lib/refutpat.ml4) for introducing an ad-hoc syntax
- "let* x::q = ...". This is now removed (too complex during
- porting to camlp4).
-
- Now, we simply use the following function that turns x::q
- into (x,q) and hence an irrefutable pattern. Yes, this adds
- a (small) cost since the intermediate pair is allocated
- (in opt, the cost might even be 0 due to inlining).
- If somebody want someday to remove this extra cost,
- (x::q) could also be turned brutally into (x,q) by Obj.magic
- (beware: be sure no floats are around in this case).
-*)
-
-let of_cons = function
- | [] -> assert false
- | x::q -> x,q
-
-exception NotInIdeal
-
-module type S = sig
-
-(* Monomials *)
-type mon = int array
-
-val mult_mon : int -> mon -> mon -> mon
-val deg : mon -> int
-val compare_mon : int -> mon -> mon -> int
-val div_mon : int -> mon -> mon -> mon
-val div_mon_test : int -> mon -> mon -> bool
-val ppcm_mon : int -> mon -> mon -> mon
-
-(* Polynomials *)
-
-type deg = int
-type coef
-type poly
-val repr : poly -> (coef * mon) list
-val polconst : deg -> coef -> poly
-val zeroP : poly
-val gen : deg -> int -> poly
-
-val equal : poly -> poly -> bool
-val name_var : string list ref
-val getvar : string list -> int -> string
-val lstringP : poly list -> string
-val printP : poly -> unit
-val lprintP : poly list -> unit
-
-val div_pol_coef : poly -> coef -> poly
-val plusP : deg -> poly -> poly -> poly
-val mult_t_pol : deg -> coef -> mon -> poly -> poly
-val selectdiv : deg -> mon -> poly list -> poly
-val oppP : poly -> poly
-val emultP : coef -> poly -> poly
-val multP : deg -> poly -> poly -> poly
-val puisP : deg -> poly -> int -> poly
-val contentP : poly -> coef
-val contentPlist : poly list -> coef
-val pgcdpos : coef -> coef -> coef
-val div_pol : deg -> poly -> poly -> coef -> coef -> mon -> poly
-val reduce2 : deg -> poly -> poly list -> coef * poly
-
-val poldepcontent : coef list ref
-val coefpoldep_find : poly -> poly -> poly
-val coefpoldep_set : poly -> poly -> poly -> unit
-val initcoefpoldep : deg -> poly list -> unit
-val reduce2_trace : deg -> poly -> poly list -> poly list -> poly list * poly
-val spol : deg -> poly -> poly -> poly
-val etrangers : deg -> poly -> poly -> bool
-val div_ppcm : deg -> poly -> poly -> poly -> bool
-
-val genPcPf : deg -> poly -> poly list -> poly list -> poly list
-val genOCPf : deg -> poly list -> poly list
-
-val is_homogeneous : poly -> bool
-
-type certificate =
- { coef : coef; power : int;
- gb_comb : poly list list; last_comb : poly list }
-val test_dans_ideal : deg -> poly -> poly list -> poly list ->
- poly list * poly * certificate
-val in_ideal : deg -> poly list -> poly -> poly list * poly * certificate
-
-end
-
-(***********************************************************************
- Global options
-*)
-let lexico = ref false
-let use_hmon = ref false
-
-(***********************************************************************
- Functor
-*)
-
-module Make (P:Polynom.S) = struct
-
- type coef = P.t
- let coef0 = P.of_num (Num.Int 0)
- let coef1 = P.of_num (Num.Int 1)
- let coefm1 = P.of_num (Num.Int (-1))
- let string_of_coef c = "["^(P.to_string c)^"]"
-
-(***********************************************************************
- Monomes
- tableau d'entiers
- le premier coefficient est le degre
- *)
-
-type mon = int array
-type deg = int
-type poly = (coef * mon) list
-
-(* d représente le nb de variables du monome *)
-
-(* Multiplication de monomes ayant le meme nb de variables = d *)
-let mult_mon d m m' =
- let m'' = Array.create (d+1) 0 in
- for i=0 to d do
- m''.(i)<- (m.(i)+m'.(i));
- done;
- m''
-
-(* Degré d'un monome *)
-let deg m = m.(0)
-
-
-let compare_mon d m m' =
- if !lexico
- then (
- (* Comparaison de monomes avec ordre du degre lexicographique
- = on commence par regarder la 1ere variable*)
- let res=ref 0 in
- let i=ref 1 in (* 1 si lexico pur 0 si degre*)
- while (!res=0) && (!i<=d) do
- res:= (compare m.(!i) m'.(!i));
- i:=!i+1;
- done;
- !res)
- else (
- (* degre lexicographique inverse *)
- match compare m.(0) m'.(0) with
- | 0 -> (* meme degre total *)
- let res=ref 0 in
- let i=ref d in
- while (!res=0) && (!i>=1) do
- res:= - (compare m.(!i) m'.(!i));
- i:=!i-1;
- done;
- !res
- | x -> x)
-
-(* Division de monome ayant le meme nb de variables *)
-let div_mon d m m' =
- let m'' = Array.create (d+1) 0 in
- for i=0 to d do
- m''.(i)<- (m.(i)-m'.(i));
- done;
- m''
-
-let div_pol_coef p c =
- List.map (fun (a,m) -> (P.divP a c,m)) p
-
-(* m' divise m *)
-let div_mon_test d m m' =
- let res=ref true in
- let i=ref 1 in
- while (!res) && (!i<=d) do
- res:= (m.(!i) >= m'.(!i));
- i:=succ !i;
- done;
- !res
-
-
-(* Mise en place du degré total du monome *)
-let set_deg d m =
- m.(0)<-0;
- for i=1 to d do
- m.(0)<- m.(i)+m.(0);
- done;
- m
-
-
-(* ppcm de monomes *)
-let ppcm_mon d m m' =
- let m'' = Array.create (d+1) 0 in
- for i=1 to d do
- m''.(i)<- (max m.(i) m'.(i));
- done;
- set_deg d m''
-
-
-
-(**********************************************************************
-
- Polynomes
- liste de couples (coefficient,monome) ordonnee en decroissant
-
- *)
-
-
-let repr p = p
-
-let equal =
- Util.list_for_all2eq
- (fun (c1,m1) (c2,m2) -> P.equal c1 c2 && m1=m2)
-
-let hash p =
- let c = List.map fst p in
- let m = List.map snd p in
- List.fold_left (fun h p -> h * 17 + P.hash p) (Hashtbl.hash m) c
-
-module Hashpol = Hashtbl.Make(
- struct
- type t = poly
- let equal = equal
- let hash = hash
- end)
-
-
-(*
- A pretty printer for polynomials, with Maple-like syntax.
- *)
-
-open Format
-
-let getvar lv i =
- try (List.nth lv i)
- with _ -> (List.fold_left (fun r x -> r^" "^x) "lv= " lv)
- ^" i="^(string_of_int i)
-
-let string_of_pol zeroP hdP tlP coefterm monterm string_of_coef
- dimmon string_of_exp lvar p =
-
- let rec string_of_mon m coefone =
- let s=ref [] in
- for i=1 to (dimmon m) do
- (match (string_of_exp m i) with
- "0" -> ()
- | "1" -> s:= (!s) @ [(getvar !lvar (i-1))]
- | e -> s:= (!s) @ [((getvar !lvar (i-1)) ^ "^" ^ e)]);
- done;
- (match !s with
- [] -> if coefone
- then "1"
- else ""
- | l -> if coefone
- then (String.concat "*" l)
- else ( "*" ^
- (String.concat "*" l)))
- and string_of_term t start = let a = coefterm t and m = monterm t in
- match (string_of_coef a) with
- "0" -> ""
- | "1" ->(match start with
- true -> string_of_mon m true
- |false -> ( "+ "^
- (string_of_mon m true)))
- | "-1" ->( "-" ^" "^(string_of_mon m true))
- | c -> if (String.get c 0)='-'
- then ( "- "^
- (String.sub c 1
- ((String.length c)-1))^
- (string_of_mon m false))
- else (match start with
- true -> ( c^(string_of_mon m false))
- |false -> ( "+ "^
- c^(string_of_mon m false)))
- and stringP p start =
- if (zeroP p)
- then (if start
- then ("0")
- else "")
- else ((string_of_term (hdP p) start)^
- " "^
- (stringP (tlP p) false))
- in
- (stringP p true)
-
-
-
-let print_pol zeroP hdP tlP coefterm monterm string_of_coef
- dimmon string_of_exp lvar p =
-
- let rec print_mon m coefone =
- let s=ref [] in
- for i=1 to (dimmon m) do
- (match (string_of_exp m i) with
- "0" -> ()
- | "1" -> s:= (!s) @ [(getvar !lvar (i-1))]
- | e -> s:= (!s) @ [((getvar !lvar (i-1)) ^ "^" ^ e)]);
- done;
- (match !s with
- [] -> if coefone
- then print_string "1"
- else ()
- | l -> if coefone
- then print_string (String.concat "*" l)
- else (print_string "*";
- print_string (String.concat "*" l)))
- and print_term t start = let a = coefterm t and m = monterm t in
- match (string_of_coef a) with
- "0" -> ()
- | "1" ->(match start with
- true -> print_mon m true
- |false -> (print_string "+ ";
- print_mon m true))
- | "-1" ->(print_string "-";print_space();print_mon m true)
- | c -> if (String.get c 0)='-'
- then (print_string "- ";
- print_string (String.sub c 1
- ((String.length c)-1));
- print_mon m false)
- else (match start with
- true -> (print_string c;print_mon m false)
- |false -> (print_string "+ ";
- print_string c;print_mon m false))
- and printP p start =
- if (zeroP p)
- then (if start
- then print_string("0")
- else ())
- else (print_term (hdP p) start;
- if start then open_hovbox 0;
- print_space();
- print_cut();
- printP (tlP p) false)
- in open_hovbox 3;
- printP p true;
- print_flush()
-
-
-let name_var= ref []
-
-let stringP = string_of_pol
- (fun p -> match p with [] -> true | _ -> false)
- (fun p -> match p with (t::p) -> t |_ -> failwith "print_pol dans dansideal")
- (fun p -> match p with (t::p) -> p |_ -> failwith "print_pol dans dansideal")
- (fun (a,m) -> a)
- (fun (a,m) -> m)
- string_of_coef
- (fun m -> (Array.length m)-1)
- (fun m i -> (string_of_int (m.(i))))
- name_var
-
-
-let stringPcut p =
- if (List.length p)>20
- then (stringP [List.hd p])^" + "^(string_of_int (List.length p))^" termes"
- else stringP p
-
-let rec lstringP l =
- match l with
- [] -> ""
- |p::l -> (stringP p)^("\n")^(lstringP l)
-
-let printP = print_pol
- (fun p -> match p with [] -> true | _ -> false)
- (fun p -> match p with (t::p) -> t |_ -> failwith "print_pol dans dansideal")
- (fun p -> match p with (t::p) -> p |_ -> failwith "print_pol dans dansideal")
- (fun (a,m) -> a)
- (fun (a,m) -> m)
- string_of_coef
- (fun m -> (Array.length m)-1)
- (fun m i -> (string_of_int (m.(i))))
- name_var
-
-
-let rec lprintP l =
- match l with
- [] -> ()
- |p::l -> printP p;print_string "\n"; lprintP l
-
-
-(* Operations
- *)
-
-let zeroP = []
-
-(* Retourne un polynome constant à d variables *)
-let polconst d c =
- let m = Array.create (d+1) 0 in
- let m = set_deg d m in
- [(c,m)]
-
-
-
-(* somme de polynomes= liste de couples (int,monomes) *)
-let plusP d p q =
- let rec plusP p q =
- match p with
- [] -> q
- |t::p' ->
- match q with
- [] -> p
- |t'::q' ->
- match compare_mon d (snd t) (snd t') with
- 1 -> t::(plusP p' q)
- |(-1) -> t'::(plusP p q')
- |_ -> let c=P.plusP (fst t) (fst t') in
- match P.equal c coef0 with
- true -> (plusP p' q')
- |false -> (c,(snd t))::(plusP p' q')
- in plusP p q
-
-
-(* multiplication d'un polynome par (a,monome) *)
-let mult_t_pol d a m p =
- let rec mult_t_pol p =
- match p with
- [] -> []
- |(b,m')::p -> ((P.multP a b),(mult_mon d m m'))::(mult_t_pol p)
- in mult_t_pol p
-
-
-(* Retourne un polynome de l dont le monome de tete divise m, [] si rien *)
-let rec selectdiv d m l =
- match l with
- [] -> []
- | q::r ->
- let m'= snd (List.hd q) in
- if div_mon_test d m m' then q else selectdiv d m r
-
-
-(* Retourne un polynome générateur 'i' à d variables *)
-let gen d i =
- let m = Array.create (d+1) 0 in
- m.(i) <- 1;
- let m = set_deg d m in
- [(coef1,m)]
-
-
-
-(* oppose d'un polynome *)
-let oppP p =
- let rec oppP p =
- match p with
- [] -> []
- |(b,m')::p -> ((P.oppP b),m')::(oppP p)
- in oppP p
-
-
-(* multiplication d'un polynome par un coefficient par 'a' *)
-let emultP a p =
- let rec emultP p =
- match p with
- [] -> []
- |(b,m')::p -> ((P.multP a b),m')::(emultP p)
- in emultP p
-
-
-let multP d p q =
- let rec aux p =
- match p with
- [] -> []
- |(a,m)::p' -> plusP d (mult_t_pol d a m q) (aux p')
- in aux p
-
-
-let puisP d p n=
- let rec puisP n =
- match n with
- 0 -> [coef1, Array.create (d+1) 0]
- | 1 -> p
- |_ -> multP d p (puisP (n-1))
- in puisP n
-
-let pgcdpos a b = P.pgcdP a b
-let pgcd1 p q =
- if P.equal p coef1 || P.equal p coefm1 then p else P.pgcdP p q
-
-let rec contentP p =
- match p with
- |[] -> coef1
- |[a,m] -> a
- |(a,m)::p1 -> pgcd1 a (contentP p1)
-
-let contentPlist lp =
- match lp with
- |[] -> coef1
- |p::l1 -> List.fold_left (fun r q -> pgcd1 r (contentP q)) (contentP p) l1
-
-(***********************************************************************
- Division de polynomes
- *)
-
-let hmon = (Hashtbl.create 1000 : (mon,poly) Hashtbl.t)
-
-let find_hmon m =
- if !use_hmon
- then Hashtbl.find hmon m
- else raise Not_found
-
-let add_hmon m q =
- if !use_hmon then Hashtbl.add hmon m q
-
-let selectdiv_cache d m l =
- try find_hmon m
- with Not_found ->
- match selectdiv d m l with
- [] -> []
- | q -> add_hmon m q; q
-
-let div_pol d p q a b m =
-(* info ".";*)
- plusP d (emultP a p) (mult_t_pol d b m q)
-
-(* si false on ne divise que le monome de tete *)
-let reduire_les_queues = false
-
-(* reste de la division de p par les polynomes de l
- rend le reste et le coefficient multiplicateur *)
-
-let reduce2 d p l =
- let rec reduce p =
- match p with
- [] -> (coef1,[])
- | (a,m)::p' ->
- let q = selectdiv_cache d m l in
- (match q with
- [] ->
- if reduire_les_queues
- then
- let (c,r)=(reduce p') in
- (c,((P.multP a c,m)::r))
- else (coef1,p)
- |(b,m')::q' ->
- let c=(pgcdpos a b) in
- let a'= (P.divP b c) in
- let b'=(P.oppP (P.divP a c)) in
- let (e,r)=reduce (div_pol d p' q' a' b'
- (div_mon d m m')) in
- (P.multP a' e,r)) in
- reduce p
-
-(* trace des divisions *)
-
-(* liste des polynomes de depart *)
-let poldep = ref []
-let poldepcontent = ref []
-
-
-module HashPolPair = Hashtbl.Make
- (struct
- type t = poly * poly
- let equal (p,q) (p',q') = equal p p' && equal q q'
- let hash (p,q) =
- let c = List.map fst p @ List.map fst q in
- let m = List.map snd p @ List.map snd q in
- List.fold_left (fun h p -> h * 17 + P.hash p) (Hashtbl.hash m) c
- end)
-
-(* table des coefficients des polynomes en fonction des polynomes de depart *)
-let coefpoldep = HashPolPair.create 51
-
-(* coefficient de q dans l expression de p = sum_i c_i*q_i *)
-let coefpoldep_find p q =
- try (HashPolPair.find coefpoldep (p,q))
- with _ -> []
-
-let coefpoldep_set p q c =
- HashPolPair.add coefpoldep (p,q) c
-
-let initcoefpoldep d lp =
- poldep:=lp;
- poldepcontent:= List.map contentP (!poldep);
- List.iter
- (fun p -> coefpoldep_set p p (polconst d coef1))
- lp
-
-(* garde la trace dans coefpoldep
- divise sans pseudodivisions *)
-
-let reduce2_trace d p l lcp =
- (* rend (lq,r), ou r = p + sum(lq) *)
- let rec reduce p =
- match p with
- [] -> ([],[])
- |t::p' -> let (a,m)=t in
- let q =
- (try Hashtbl.find hmon m
- with Not_found ->
- let q = selectdiv d m l in
- match q with
- t'::q' -> (Hashtbl.add hmon m q;q)
- |[] -> q) in
- match q with
- [] ->
- if reduire_les_queues
- then
- let (lq,r)=(reduce p') in
- (lq,((a,m)::r))
- else ([],p)
- |(b,m')::q' ->
- let b' = P.oppP (P.divP a b) in
- let m''= div_mon d m m' in
- let p1=plusP d p' (mult_t_pol d b' m'' q') in
- let (lq,r)=reduce p1 in
- ((b',m'',q)::lq, r)
- in let (lq,r) = reduce p in
- (*info "reduce2_trace:\n";
- List.iter
- (fun (a,m,s) ->
- let x = mult_t_pol d a m s in
- info ((stringP x)^"\n"))
- lq;
- info "ok\n";*)
- (List.map2
- (fun c0 q ->
- let c =
- List.fold_left
- (fun x (a,m,s) ->
- if equal s q
- then
- plusP d x (mult_t_pol d a m (polconst d coef1))
- else x)
- c0
- lq in
- c)
- lcp
- !poldep,
- r)
-
-(***********************************************************************
- Algorithme de Janet (V.P.Gerdt Involutive algorithms...)
-*)
-
-(***********************************
- deuxieme version, qui elimine des poly inutiles
-*)
-let homogeneous = ref false
-let pol_courant = ref []
-
-
-type pol3 =
- {pol : poly;
- anc : poly;
- nmp : mon}
-
-let fst_mon q = snd (List.hd q.pol)
-let lm_anc q = snd (List.hd q.anc)
-
-let pol_to_pol3 d p =
- {pol = p; anc = p; nmp = Array.create (d+1) 0}
-
-let is_multiplicative u s i =
- if i=1
- then List.for_all (fun q -> (fst_mon q).(1) <= u.(1)) s
- else
- List.for_all
- (fun q ->
- let v = fst_mon q in
- let res = ref true in
- let j = ref 1 in
- while !j < i && !res do
- res:= v.(!j) = u.(!j);
- j:= !j + 1;
- done;
- if !res
- then v.(i) <= u.(i)
- else true)
- s
-
-let is_multiplicative_rev d u s i =
- if i=1
- then
- List.for_all
- (fun q -> (fst_mon q).(d+1-1) <= u.(d+1-1))
- s
- else
- List.for_all
- (fun q ->
- let v = fst_mon q in
- let res = ref true in
- let j = ref 1 in
- while !j < i && !res do
- res:= v.(d+1- !j) = u.(d+1- !j);
- j:= !j + 1;
- done;
- if !res
- then v.(d+1-i) <= u.(d+1-i)
- else true)
- s
-
-let monom_multiplicative d u s =
- let m = Array.create (d+1) 0 in
- for i=1 to d do
- if is_multiplicative u s i
- then m.(i)<- 1;
- done;
- m
-
-(* mu monome des variables multiplicative de u *)
-let janet_div_mon d u mu v =
- let res = ref true in
- let i = ref 1 in
- while !i <= d && !res do
- if mu.(!i) = 0
- then res := u.(!i) = v.(!i)
- else res := u.(!i) <= v.(!i);
- i:= !i + 1;
- done;
- !res
-
-let find_multiplicative p mg =
- try Hashpol.find mg p.pol
- with Not_found -> (info "\nPROBLEME DANS LA TABLE DES VAR MULT";
- info (stringPcut p.pol);
- failwith "aie")
-
-(* mg hashtable de p -> monome_multiplicatif de g *)
-
-let hashtbl_reductor = ref (Hashtbl.create 51 : (mon,pol3) Hashtbl.t)
-
-(* suppose que la table a été initialisée *)
-let find_reductor d v lt mt =
- try Hashtbl.find !hashtbl_reductor v
- with Not_found ->
- let r =
- List.find
- (fun q ->
- let u = fst_mon q in
- let mu = find_multiplicative q mt in
- janet_div_mon d u mu v
- )
- lt in
- Hashtbl.add !hashtbl_reductor v r;
- r
-
-let normal_form d p g mg onlyhead =
- let rec aux = function
- [] -> (coef1,p)
- | (a,v)::p' ->
- (try
- let q = find_reductor d v g mg in
- let (b,u),q' = of_cons 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' = of_cons 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' = of_cons !h in
- let (b,u),q' = of_cons (!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 = of_cons !ls in
- ls := ls1;
- if !homogeneous && p.pol<>[] && deg_hom p.pol > deg_hom !pol_courant
- then info "h"
- else
- match head_normal_form d p lt mt with
- (_,v)::_ as h ->
- if fst_mon p <> v
- then lq := (!lq) @ [{pol = h; anc = h; nmp = Array.create (d+1) 0}]
- else lq := (!lq) @ [p]
- | [] ->
- (* info "*";*)
- if fst_mon p = lm_anc p
- then ls := List.filter (fun q -> not (equal q.anc p.anc)) !ls
- done;
- (*info ("temps de head_reduce: "
- ^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));*)
- !lq
-
-let choose_irreductible d lf =
- List.hd lf
-(* bien plus lent
- (List.sort (fun p q -> compare_mon d (fst_mon p.pol) (fst_mon q.pol)) lf)
-*)
-
-
-let hashtbl_multiplicative d lf =
- let mg = Hashpol.create 51 in
- hashtbl_reductor := Hashtbl.create 51;
- List.iter
- (fun g ->
- let (_,u) = List.hd g.pol in
- Hashpol.add mg g.pol (monom_multiplicative d u lf))
- lf;
- (*info ("temps de hashtbl_multiplicative: "
- ^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));*)
- mg
-
-let list_diff l x =
- List.filter (fun y -> y <> x) l
-
-let janet2 d lf p0 =
- hashtbl_reductor := Hashtbl.create 51;
- let t1 = Unix.gettimeofday() in
- info ("------------- Janet2 ------------------\n");
- pol_courant := p0;
- let r = ref p0 in
- let lf = List.map (pol_to_pol3 d) lf in
- let f = choose_irreductible d lf in
- let lt = ref [f] in
- let mt = ref (hashtbl_multiplicative d !lt) in
- let lq = ref (list_diff lf f) in
- lq := head_reduce d !lq !lt !mt;
- r := (* lent head_janet_normal_form d !r !lt !mt ; *)
- reduce_rem d !r !lt !lq ;
- info ("reste: "^(stringPcut !r)^"\n");
- while !lq <> [] && !r <> [] do
- let p = choose_irreductible d !lq in
- lq := list_diff !lq p;
- if p.pol = p.anc
- then ( (* on enleve de lt les pol divisibles par p et on les met dans lq *)
- let m = fst_mon p in
- let lt1 = !lt in
- List.iter
- (fun q ->
- let m'= fst_mon q in
- if div_strict d m m'
- then (
- lq := (!lq) @ [q];
- lt := list_diff !lt q))
- lt1;
- mt := hashtbl_multiplicative d !lt;
- );
- let h = tail_normal_form d p.pol !lt !mt in
- if h = []
- then info "************ reduction a zero, ne devrait pas arriver *\n"
- else (
- lt := (!lt) @ [{pol = h; anc = p.anc; nmp = Array.copy p.nmp}];
- info ("nouveau polynome: "^(stringPcut h)^"\n");
- mt := hashtbl_multiplicative d !lt;
- r := (* lent head_janet_normal_form d !r !lt !mt ; *)
- reduce_rem d !r !lt !lq ;
- info ("reste: "^(stringPcut !r)^"\n");
- if !r <> []
- then (
- List.iter
- (fun q ->
- let mq = find_multiplicative q !mt in
- for i=1 to d do
- if mq.(i) = 1
- then q.nmp.(i)<- 0
- else
- if q.nmp.(i) = 0
- then (
- (* info "+";*)
- lq := (!lq) @
- [{pol = multP d (gen d i) q.pol;
- anc = q.anc;
- nmp = Array.create (d+1) 0}];
- q.nmp.(i)<-1;
- );
- done;
- )
- !lt;
- lq := head_reduce d !lq !lt !mt;
- info ("["^(string_of_int (List.length !lt))^";"
- ^(string_of_int (List.length !lq))^"]");
- ));
- done;
- info ("--- Janet2 donne:\n");
- (* List.iter (fun p -> info ("polynome: "^(stringPcut p.pol)^"\n")) !lt; *)
- info ("reste: "^(stringPcut !r)^"\n");
- info ("--- fin Janet2\n");
- info ("temps: "^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));
- List. map (fun q -> q.pol) !lt
-
-(**********************************************************************
- version 3 *)
-
-let head_normal_form3 d p lt mt =
- let h = ref (p.pol) in
- let res =
- try (
- let v = snd(List.hd !h) in
- let g = ref (find_reductor d v lt mt) in
- if snd(List.hd !h) <> lm_anc p && criteria d p !g lt
- then ((* info "=";*) [])
- else (
- while !h <> [] && (!g).pol <> [] do
- let (a,v),p' = of_cons !h in
- let (b,u),q' = of_cons (!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 = of_cons 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 = of_cons !lq in
- lq := lq1;
-(*
- if p.pol = p.anc
- then ( (* on enleve de lt les pol divisibles par p et on les met dans lq *)
- let m = fst_mon (p.pol) in
- let lt1 = !lt in
- List.iter
- (fun q ->
- let m'= fst_mon (q.pol) in
- if div_strict d m m'
- then (
- lq := (!lq) @ [q];
- lt := list_diff !lt q))
- lt1;
- mt := hashtbl_multiplicative d !lt;
- );
-*)
- let h = head_normal_form3 d p !lt !mt in
- if h = []
- then (
- info "*";
- if fst_mon p = lm_anc p
- then
- lq := List.filter (fun q -> not (equal q.anc p.anc)) !lq;
- )
- else (
- let q =
- if fst_mon p <> snd(List.hd h)
- then {pol = h; anc = h; nmp = Array.create (d+1) 0}
- else {pol = h; anc = p.anc; nmp = Array.copy p.nmp}
- in
- lt:= (!lt) @ [q];
- info ("nouveau polynome: "^(stringPcut q.pol)^"\n");
- mt := hashtbl_multiplicative d !lt;
- r := reduce_rem d !r !lt !lq ; (* janet_normal_form d !r !lt !mt ;*)
- info ("reste: "^(stringPcut !r)^"\n");
- if !r <> []
- then (
- List.iter
- (fun q ->
- let mq = find_multiplicative q !mt in
- for i=1 to d do
- if mq.(i) = 1
- then q.nmp.(i)<- 0
- else
- if q.nmp.(i) = 0
- then (
- (* info "+";*)
- lq := (!lq) @
- [{pol = multP d (gen d i) q.pol;
- anc = q.anc;
- nmp = Array.create (d+1) 0}];
- q.nmp.(i)<-1;
- );
- done;
- )
- !lt;
- info ("["^(string_of_int (List.length !lt))^";"
- ^(string_of_int (List.length !lq))^"]");
- ));
- done;
- info ("--- Janet3 donne:\n");
- (* List.iter (fun p -> info ("polynome: "^(stringPcut p.pol)^"\n")) !lt; *)
- info ("reste: "^(stringPcut !r)^"\n");
- info ("--- fin Janet3\n");
- info ("temps: "^(Format.sprintf "@[%10.3f@]s\n" ((Unix.gettimeofday ())-.t1)));
- List. map (fun q -> q.pol) !lt
-
-
-(***********************************************************************
- Completion
- *)
-
-let sugar_flag = ref true
-
-let sugartbl = (Hashpol.create 1000 : int Hashpol.t)
-
-let compute_sugar p =
- List.fold_left (fun s (a,m) -> max s m.(0)) 0 p
-
-let sugar p =
- try Hashpol.find sugartbl p
- with Not_found ->
- let s = compute_sugar p in
- Hashpol.add sugartbl p s;
- s
-
-let spol d p q=
- let m = snd (List.hd p) in
- let m'= snd (List.hd q) in
- let a = fst (List.hd p) in
- let b = fst (List.hd q) in
- let p'= List.tl p in
- let q'= List.tl q in
- let c=(pgcdpos a b) in
- let m''=(ppcm_mon d m m') in
- let m1 = div_mon d m'' m in
- let m2 = div_mon d m'' m' in
- let fsp p' q' =
- plusP d
- (mult_t_pol d (P.divP b c) m1 p')
- (mult_t_pol d (P.oppP (P.divP a c)) m2 q') in
- let sp = fsp p' q' in
- coefpoldep_set sp p (fsp (polconst d coef1) []);
- coefpoldep_set sp q (fsp [] (polconst d coef1));
- Hashpol.add sugartbl sp
- (max (m1.(0) + (sugar p)) (m2.(0) + (sugar q)));
- sp
-
-let etrangers d p p'=
- let m = snd (List.hd p) in
- let m'= snd (List.hd p') in
- let res=ref true in
- let i=ref 1 in
- while (!res) && (!i<=d) do
- res:= (m.(!i) = 0) || (m'.(!i)=0);
- i:=!i+1;
- done;
- !res
-
-
-(* teste si le monome dominant de p''
- divise le ppcm des monomes dominants de p et p' *)
-
-let div_ppcm d p p' p'' =
- let m = snd (List.hd p) in
- let m'= snd (List.hd p') in
- let m''= snd (List.hd p'') in
- let res=ref true in
- let i=ref 1 in
- while (!res) && (!i<=d) do
- res:= ((max m.(!i) m'.(!i)) >= m''.(!i));
- i:=!i+1;
- done;
- !res
-
-(***********************************************************************
- Code extrait de la preuve de L.Thery en Coq
- *)
-type 'poly cpRes =
- Keep of ('poly list)
- | DontKeep of ('poly list)
-
-let addRes i = function
- Keep h'0 -> Keep (i::h'0)
- | DontKeep h'0 -> DontKeep (i::h'0)
-
-let rec slice d i a = function
- [] -> if etrangers d i a then DontKeep [] else Keep []
- | b::q1 ->
- if div_ppcm d i a b then DontKeep (b::q1)
- else if div_ppcm d i b a then slice d i a q1
- else addRes b (slice d i a q1)
-
-let rec addS x l = l @[x]
-
-let addSugar x l =
- if !sugar_flag
- then
- let sx = sugar x in
- let rec insere l =
- match l with
- | [] -> [x]
- | y::l1 ->
- if sx <= (sugar y)
- then x::l
- else y::(insere l1)
- in insere l
- else addS x l
-
-(* ajoute les spolynomes de i avec la liste de polynomes aP,
- a la liste q *)
-
-let rec genPcPf d i aP q =
- match aP with
- [] -> q
- | a::l1 ->
- (match slice d i a l1 with
- Keep l2 -> addSugar (spol d i a) (genPcPf d i l2 q)
- | DontKeep l2 -> genPcPf d i l2 q)
-
-let rec genOCPf d = function
- [] -> []
- | a::l -> genPcPf d a l (genOCPf d l)
-
-let step = ref 0
-
-let infobuch p q =
- if !step = 0
- then (info ("[" ^ (string_of_int (List.length p))
- ^ "," ^ (string_of_int (List.length q))
- ^ "]"))
-
-(***********************************************************************)
-(* dans lp les nouveaux sont en queue *)
-
-let coef_courant = ref coef1
-
-
-type certificate =
- { coef : coef; power : int;
- gb_comb : poly list list; last_comb : poly list }
-
-let test_dans_ideal d p lp lp0 =
- let (c,r) = reduce2 d !pol_courant lp in
- info ("reste: "^(stringPcut r)^"\n");
- coef_courant:= P.multP !coef_courant c;
- pol_courant:= r;
- if r=[]
- then (info "polynome reduit a 0\n";
- let lcp = List.map (fun q -> []) !poldep in
- let c = !coef_courant in
- let (lcq,r) = reduce2_trace d (emultP c p) lp lcp in
- info ("r: "^(stringP r)^"\n");
- let res=ref (emultP c p) in
- List.iter2
- (fun cq q -> res:=plusP d (!res) (multP d cq q);
- )
- lcq !poldep;
- info ("verif somme: "^(stringP (!res))^"\n");
- info ("coefficient: "^(stringP (polconst d c))^"\n");
- let rec aux lp =
- match lp with
- |[] -> []
- |p::lp ->
- (List.map
- (fun q -> coefpoldep_find p q)
- lp)::(aux lp)
- in
- let liste_polynomes_de_depart = List.rev lp0 in
- let polynome_a_tester = p in
- let liste_des_coefficients_intermediaires =
- (let lci = List.rev (aux (List.rev lp)) in
- let lci = ref lci (* (List.map List.rev lci) *) in
- List.iter (fun x -> lci := List.tl (!lci)) lp0;
- !lci) in
- let liste_des_coefficients =
- List.map
- (fun cq -> emultP coefm1 cq)
- (List.rev lcq) in
- (liste_polynomes_de_depart,
- polynome_a_tester,
- {coef=c;
- power=1;
- gb_comb = liste_des_coefficients_intermediaires;
- last_comb = liste_des_coefficients}))
- else ((*info "polynome non reduit a 0\n";
- info ("\nreste: "^(stringPcut r)^"\n");*)
- raise NotInIdeal)
-
-let divide_rem_with_critical_pair = ref false
-
-let choix_spol d p l =
- if !divide_rem_with_critical_pair
- then (
- let (_,m) = List.hd p in
- try (
- let q =
- List.find
- (fun q ->
- let (_,m') = List.hd q in
- div_mon_test d m m')
- l in
- q::(list_diff l q))
- with _ -> l)
- else l
-
-let pbuchf d pq p lp0=
- info "calcul de la base de Groebner\n";
- step:=0;
- Hashtbl.clear hmon;
- let rec pbuchf lp lpc =
- infobuch lp lpc;
-(* step:=(!step+1)mod 10;*)
- match lpc with
- [] -> test_dans_ideal d p lp lp0
- | _ ->
- let a,lpc2 = of_cons (choix_spol d !pol_courant lpc) in
- if !homogeneous && a<>[] && deg_hom a > deg_hom !pol_courant
- then (info "h";pbuchf lp lpc2)
- else
- let sa = sugar a in
- let (ca,a0)= reduce2 d a lp in
- match a0 with
- [] -> info "0";pbuchf lp lpc2
- | _ ->
- let a1 = emultP ca a in
- List.iter
- (fun q ->
- coefpoldep_set a1 q (emultP ca (coefpoldep_find a q)))
- !poldep;
- let (lca,a0) = reduce2_trace d a1 lp
- (List.map (fun q -> coefpoldep_find a1 q) !poldep) in
- List.iter2 (fun c q -> coefpoldep_set a0 q c) lca !poldep;
- info ("\nnouveau polynome: "^(stringPcut a0)^"\n");
- let ct = coef1 (* contentP a0 *) in
- (*info ("contenu: "^(string_of_coef ct)^"\n");*)
- Hashpol.add sugartbl a0 sa;
- poldep:=addS a0 lp;
- poldepcontent:=addS ct (!poldepcontent);
- try test_dans_ideal d p (addS a0 lp) lp0
- with NotInIdeal -> pbuchf (addS a0 lp) (genPcPf d a0 lp lpc2)
- in pbuchf (fst pq) (snd pq)
-
-let is_homogeneous p =
- match p with
- | [] -> true
- | (a,m)::p1 -> let d = m.(0) in
- List.for_all (fun (b,m') -> m'.(0)=d) p1
-
-(* rend
- c
- lp = [pn;...;p1]
- p
- lci = [[a(n+1,n);...;a(n+1,1)];
- [a(n+2,n+1);...;a(n+2,1)];
- ...
- [a(n+m,n+m-1);...;a(n+m,1)]]
- lc = [qn+m; ... q1]
-
- tels que
- c*p = sum qi*pi
- ou pn+k = a(n+k,n+k-1)*pn+k-1 + ... + a(n+k,1)* p1
- *)
-
-let in_ideal d lp p =
- homogeneous := List.for_all is_homogeneous (p::lp);
- if !homogeneous then info "polynomes homogenes\n";
- (* janet2 d lp p;*)
- info ("p: "^stringPcut p^"\n");
- info ("lp:\n"^List.fold_left (fun r p -> r^stringPcut p^"\n") "" lp);
- initcoefpoldep d lp;
- coef_courant:=coef1;
- pol_courant:=p;
- try test_dans_ideal d p lp lp
- with NotInIdeal -> pbuchf d (lp, genOCPf d lp) p lp
-
-end
diff --git a/plugins/groebner/ideal.mli b/plugins/groebner/ideal.mli
deleted file mode 100644
index 9c8f43d7d..000000000
--- a/plugins/groebner/ideal.mli
+++ /dev/null
@@ -1,80 +0,0 @@
-(************************************************************************)
-(* 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 *)
-(************************************************************************)
-
-exception NotInIdeal
-val lexico : bool ref
-val use_hmon : bool ref
-
-module type S = sig
-
-(* Monomials *)
-type mon = int array
-
-val mult_mon : int -> mon -> mon -> mon
-val deg : mon -> int
-val compare_mon : int -> mon -> mon -> int
-val div_mon : int -> mon -> mon -> mon
-val div_mon_test : int -> mon -> mon -> bool
-val ppcm_mon : int -> mon -> mon -> mon
-
-(* Polynomials *)
-
-type deg = int
-type coef
-type poly
-val repr : poly -> (coef * mon) list
-val polconst : deg -> coef -> poly
-val zeroP : poly
-val gen : deg -> int -> poly
-
-val equal : poly -> poly -> bool
-val name_var : string list ref
-val getvar : string list -> int -> string
-val lstringP : poly list -> string
-val printP : poly -> unit
-val lprintP : poly list -> unit
-
-val div_pol_coef : poly -> coef -> poly
-val plusP : deg -> poly -> poly -> poly
-val mult_t_pol : deg -> coef -> mon -> poly -> poly
-val selectdiv : deg -> mon -> poly list -> poly
-val oppP : poly -> poly
-val emultP : coef -> poly -> poly
-val multP : deg -> poly -> poly -> poly
-val puisP : deg -> poly -> int -> poly
-val contentP : poly -> coef
-val contentPlist : poly list -> coef
-val pgcdpos : coef -> coef -> coef
-val div_pol : deg -> poly -> poly -> coef -> coef -> mon -> poly
-val reduce2 : deg -> poly -> poly list -> coef * poly
-
-val poldepcontent : coef list ref
-val coefpoldep_find : poly -> poly -> poly
-val coefpoldep_set : poly -> poly -> poly -> unit
-val initcoefpoldep : deg -> poly list -> unit
-val reduce2_trace : deg -> poly -> poly list -> poly list -> poly list * poly
-val spol : deg -> poly -> poly -> poly
-val etrangers : deg -> poly -> poly -> bool
-val div_ppcm : deg -> poly -> poly -> poly -> bool
-
-val genPcPf : deg -> poly -> poly list -> poly list -> poly list
-val genOCPf : deg -> poly list -> poly list
-
-val is_homogeneous : poly -> bool
-
-type certificate =
- { coef : coef; power : int;
- gb_comb : poly list list; last_comb : poly list }
-val test_dans_ideal : deg -> poly -> poly list -> poly list ->
- poly list * poly * certificate
-
-val in_ideal : deg -> poly list -> poly -> poly list * poly * certificate
-
-end
-
-module Make (P:Polynom.S) : S with type coef = P.t
diff --git a/plugins/groebner/polynom.ml b/plugins/groebner/polynom.ml
deleted file mode 100644
index 0a9c3e270..000000000
--- a/plugins/groebner/polynom.ml
+++ /dev/null
@@ -1,1051 +0,0 @@
-(************************************************************************)
-(* 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 *)
-(************************************************************************)
-
-(*
- Polynomes récursifs: Z[x1]...[xn].
-*)
-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
-
-(***********************************************************************
- 2. Le type des polynômes, 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 (* 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:
- - les variables sont des entiers strictement positifs.
- - les coefficients d'un polynome en x ne font intervenir que des variables < x.
- - pas de coefficient nul en tête.
- - pas de Prec(x,a) ou a n'a qu'un element (constant en x).
-*)
-
-(* Polynômes constant *)
-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
- Pint _ -> 0
- |Prec(x,_) -> x
-
-
-(* p n'est pas forcément normalisé *)
-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
-
-
-(* variable max d'une liste de polynômes *)
-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 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
-
-(* vire les zéros de tête d'un polynôme non normalisé, dont les coefficients
- sont supposés normalisés.
- si constant, rend le coef constant.
-*)
-
-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))
-
-
-(* degré en la variable v du polynome p, v >= max var de p *)
-let rec deg v p =
- match p with
- Prec(x,p1) when x=v -> Array.length p1 -1
- |_ -> 0
-
-
-(* degré total *)
-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
-
-
-(* copie le polynome *)
-let rec copyP p =
- match p with
- Pint i -> Pint i
- |Prec(x,q) -> Prec(x,Array.map copyP q)
-
-
-(* coefficient de degre i en v, v >= max var de 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
-
-
-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 x<y then (let q2=Array.map copyP q1 in
- q2.(0)<- plusP p q1.(0);
- Prec (y,q2))
- else if x>y 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
-
-
-(* contenu entier positif *)
-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)
-
-
-(* divise tous les coefficients de p par l'entier a*)
-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)
-
-(* divise p par son contenu entier positif. *)
-let vire_contenu p =
- let c = content p in
- if C.equal c coef0 then p else div_int p c
-
-
-(* liste triee des variables impliquees dans un poly *)
-let rec vars=function
- Pint _->[]
- | Prec (x,l)->(List.flatten ([x]::(List.map vars (Array.to_list l))))
-
-
-(* multiplie p par 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))
-
-
-(* produit de 2 polynomes *)
-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<y
- then (let q2 = Array.map (fun z-> 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 par rapport a la 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
-
-
-(* opposé de p *)
-let rec oppP p =
- match p with
- Pint a -> Pint (C.opp a)
- |Prec(x,p1) -> Prec(x,Array.map oppP p1)
-
-
-(* différence de deux polynômes. *)
-let moinsP p q=plusP p (oppP q)
-
-let rec puisP p n = match n with
- 0 -> cf1
- |_ -> (multP p (puisP p (n-1)))
-
-
-(* notations infixes...*)
-(*let (++) a b = plusP a b
-*)
-let (@@) a b = multP a b
-
-let (--) a b = moinsP a b
-
-let (^^) a b = puisP a b
-
-
-(* coefficient dominant de p, v>= max_var p *)
-
-let coefDom v p= coef v (deg v p) p
-
-
-let coefConst v p = coef v 0 p
-
-(* queue d'un polynôme *)
-let remP v p =
- moinsP p (multP (coefDom v p) (puisP (x v) (deg v p)))
-
-
-(* premier coef entier de 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)
-
-
-(* divise par le contenu, et rend positif le premier coefficient entier *)
-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
-
-
-(*coef constant d'un polynome normalise*)
-let rec coef_constant p =
- match p with
- Pint a->a
- |Prec(_,q)->coef_constant q.(0)
-
-
-(***********************************************************************
- 3. Affichage des polynômes.
-*)
-
-(* si univ=false, on utilise x,y,z,a,b,c,d... comme noms de variables,
- sinon, x1,x2,...
-*)
-let univ=ref true
-
-(* joli jusqu'a trois variables -- sinon changer le 'w' *)
-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)
-
-(* #install_printer printP *)
-
-(***********************************************************************
- 4. Division exacte de polynômes.
-*)
-
-(* rend (s,r) tel que 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<m
- then continue:=false
- else (
- let a = coefDom x !r in
- let p1 = remP x !r in (* r = a*x^n+p1 *)
- let c = div_pol a b (x-1) in (* a = c*b *)
- let s1 = c @@ ((monome x (n-m))) in
- s:= plusP (!s) s1;
- r:= p1 -- (s1 @@ q1);
- )
- done;
- (!s,!r)
-
-(* 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 equal r cf0
- then s
- else failwith ("div_pol:\n"
- ^"p:"^(to_string p)^"\n"
- ^"q:"^(to_string q)^"\n"
- ^"r:"^(to_string r)^"\n"
- ^"x:"^(string_of_int x)^"\n"
- )
-
-
-(* test de division exacte de p par q mais constantes rationnels
- à vérifier *)
-let divP p q=
- let x = max (max_var_pol p) (max_var_pol q) in
- div_pol p q x
-
-(* test de division exacte de p par q mais constantes rationnels
- à vérifier *)
-let div_pol_rat p q=
- let x = max (max_var_pol p) (max_var_pol q) in
- try (let s = div_pol (multP p (puisP (Pint(coef_int_tete q))
- (1+(deg x p) - (deg x q))))
- q x in
- (*degueulasse, mais c 'est pour enlever un warning *)
- if s==s then true else true)
- with _ -> false
-
-
-
-
-(***********************************************************************
- 5. Pseudo-division et pgcd par les sous-résultants.
-*)
-
-(* pseudo division :
- q = c*x^m+q1
- rend (r,c,d,s) tels que 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)
- )
-
-
-(* pgcd de polynômes par les sous-résultants *)
-
-
-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
- )
-
-(* Sous-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<d'
- then gcd_sub_res q p x
- else
- let delta = d-d' in
- let c' = coefDom x q in
- let r = snd (quo_rem_pol (((oppP c')^^(delta+1))@@p) (oppP q) x) in
- gcd_sub_res_rec q r (c'^^delta) c' d' x
-
-and gcd_sub_res_rec p q s c d x =
- if equal q cf0
- then p
- else (
- let d' = deg x q in
- let c' = coefDom x q in
- let delta = d-d' in
- let r = snd (quo_rem_pol (((oppP c')^^(delta+1))@@p) (oppP q) x) in
- let s'= lazard_power c' s delta x in
- gcd_sub_res_rec q (div_pol r (c @@ (s^^delta)) x) s' c' d' x
- )
-
-and lazard_power c s d x =
- let res = ref c in
- for i=1 to d-1 do
- res:= div_pol ((!res)@@c) s x;
- done;
- !res
-
-
-
-(***********************************************************************
- 6. Décomposition sans carré, factorisation.
-*)
-
-(*
- p = f1 f2^2 ... fn^r
- p/\p'= f2 f3^2...fn^(r-1)
- sans_carré(p)= p/p/\p '= f1 f2 ... fn
-*)
-
-let sans_carre p x =
- if (deg x p) <= 1 then p
- else
- let p' = deriv x p in
- div_pol p (pgcd_pol p p' x) x
-
-
-(* liste [f1;...;fn] *)
-let facteurs p x =
- let rec facteurs_rec p q =
- if (deg x p)=0 then []
- else
- let p2 = div_pol p q x in
- let q2 = sans_carre p2 x in
- (div_pol q q2 x)::(facteurs_rec p2 q2)
- in facteurs_rec p (sans_carre p x)
-
-
-(* liste [f1;f3;...] *)
-let facteurs_impairs p x =
- let lf = Array.of_list (facteurs p x) in
- let r = ref [] in
- Array.iteri (fun i f ->
- if ((i+1) mod 2)=1
- then r:=(!r)@[f])
- lf;
- !r
-
-
-(* décomposition sans carrés en toutes les variables *)
-
-
-let hcontentP = (Hashtbl.create 51 : (string,t) Hashtbl.t)
-
-let prcontentP () =
- Hashtbl.iter (fun s c ->
- prn (s^" -> "^(to_string c)))
- hcontentP
-
-
-let contentP =
- 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 *)
-
-(* fonction de hachage des polynômes *)
-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 = equal
- let hash = hash
- end)
-
-
-let memoP s memoire fonction x =
- try (let v = Hashpol.find memoire x in pr s;v)
- with _ -> (pr "#";
- let v = fonction x in
- Hashpol.add memoire x v;
- v)
-
-
-let hfactorise = Hashpol.create 51
-
-let prfactorise () =
- Hashpol.iter (fun p c ->
- prn ((to_string p)^" -> ");
- print_lpoly (List.flatten c))
- hfactorise
-
-let factorise =
- memoP "f" hfactorise
- (fun p ->
- let rec fact p x =
- if x=0
- then []
- else
- let c = contentP (p,x) in
- let q = div_pol p c x in
- (facteurs q x)::(fact c (x-1))
- in fact p (max_var_pol p))
-
-
-(* liste des facteurs sans carré non constants,
- avec coef entier de tête positif *)
-let facteurs2 p =
- List.map normc
- (List.filter (fun q -> deg_total q >0)
- (List.flatten (factorise (normc p))))
-
-
-(* produit des facteurs non constants d'une décomposition sans carré *)
-let pol_de_factorisation lf =
- let p = ref cf1 in
- List.iter (fun lq ->
- Array.iteri (fun i q -> if (deg_total q)>0 then p:=(!p)@@(q^^(i+1)))
- (Array.of_list lq))
- lf;
- !p
-
-
-let set_of_array_facteurs tf =
- let h = Hashpol.create 51 in
- Array.iter (fun lf ->
- List.iter (fun p -> if not (Hashpol.mem h p)
- then Hashpol.add h p true)
- lf)
- tf;
- let res = ref [] in
- Hashpol.iter (fun p _ -> res:=p::(!res)) h;
- !res
-
-
-(* Factorise un tableau de polynômes f, et rend:
- - un tableau p de facteurs (degré>0, contenu entier 1,
- coefficient de tête >0) obtenu par décomposition sans carrés
- puis par division mutuelle
- - un tableau l de couples (constante, listes d'indices l)
- tels que f.(i) = l.(i)_1*Produit(p.(j), j dans l.(i)_2)
-*)
-
-(* on donne le tableau des facteurs de chaque poly de f *)
-let factorise_tableauP2 f l1 =
- let x = max_var f in
- (* liste des facteurs sans carré des polynômes de f *)
- pr"<";
- let l1 = set_of_array_facteurs l1 in
- (* on les divise entre eux pour éventuellement trouver
- de nouveaux facteurs *)
- pr "~";
- let l1 = Sort.list (fun p q -> (deg_total p)<(deg_total q)) l1 in
- let l1 = Array.of_list (facteurs_liste (fun a b -> div_pol a b x)
- (fun p -> (deg_total p)<1)
- l1) in
- (* 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 -> equal p cf0)
- cf0
- f l1 in
- pr ">";
- res
-
-let factorise_tableauP f =
- factorise_tableauP2 f (Array.map facteurs2 f)
-
-
-(***********************************************************************
- 7. Pseudo-division avec reste de même signe,
- en utilisant des polynômes non nuls pour diviser le reste.
-*)
-
-(* polynôme pair et coefficients positifs *)
-let rec is_positif p =
-
- let res =
- match p with
- 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 (equal c cf0)
- then failwith "pas pair")
- p1;
- true)
- with Failure _ -> false)
- in
- res
-
-
-
-let is_negatif p = is_positif (oppP p)
-
-
-(* rend r tel que deg r < deg q et r a le signe de p en les racines de q.
- le coefficient dominant de q est non nul
- quand les polynômes de coef_non_nuls le sont.
- (rs,cs,ds,ss,crs,lpos,lpol)= pseudo_euclide coef_non_nuls vect.(s-1) res.(s-1) v
-*)
-let pseudo_euclide coef_non_nuls p q x =
- let (r,c,d,s) = pseudo_div p q x in
- (*
- c^d * p = s*q + r, c = coef dominant de q
- *)
- (* vérification de la pseudo-division:
- let verif = ((c^^d)@@p)--(plusP (s@@q) r) in
- 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");
- *)
-
- (* pour autoriser des c pas positifs, il faut modifier algo14 et preuve3*)
- let r = if d mod 2 = 1 then c@@r else r in
- let s = if d mod 2 = 1 then c@@s else s in
- 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 equal r cf0
- then ((*pr "reste nul"; *) (r,c,d,s,cf1,[],[]))
- else (
- let r1 = vire_contenu r in
- let cr = div_pol r r1 x in
- let r = ref r1 in
- (* r a maintenant le même signe que p en les racines de q.*)
- (* on tente de diviser r par les polynômes de coef_non_nuls *)
- let lf = ref [] in (* liste de (facteur, puissance) *)
- List.iter (fun f ->
- if (deg_total f)>0 && (max_var_pol f) < x
- then (
- let k = ref 0 in
- (try (while true do
- let rd = div_pol !r f x in
- (* verification de la division
- if not (equal cf0 ((!r)--(f@@rd)))
- then failwith "erreur dans la division";
- *)
- k:=(!k)+1;
- r:=rd;
- (*pr "+";*)
- done)
- with _ -> ());
- lf:=(f,!k)::(!lf)))
- coef_non_nuls;
- (* il faut éventuellement remultiplier pour garder le signe de r *)
- let lpos = ref [] in
- let lpol = ref [] in
- List.iter (fun (f,k) ->
- if k>0
- then (
- if (is_positif f)
- (* f est positif, tout va bien *)
- then lpos:=(f,k)::(!lpos)
- else if (is_negatif f)
- (* f est négatif *)
- then (if k mod 2 = 1
- (* k est impair *)
- then (lpos:=(oppP f,k)::(!lpos);
- r:=oppP (!r))
- else lpos:=(f,k)::(!lpos))
- (* on ne connaît pas le signe de f *)
- else if k mod 2 = 0
- (* k est pair, tout va bien *)
- then lpol:=(f,k)::(!lpol)
- (* k est impair *)
- else (lpol:=(f,k-1)::(!lpol);
- r:=multP (!r) f)))
- !lf;
- (*
- pr ((* "------reste de même signe: "
- ^(to_string c)
- ^" variable: "^(string_of_int x)
- ^" c:"^(string_of_int (deg_total c))
- ^" d:"^(string_of_int d)
- ^" deg(r)-deg(r0):"
- ^*)(string_of_int ((deg_total !r)-(deg_total r0))));
- *)
- (* lpos = liste de (f,k) ou f est non nul positif, et f^k divise r0
- lpol = liste de (f,k) ou f non nul, k est pair et f^k divise r0
- on c^d * p = s*q + r0
- avec d pair
- r0 = cr * r * PI_lpos f^k * PI_lpol g^k
- cr non nul positif
- *)
- (!r,c,d,s,cr,!lpos,!lpol))
-
-
-
-(* teste si la non-nullité des polynômes de lp entraîne celle de p:
- chacun des facteurs de la décomposition sans carrés de p
- divise un des polynômes de lp (dans Q[x1...xn]) *)
-
-let implique_non_nul lp p =
- if equal p cf0 then false
- else(
- pr "[";
- let lf = facteurs2 p in
- let r =(
- try (List.iter (fun f ->
- if (try (List.iter (fun q ->
- if div_pol_rat q f
- then failwith "divise")
- lp;
- true)
- with _ -> false)
- then failwith "raté")
- lf;
- true)
- with _ -> false)
- in pr "]";r)
-
-
-let ajoute_non_nul p lp =
- if (deg_total p) >0
- then(
- let p = normc p in
- let lf = facteurs2 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
deleted file mode 100644
index 79680262a..000000000
--- a/plugins/groebner/polynom.mli
+++ /dev/null
@@ -1,120 +0,0 @@
-(************************************************************************)
-(* 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 *)
-(************************************************************************)
-
-(* Building recursive polynom operations from a type of coefficients *)
-
-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
deleted file mode 100644
index 40644489b..000000000
--- a/plugins/groebner/utile.ml
+++ /dev/null
@@ -1,161 +0,0 @@
-(***********************************************************************
- Utilitaires.
-*)
-
-(* Impression *)
-
-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
-
-(**********************************************************************
- Listes
-*)
-
-(* appartenance à une liste , on donne l'égalité *)
-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)
-
-(* vire les repetitions d'une liste, on donne l'égalité *)
-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
-
-
-(***********************************************************************
- Un outil pour faire une mémo-fonction:
- fonction est la fonction(!)
- memoire est une référence au graphe déjà calculé
- (liste de couples, c'est une variable globale)
- egal est l'égalité sur les arguments
- valeur est une valeur possible de la fonction (sert uniquement pour le typage)
-*)
-
-let memo memoire egal valeur fonction x =
- let res = ref valeur in
- try (List.iter (fun (y,v) -> if egal y x
- then (res:=v;
- failwith "trouve"))
- !memoire;
- let v = fonction x in
- memoire:=(x,v)::(!memoire);
- v)
- with _ -> !res
-
-
-(* un autre plus efficace,
- utilisant une fonction intermediaire (utile si on n'a pas
- l'égalité = sur les arguments de fonction)
- s chaîne imprimée s'il n'y a pas calcul *)
-
-let memos s memoire print fonction x =
- try (let v = Hashtbl.find memoire (print x) in pr s;v)
- with _ -> (pr "#";
- let v = fonction x in
- Hashtbl.add memoire (print 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/groebner/utile.mli b/plugins/groebner/utile.mli
deleted file mode 100644
index 209258dd6..000000000
--- a/plugins/groebner/utile.mli
+++ /dev/null
@@ -1,24 +0,0 @@
-
-(* 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 memo :
- ('a * 'b) list ref -> ('a -> 'a -> bool) -> 'b -> ('a -> 'b) -> 'a -> 'b
-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/groebner/vo.itarget b/plugins/groebner/vo.itarget
deleted file mode 100644
index ad9d0e4ef..000000000
--- a/plugins/groebner/vo.itarget
+++ /dev/null
@@ -1,2 +0,0 @@
-GroebnerR.vo
-GroebnerZ.vo
diff --git a/plugins/pluginsbyte.itarget b/plugins/pluginsbyte.itarget
index c76c19942..04cbdccb0 100644
--- a/plugins/pluginsbyte.itarget
+++ b/plugins/pluginsbyte.itarget
@@ -13,7 +13,7 @@ xml/xml_plugin.cma
subtac/subtac_plugin.cma
ring/ring_plugin.cma
cc/cc_plugin.cma
-groebner/groebner_plugin.cma
+nsatz/nsatz_plugin.cma
funind/recdef_plugin.cma
syntax/ascii_syntax_plugin.cma
syntax/nat_syntax_plugin.cma
diff --git a/plugins/pluginsdyn.itarget b/plugins/pluginsdyn.itarget
index 559adfb6e..bbadfe696 100644
--- a/plugins/pluginsdyn.itarget
+++ b/plugins/pluginsdyn.itarget
@@ -13,7 +13,7 @@ xml/xml_plugin.cmxs
subtac/subtac_plugin.cmxs
ring/ring_plugin.cmxs
cc/cc_plugin.cmxs
-groebner/groebner_plugin.cmxs
+nsatz/nsatz_plugin.cmxs
funind/recdef_plugin.cmxs
syntax/ascii_syntax_plugin.cmxs
syntax/nat_syntax_plugin.cmxs
diff --git a/plugins/pluginsopt.itarget b/plugins/pluginsopt.itarget
index 62868a2c4..74b3f5275 100644
--- a/plugins/pluginsopt.itarget
+++ b/plugins/pluginsopt.itarget
@@ -13,7 +13,7 @@ xml/xml_plugin.cmxa
subtac/subtac_plugin.cmxa
ring/ring_plugin.cmxa
cc/cc_plugin.cmxa
-groebner/groebner_plugin.cmxa
+nsatz/nsatz_plugin.cmxa
funind/recdef_plugin.cmxa
syntax/ascii_syntax_plugin.cmxa
syntax/nat_syntax_plugin.cmxa
diff --git a/plugins/pluginsvo.itarget b/plugins/pluginsvo.itarget
index 48cbae3ec..033d1a199 100644
--- a/plugins/pluginsvo.itarget
+++ b/plugins/pluginsvo.itarget
@@ -2,7 +2,7 @@ dp/vo.otarget
field/vo.otarget
fourier/vo.otarget
funind/vo.otarget
-groebner/vo.otarget
+nsatz/vo.otarget
micromega/vo.otarget
omega/vo.otarget
quote/vo.otarget
@@ -10,4 +10,3 @@ ring/vo.otarget
romega/vo.otarget
rtauto/vo.otarget
setoid_ring/vo.otarget
-extraction/vo.itarget \ No newline at end of file
diff --git a/test-suite/success/Nsatz.v b/test-suite/success/Nsatz.v
new file mode 100644
index 000000000..fde9f4700
--- /dev/null
+++ b/test-suite/success/Nsatz.v
@@ -0,0 +1,216 @@
+Require Import NsatzR ZArith Reals List Ring_polynom.
+
+Section Examples.
+
+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.
+
+Lemma example1 : forall x y,
+ x+y=0 ->
+ x*y=0 ->
+ x^2=0.
+Proof.
+ nsatzR.
+Qed.
+
+Lemma example2 : forall x, x^2=0 -> x=0.
+Proof.
+ nsatzR.
+Qed.
+
+(*
+Notation X := (PEX Z 3).
+Notation Y := (PEX Z 2).
+Notation Z_ := (PEX Z 1).
+*)
+Lemma example3 : forall x y z,
+ x+y+z=0 ->
+ x*y+x*z+y*z=0->
+ x*y*z=0 -> x^3=0.
+Proof.
+Time nsatzR.
+Qed.
+
+(*
+Notation X := (PEX Z 4).
+Notation Y := (PEX Z 3).
+Notation Z_ := (PEX Z 2).
+Notation U := (PEX Z 1).
+*)
+Lemma example4 : 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.
+Proof.
+Time nsatzR.
+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").
+*)
+
+Lemma example5 : 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->
+ x*y*z+x*y*u+x*y*v+x*z*u+x*z*v+x*u*v+y*z*u+y*z*v+y*u*v+z*u*v=0->
+ x*y*z*u+y*z*u*v+z*u*v*x+u*v*x*y+v*x*y*z=0 ->
+ x*y*z*u*v=0 -> x^5=0.
+Proof.
+Time nsatzR.
+Qed.
+
+End Examples.
+
+Section Geometry.
+Require Export Reals NsatzR.
+Open Scope R_scope.
+
+Record point:Type:={
+ X:R;
+ Y:R}.
+
+Definition collinear(A B C:point):=
+ (X A - X B)*(Y C - Y B)-(Y A - Y B)*(X C - X B)=0.
+
+Definition parallel (A B C D:point):=
+ ((X A)-(X B))*((Y C)-(Y D))=((Y A)-(Y B))*((X C)-(X D)).
+
+Definition notparallel (A B C D:point)(x:R):=
+ x*(((X A)-(X B))*((Y C)-(Y D))-((Y A)-(Y B))*((X C)-(X D)))=1.
+
+Definition orthogonal (A B C D:point):=
+ ((X A)-(X B))*((X C)-(X D))+((Y A)-(Y B))*((Y C)-(Y D))=0.
+
+Definition equal2(A B:point):=
+ (X A)=(X B) /\ (Y A)=(Y B).
+
+Definition equal3(A B:point):=
+ ((X A)-(X B))^2+((Y A)-(Y B))^2 = 0.
+
+Definition nequal2(A B:point):=
+ (X A)<>(X B) \/ (Y A)<>(Y B).
+
+Definition nequal3(A B:point):=
+ not (((X A)-(X B))^2+((Y A)-(Y B))^2 = 0).
+
+Definition middle(A B I:point):=
+ 2*(X I)=(X A)+(X B) /\ 2*(Y I)=(Y A)+(Y B).
+
+Definition distance2(A B:point):=
+ (X B - X A)^2 + (Y B - Y A)^2.
+
+(* AB = CD *)
+Definition samedistance2(A B C D:point):=
+ (X B - X A)^2 + (Y B - Y A)^2 = (X D - X C)^2 + (Y D - Y C)^2.
+Definition determinant(A O B:point):=
+ (X A - X O)*(Y B - Y O) - (Y A - Y O)*(X B - X O).
+Definition scalarproduct(A O B:point):=
+ (X A - X O)*(X B - X O) + (Y A - Y O)*(Y B - Y O).
+Definition norm2(A O B:point):=
+ ((X A - X O)^2+(Y A - Y O)^2)*((X B - X O)^2+(Y B - Y O)^2).
+
+
+Lemma a1:forall A B C:Prop, ((A\/B)/\(A\/C)) -> (A\/(B/\C)).
+intuition.
+Qed.
+
+Lemma a2:forall A B C:Prop, ((A\/C)/\(B\/C)) -> ((A/\B)\/C).
+intuition.
+Qed.
+
+Lemma a3:forall a b c d:R, (a-b)*(c-d)=0 -> (a=b \/ c=d).
+intros.
+assert ( (a-b = 0) \/ (c-d = 0)).
+apply Rmult_integral.
+trivial.
+destruct H0.
+left; nsatz.
+right; nsatz.
+Qed.
+
+Ltac geo_unfold :=
+ unfold collinear; unfold parallel; unfold notparallel; unfold orthogonal;
+ unfold equal2; unfold equal3; unfold nequal2; unfold nequal3;
+ unfold middle; unfold samedistance2;
+ unfold determinant; unfold scalarproduct; unfold norm2; unfold distance2.
+
+Ltac geo_end :=
+ repeat (
+ repeat (match goal with h:_/\_ |- _ => decompose [and] h; clear h end);
+ repeat (apply a1 || apply a2 || apply a3);
+ repeat split).
+
+Ltac geo_rewrite_hyps:=
+ repeat (match goal with
+ | h:X _ = _ |- _ => rewrite h in *; clear h
+ | h:Y _ = _ |- _ => rewrite h in *; clear h
+ end).
+
+Ltac geo_begin:=
+ geo_unfold;
+ intros;
+ geo_rewrite_hyps;
+ geo_end.
+
+(* Examples *)
+
+Lemma Thales: forall O A B C D:point,
+ collinear O A C -> collinear O B D ->
+ parallel A B C D ->
+ (distance2 O B * distance2 O C = distance2 O D * distance2 O A
+ /\ distance2 O B * distance2 C D = distance2 O D * distance2 A B)
+ \/ collinear O A B.
+repeat geo_begin.
+(*
+Time nsatz.
+*)
+Time nsatz without sugar.
+(*
+Time nsatz with lexico sugar.
+Time nsatz with lexico.
+*)
+(*
+Time nsatzRpv 1%N 1%Z (@nil R) (@nil R). (* revlex, sugar, no div *)
+(*Finished transaction in 1. secs (0.479927u,0.s)*)
+Time nsatzRpv 1%N 0%Z (@nil R) (@nil R). (* revlex, no sugar, no div *)
+(*Finished transaction in 0. secs (0.543917u,0.s)*)
+Time nsatzRpv 1%N 2%Z (@nil R) (@nil R). (* lex, no sugar, no div *)
+(*Finished transaction in 0. secs (0.586911u,0.s)*)
+Time nsatzRpv 1%N 3%Z (@nil R) (@nil R). (* lex, sugar, no div *)
+(*Finished transaction in 0. secs (0.481927u,0.s)*)
+Time nsatzRpv 1%N 5%Z (@nil R) (@nil R). (* revlex, sugar, div *)
+(*Finished transaction in 1. secs (0.601909u,0.s)*)
+*)
+Time nsatz.
+Qed.
+
+Lemma hauteurs:forall A B C A1 B1 C1 H:point,
+ collinear B C A1 -> orthogonal A A1 B C ->
+ collinear A C B1 -> orthogonal B B1 A C ->
+ collinear A B C1 -> orthogonal C C1 A B ->
+ collinear A A1 H -> collinear B B1 H ->
+
+ collinear C C1 H
+ \/ collinear A B C.
+
+geo_begin.
+Time nsatz.
+(*Finished transaction in 3. secs (2.43263u,0.010998s)*)
+Qed.
+
+End Geometry.
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-21 18:44:51 +0000