diff options
| author | barras | 2009-03-05 21:20:08 +0000 |
|---|---|---|
| committer | barras | 2009-03-05 21:20:08 +0000 |
| commit | 55c214eb0651974dc5a9f47f9bdae9fbd6960d57 (patch) | |
| tree | 5677d763c8663880d528089b8b7b1e2054e7ca4b | |
| parent | f8ae76ca7965e176094496bf17dd77bff1acae10 (diff) | |
ajout de la tactique groebner de Loic Pottier
git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/coq/trunk@11964 85f007b7-540e-0410-9357-904b9bb8a0f7
| -rw-r--r-- | Makefile.build | 26 | ||||
| -rw-r--r-- | Makefile.common | 423 | ||||
| -rw-r--r-- | contrib/groebner/GroebnerR.v | 465 | ||||
| -rw-r--r-- | contrib/groebner/GroebnerZ.v | 67 | ||||
| -rw-r--r-- | contrib/groebner/groebner.ml4 | 450 | ||||
| -rw-r--r-- | contrib/groebner/ideal.ml | 792 | ||||
| -rw-r--r-- | contrib/groebner/ideal.mli | 82 | ||||
| -rw-r--r-- | contrib/groebner/polynom.ml | 1055 | ||||
| -rw-r--r-- | contrib/groebner/polynom.mli | 111 | ||||
| -rw-r--r-- | contrib/groebner/utile.ml | 228 | ||||
| -rw-r--r-- | contrib/groebner/utile.mli | 34 |
11 files changed, 3485 insertions, 248 deletions
diff --git a/Makefile.build b/Makefile.build index 14e0cd211a..5513c854a0 100644 --- a/Makefile.build +++ b/Makefile.build @@ -62,18 +62,7 @@ else HIDE = @ endif -LOCALINCLUDES=-I config -I tools -I tools/coqdoc \ - -I scripts -I lib -I kernel -I kernel/byterun -I library \ - -I proofs -I tactics -I pretyping \ - -I interp -I toplevel -I parsing -I ide/utils -I ide \ - -I contrib/omega -I contrib/romega -I contrib/micromega \ - -I contrib/quote \ - -I contrib/ring -I contrib/dp -I contrib/setoid_ring \ - -I contrib/xml -I contrib/extraction \ - -I contrib/interface -I contrib/fourier -I contrib/cc \ - -I contrib/funind -I contrib/firstorder \ - -I contrib/field -I contrib/subtac -I contrib/rtauto - +LOCALINCLUDES=$(addprefix -I , $(SRCDIRS) ) MLINCLUDES=$(LOCALINCLUDES) -I $(MYCAMLP4LIB) OCAMLC += $(CAMLFLAGS) @@ -87,12 +76,12 @@ CAMLP4EXTENDFLAGS=-I . #grammar dependencies are now in camlp4use statements CAMLP4DEPS=sed -n -e 's@^(\*.*camlp4deps: "\(.*\)".*\*)@\1@p' CAMLP4USE=sed -n -e 's/pa_macro.cmo/pa_macro.cmo -D$(CAMLVERSION)/' -e 's@^(\*.*camlp4use: "\(.*\)".*\*)@\1@p' -COQINCLUDES= # coqtop includes itself the needed paths +COQINCLUDES= # coqtop includes itself the needed paths COQ_XML= # is "-xml" when building XML library -VM= # is "-no-vm" to not use the vm" -UNBOXEDVALUES= # is "-unboxed-values" to use unboxed values +VM= # is "-no-vm" to not use the vm" +UNBOXEDVALUES= # is "-unboxed-values" to use unboxed values COQOPTS=$(COQ_XML) $(VM) $(UNBOXEDVALUES) -TIME= # is "'time -p'" to get compilation time of .v +TIME= # is "'time -p'" to get compilation time of .v BOOTCOQTOP:=$(TIME) $(BESTCOQTOP) -boot $(COQOPTS) @@ -190,7 +179,7 @@ $(COQTOPBYTE): $(COQMKTOP) $(LINKCMO) $(LIBCOQRUN) $(COQTOPEXE): $(ORDER_ONLY_SEP) $(BESTCOQTOP) cd bin; ln -sf coqtop.$(BEST)$(EXE) coqtop$(EXE) -LOCALCHKLIBS:=-I checker -I lib -I config -I kernel +LOCALCHKLIBS:=$(addprefix -I , $(CHKSRCDIRS) ) CHKLIBS:=$(LOCALCHKLIBS) -I $(MYCAMLP4LIB) CHKBYTEFLAGS:=$(CHKLIBS) $(CAMLDEBUG) $(USERFLAGS) CHKOPTFLAGS:=$(CHKLIBS) $(CAMLDEBUGOPT) $(CAMLTIMEPROF) $(USERFLAGS) @@ -345,6 +334,9 @@ $(FUNINDCMA:.cma=.cmxa): $(FUNINDCMO:.cmo=.cmx) $(SUBTACCMA): $(SUBTACCMO) $(SUBTACCMA:.cma=.cmxa): $(SUBTACCMO:.cmo=.cmx) +$(GBCMA): $(GBCMO) +$(GBCMA:.cma=.cmxa): $(GBCMO:.cmo=.cmx) + %.cma: $(SHOW)'OCAMLC -a -o $@' $(HIDE)$(OCAMLC) $(BYTEFLAGS) -a -o $@ $^ diff --git a/Makefile.common b/Makefile.common index 636bf323f2..db78a65ed3 100644 --- a/Makefile.common +++ b/Makefile.common @@ -55,6 +55,21 @@ OTHERBINARIES:=$(COQMKTOPBYTE) $(COQCBYTE) CSDPCERT:=contrib/micromega/csdpcert$(EXE) +SRCDIRS:=\ + config tools tools/coqdoc scripts lib \ + kernel kernel/byterun library proofs tactics \ + pretyping interp toplevel parsing ide/utils \ + ide \ + $(addprefix contrib/, \ + omega romega micromega quote ring dp \ + setoid_ring xml extraction interface fourier \ + cc funind firstorder field subtac \ + rtauto groebner ) + +# Order is relevent here because kernel and checker contain files +# with the same name +CHKSRCDIRS:= checker lib config kernel + ########################################################################### # tools ########################################################################### @@ -85,25 +100,22 @@ COQTEXOPTS:=-n 72 -image "$(COQSRC)/$(COQTOPEXE) -boot" -sl -small DOCCOMMON:=doc/common/version.tex doc/common/title.tex doc/common/macros.tex -REFMANCOQTEXFILES:=\ - doc/refman/RefMan-gal.v.tex doc/refman/RefMan-ext.v.tex \ - doc/refman/RefMan-mod.v.tex doc/refman/RefMan-tac.v.tex \ - doc/refman/RefMan-cic.v.tex doc/refman/RefMan-lib.v.tex \ - doc/refman/RefMan-tacex.v.tex doc/refman/RefMan-syn.v.tex \ - doc/refman/RefMan-oth.v.tex doc/refman/RefMan-ltac.v.tex \ - doc/refman/RefMan-decl.v.tex \ - doc/refman/Cases.v.tex doc/refman/Coercion.v.tex doc/refman/Extraction.v.tex \ - doc/refman/Program.v.tex doc/refman/Omega.v.tex doc/refman/Micromega.v.tex doc/refman/Polynom.v.tex \ - doc/refman/Setoid.v.tex doc/refman/Helm.tex doc/refman/Classes.v.tex - -REFMANTEXFILES:=\ - doc/refman/headers.sty \ - doc/refman/Reference-Manual.tex doc/refman/RefMan-pre.tex \ - doc/refman/RefMan-int.tex doc/refman/RefMan-pro.tex \ - doc/refman/RefMan-com.tex \ - doc/refman/RefMan-uti.tex doc/refman/RefMan-ide.tex \ - doc/refman/RefMan-add.tex doc/refman/RefMan-modr.tex \ - doc/refman/ExternalProvers.tex \ +REFMANCOQTEXFILES:=$(addprefix doc/refman/, \ + RefMan-gal.v.tex RefMan-ext.v.tex \ + RefMan-mod.v.tex RefMan-tac.v.tex \ + RefMan-cic.v.tex RefMan-lib.v.tex \ + RefMan-tacex.v.tex RefMan-syn.v.tex \ + 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 \ + Setoid.v.tex Helm.tex Classes.v.tex ) + +REFMANTEXFILES:=$(addprefix doc/refman/, \ + headers.sty Reference-Manual.tex \ + RefMan-pre.tex RefMan-int.tex RefMan-pro.tex RefMan-com.tex \ + RefMan-uti.tex RefMan-ide.tex RefMan-add.tex RefMan-modr.tex \ + ExternalProvers.tex ) \ $(REFMANCOQTEXFILES) \ REFMANEPSFILES:=doc/refman/coqide.eps doc/refman/coqide-queries.eps @@ -129,56 +141,42 @@ CAMLP4OBJS:=gramlib.cma CONFIG:=\ config/coq_config.cmo -LIBREP:=\ - lib/pp_control.cmo lib/pp.cmo lib/compat.cmo lib/flags.cmo \ - lib/util.cmo lib/bigint.cmo lib/hashcons.cmo lib/dyn.cmo lib/system.cmo \ - lib/envars.cmo lib/bstack.cmo lib/edit.cmo lib/gset.cmo lib/gmap.cmo \ - lib/tlm.cmo lib/gmapl.cmo lib/profile.cmo lib/explore.cmo \ - lib/predicate.cmo lib/rtree.cmo lib/heap.cmo lib/option.cmo lib/dnet.cmo +LIBREP:=$(addprefix lib/, \ + pp_control.cmo pp.cmo compat.cmo flags.cmo \ + util.cmo bigint.cmo hashcons.cmo dyn.cmo \ + system.cmo envars.cmo bstack.cmo edit.cmo \ + gset.cmo gmap.cmo tlm.cmo gmapl.cmo \ + profile.cmo explore.cmo predicate.cmo rtree.cmo \ + heap.cmo option.cmo dnet.cmo ) # Rem: Cygwin already uses variable LIB -BYTERUN:=\ - kernel/byterun/coq_fix_code.o kernel/byterun/coq_memory.o \ - kernel/byterun/coq_values.o kernel/byterun/coq_interp.o - -KERNEL:=\ - kernel/names.cmo kernel/univ.cmo \ - kernel/esubst.cmo kernel/term.cmo \ - kernel/mod_subst.cmo kernel/sign.cmo \ - kernel/cbytecodes.cmo kernel/copcodes.cmo \ - kernel/cemitcodes.cmo kernel/vm.cmo \ - kernel/declarations.cmo \ - kernel/retroknowledge.cmo kernel/pre_env.cmo \ - kernel/cbytegen.cmo kernel/environ.cmo \ - kernel/csymtable.cmo kernel/conv_oracle.cmo \ - kernel/closure.cmo kernel/reduction.cmo kernel/type_errors.cmo \ - kernel/entries.cmo kernel/modops.cmo \ - kernel/inductive.cmo kernel/vconv.cmo kernel/typeops.cmo \ - kernel/indtypes.cmo kernel/cooking.cmo kernel/term_typing.cmo \ - kernel/subtyping.cmo kernel/mod_typing.cmo kernel/safe_typing.cmo - -LIBRARY:=\ - library/nameops.cmo library/libnames.cmo library/libobject.cmo \ - library/summary.cmo library/nametab.cmo library/global.cmo library/lib.cmo \ - library/declaremods.cmo library/library.cmo \ - library/states.cmo \ - library/decl_kinds.cmo library/dischargedhypsmap.cmo library/goptions.cmo \ - library/decls.cmo library/heads.cmo - -PRETYPING:=\ - pretyping/termops.cmo pretyping/evd.cmo \ - pretyping/reductionops.cmo pretyping/vnorm.cmo pretyping/inductiveops.cmo \ - pretyping/retyping.cmo pretyping/cbv.cmo \ - pretyping/pretype_errors.cmo \ - pretyping/typing.cmo pretyping/evarutil.cmo pretyping/term_dnet.cmo \ - pretyping/recordops.cmo \ - pretyping/tacred.cmo pretyping/evarconv.cmo \ - pretyping/typeclasses_errors.cmo pretyping/typeclasses.cmo \ - pretyping/classops.cmo pretyping/coercion.cmo \ - pretyping/unification.cmo pretyping/clenv.cmo \ - pretyping/rawterm.cmo pretyping/pattern.cmo \ - pretyping/detyping.cmo pretyping/indrec.cmo\ - pretyping/cases.cmo pretyping/pretyping.cmo pretyping/matching.cmo +BYTERUN:=$(addprefix kernel/byterun/, \ + coq_fix_code.o coq_memory.o coq_values.o coq_interp.o ) + +KERNEL:=$(addprefix kernel/, \ + names.cmo univ.cmo esubst.cmo term.cmo \ + mod_subst.cmo sign.cmo cbytecodes.cmo copcodes.cmo \ + cemitcodes.cmo vm.cmo declarations.cmo \ + retroknowledge.cmo pre_env.cmo cbytegen.cmo environ.cmo \ + csymtable.cmo conv_oracle.cmo closure.cmo reduction.cmo \ + type_errors.cmo entries.cmo modops.cmo inductive.cmo \ + vconv.cmo typeops.cmo indtypes.cmo cooking.cmo \ + term_typing.cmo subtyping.cmo mod_typing.cmo safe_typing.cmo ) + +LIBRARY:=$(addprefix library/, \ + nameops.cmo libnames.cmo libobject.cmo summary.cmo \ + nametab.cmo global.cmo lib.cmo declaremods.cmo \ + library.cmo states.cmo decl_kinds.cmo dischargedhypsmap.cmo \ + goptions.cmo decls.cmo heads.cmo ) + +PRETYPING:=$(addprefix pretyping/, \ + termops.cmo evd.cmo reductionops.cmo vnorm.cmo \ + inductiveops.cmo retyping.cmo cbv.cmo pretype_errors.cmo \ + typing.cmo evarutil.cmo term_dnet.cmo recordops.cmo \ + tacred.cmo evarconv.cmo typeclasses_errors.cmo \ + typeclasses.cmo classops.cmo coercion.cmo unification.cmo \ + clenv.cmo rawterm.cmo pattern.cmo detyping.cmo \ + indrec.cmo cases.cmo pretyping.cmo matching.cmo ) INTERP:=\ parsing/lexer.cmo interp/topconstr.cmo interp/ppextend.cmo \ @@ -188,26 +186,21 @@ INTERP:=\ interp/modintern.cmo interp/constrextern.cmo interp/coqlib.cmo \ toplevel/discharge.cmo library/declare.cmo -PROOFS:=\ - proofs/tacexpr.cmo proofs/proof_type.cmo proofs/redexpr.cmo \ - proofs/proof_trees.cmo proofs/logic.cmo \ - proofs/refiner.cmo proofs/evar_refiner.cmo proofs/tacmach.cmo \ - proofs/pfedit.cmo proofs/tactic_debug.cmo \ - proofs/clenvtac.cmo proofs/decl_mode.cmo - -PARSING:=\ - parsing/extend.cmo \ - parsing/pcoq.cmo parsing/egrammar.cmo parsing/g_xml.cmo \ - parsing/ppconstr.cmo parsing/printer.cmo \ - parsing/pptactic.cmo parsing/ppdecl_proof.cmo parsing/tactic_printer.cmo \ - parsing/printmod.cmo parsing/prettyp.cmo - -HIGHPARSING:=\ - parsing/g_constr.cmo parsing/g_vernac.cmo parsing/g_prim.cmo \ - parsing/g_proofs.cmo parsing/g_tactic.cmo parsing/g_ltac.cmo \ - parsing/g_natsyntax.cmo parsing/g_zsyntax.cmo parsing/g_rsyntax.cmo \ - parsing/g_ascii_syntax.cmo parsing/g_string_syntax.cmo \ - parsing/g_decl_mode.cmo parsing/g_intsyntax.cmo +PROOFS:=$(addprefix proofs/, \ + tacexpr.cmo proof_type.cmo redexpr.cmo proof_trees.cmo \ + logic.cmo refiner.cmo evar_refiner.cmo tacmach.cmo \ + pfedit.cmo tactic_debug.cmo clenvtac.cmo decl_mode.cmo ) + +PARSING:=$(addprefix parsing/, \ + extend.cmo pcoq.cmo egrammar.cmo g_xml.cmo \ + ppconstr.cmo printer.cmo pptactic.cmo ppdecl_proof.cmo \ + tactic_printer.cmo printmod.cmo prettyp.cmo ) + +HIGHPARSING:=$(addprefix parsing/, \ + g_constr.cmo g_vernac.cmo g_prim.cmo g_proofs.cmo \ + g_tactic.cmo g_ltac.cmo g_natsyntax.cmo g_zsyntax.cmo \ + g_rsyntax.cmo g_ascii_syntax.cmo g_string_syntax.cmo \ + g_decl_mode.cmo g_intsyntax.cmo ) TACTICS:=\ tactics/dn.cmo tactics/termdn.cmo tactics/btermdn.cmo \ @@ -233,27 +226,23 @@ TOPLEVEL:=\ toplevel/toplevel.cmo $(REVISIONCMO) toplevel/usage.cmo \ toplevel/coqinit.cmo toplevel/coqtop.cmo -HIGHTACTICS:=\ - tactics/refine.cmo tactics/extraargs.cmo \ - tactics/extratactics.cmo tactics/eauto.cmo tactics/class_tactics.cmo \ - tactics/tauto.cmo tactics/eqdecide.cmo +HIGHTACTICS:=$(addprefix tactics/, \ + refine.cmo extraargs.cmo extratactics.cmo \ + eauto.cmo class_tactics.cmo tauto.cmo \ + eqdecide.cmo ) OMEGACMA:=contrib/omega/omega_plugin.cma -OMEGACMO:=\ - contrib/omega/omega.cmo contrib/omega/coq_omega.cmo \ - contrib/omega/g_omega.cmo +OMEGACMO:=$(addprefix contrib/omega/, \ + omega.cmo coq_omega.cmo g_omega.cmo ) ROMEGACMA:=contrib/romega/romega_plugin.cma -ROMEGACMO:=\ - contrib/romega/const_omega.cmo \ - contrib/romega/refl_omega.cmo contrib/romega/g_romega.cmo +ROMEGACMO:=$(addprefix contrib/romega/, \ + const_omega.cmo refl_omega.cmo g_romega.cmo ) MICROMEGACMA:=contrib/micromega/micromega_plugin.cma -MICROMEGACMO:=\ - contrib/micromega/mutils.cmo \ - contrib/micromega/micromega.cmo contrib/micromega/mfourier.cmo \ - contrib/micromega/certificate.cmo \ - contrib/micromega/coq_micromega.cmo contrib/micromega/g_micromega.cmo +MICROMEGACMO:=$(addprefix contrib/micromega/, \ + mutils.cmo micromega.cmo mfourier.cmo certificate.cmo \ + coq_micromega.cmo g_micromega.cmo ) QUOTECMA:=contrib/quote/quote_plugin.cma QUOTECMO:=\ @@ -266,82 +255,69 @@ RINGCMO:=\ NEWRINGCMA:=contrib/setoid_ring/newring_plugin.cma NEWRINGCMO:=contrib/setoid_ring/newring.cmo +GBCMA:=contrib/groebner/groebner_plugin.cma +GBCMO:=$(addprefix contrib/groebner/, \ + utile.cmo polynom.cmo ideal.cmo groebner.cmo ) + DPCMA:=contrib/dp/dp_plugin.cma -DPCMO:=contrib/dp/dp_why.cmo contrib/dp/dp_zenon.cmo \ - contrib/dp/dp.cmo contrib/dp/dp_gappa.cmo contrib/dp/g_dp.cmo +DPCMO:=$(addprefix contrib/dp/, \ + dp_why.cmo dp_zenon.cmo dp.cmo dp_gappa.cmo g_dp.cmo ) FIELDCMO:=contrib/field/field.cmo FIELDCMA:=contrib/field/field_plugin.cma XMLCMA:=contrib/xml/xml_plugin.cma -XMLCMO:=\ - contrib/xml/unshare.cmo contrib/xml/xml.cmo contrib/xml/acic.cmo \ - contrib/xml/doubleTypeInference.cmo \ - contrib/xml/cic2acic.cmo contrib/xml/acic2Xml.cmo \ - contrib/xml/proof2aproof.cmo \ - contrib/xml/xmlcommand.cmo contrib/xml/proofTree2Xml.cmo \ - contrib/xml/xmlentries.cmo contrib/xml/cic2Xml.cmo \ - contrib/xml/dumptree.cmo +XMLCMO:=$(addprefix contrib/xml/, \ + unshare.cmo xml.cmo acic.cmo doubleTypeInference.cmo \ + cic2acic.cmo acic2Xml.cmo proof2aproof.cmo \ + xmlcommand.cmo proofTree2Xml.cmo xmlentries.cmo cic2Xml.cmo \ + dumptree.cmo ) FOURIERCMA:=contrib/fourier/fourier_plugin.cma -FOURIERCMO:=\ - contrib/fourier/fourier.cmo contrib/fourier/fourierR.cmo \ - contrib/fourier/g_fourier.cmo +FOURIERCMO:=$(addprefix contrib/fourier/, \ + fourier.cmo fourierR.cmo g_fourier.cmo ) EXTRACTIONCMA:=contrib/extraction/extraction_plugin.cma -EXTRACTIONCMO:=\ - contrib/extraction/table.cmo\ - contrib/extraction/mlutil.cmo\ - contrib/extraction/modutil.cmo \ - contrib/extraction/extraction.cmo \ - contrib/extraction/common.cmo \ - contrib/extraction/ocaml.cmo \ - contrib/extraction/haskell.cmo \ - contrib/extraction/scheme.cmo \ - contrib/extraction/extract_env.cmo \ - contrib/extraction/g_extraction.cmo +EXTRACTIONCMO:=$(addprefix contrib/extraction/, \ + table.cmo mlutil.cmo modutil.cmo extraction.cmo common.cmo \ + ocaml.cmo haskell.cmo scheme.cmo extract_env.cmo g_extraction.cmo ) FUNINDCMA:=contrib/funind/recdef_plugin.cma -FUNINDCMO:=\ - contrib/funind/indfun_common.cmo contrib/funind/rawtermops.cmo \ - contrib/funind/recdef.cmo \ - contrib/funind/rawterm_to_relation.cmo \ - contrib/funind/functional_principles_proofs.cmo \ - contrib/funind/functional_principles_types.cmo \ - contrib/funind/invfun.cmo contrib/funind/indfun.cmo \ - contrib/funind/merge.cmo contrib/funind/g_indfun.cmo +FUNINDCMO:=$(addprefix contrib/funind/, \ + indfun_common.cmo rawtermops.cmo \ + recdef.cmo rawterm_to_relation.cmo \ + functional_principles_proofs.cmo \ + functional_principles_types.cmo \ + invfun.cmo indfun.cmo merge.cmo g_indfun.cmo ) FOCMA:=contrib/firstorder/ground_plugin.cma -FOCMO:=\ - contrib/firstorder/formula.cmo contrib/firstorder/unify.cmo \ - contrib/firstorder/sequent.cmo contrib/firstorder/rules.cmo \ - contrib/firstorder/instances.cmo contrib/firstorder/ground.cmo \ - contrib/firstorder/g_ground.cmo +FOCMO:=$(addprefix contrib/firstorder/, \ + formula.cmo unify.cmo sequent.cmo rules.cmo instances.cmo ground.cmo \ + g_ground.cmo ) CCCMA:=contrib/cc/cc_plugin.cma -CCCMO:=contrib/cc/ccalgo.cmo contrib/cc/ccproof.cmo contrib/cc/cctac.cmo \ - contrib/cc/g_congruence.cmo +CCCMO:=$(addprefix contrib/cc/, \ + ccalgo.cmo ccproof.cmo cctac.cmo g_congruence.cmo ) SUBTACCMA:=contrib/subtac/subtac_plugin.cma -SUBTACCMO:=contrib/subtac/subtac_utils.cmo contrib/subtac/eterm.cmo \ - contrib/subtac/g_eterm.cmo \ - contrib/subtac/subtac_errors.cmo contrib/subtac/subtac_coercion.cmo \ - contrib/subtac/subtac_obligations.cmo contrib/subtac/subtac_cases.cmo \ - contrib/subtac/subtac_pretyping_F.cmo contrib/subtac/subtac_pretyping.cmo \ - contrib/subtac/subtac_command.cmo contrib/subtac/subtac_classes.cmo \ - contrib/subtac/subtac.cmo contrib/subtac/g_subtac.cmo \ - contrib/subtac/equations.cmo +SUBTACCMO:=$(addprefix contrib/subtac/, \ + subtac_utils.cmo eterm.cmo g_eterm.cmo \ + subtac_errors.cmo subtac_coercion.cmo \ + subtac_obligations.cmo subtac_cases.cmo \ + subtac_pretyping_F.cmo subtac_pretyping.cmo \ + subtac_command.cmo subtac_classes.cmo \ + subtac.cmo g_subtac.cmo equations.cmo ) RTAUTOCMA:=contrib/rtauto/rtauto_plugin.cma -RTAUTOCMO:=contrib/rtauto/proof_search.cmo contrib/rtauto/refl_tauto.cmo \ - contrib/rtauto/g_rtauto.cmo +RTAUTOCMO:=$(addprefix contrib/rtauto/, \ + proof_search.cmo refl_tauto.cmo g_rtauto.cmo ) CONTRIBSTATIC:=\ $(OMEGACMO) $(ROMEGACMO) $(MICROMEGACMO) \ $(QUOTECMO) $(RINGCMO) $(NEWRINGCMO) $(DPCMO) $(FIELDCMO) \ $(FOURIERCMO) $(EXTRACTIONCMO) $(XMLCMO) \ $(CCCMO) $(FOCMO) $(SUBTACCMO) $(RTAUTOCMO) \ - $(FUNINDCMO) + $(FUNINDCMO) $(GBCMO) CONTRIBCMA:=contrib/contrib.cma @@ -353,7 +329,7 @@ ifneq ($(HASNATDYNLINK),false) PLUGINS:=$(INITPLUGINS) \ $(OMEGACMA) $(ROMEGACMA) $(MICROMEGACMA) \ $(QUOTECMA) $(RINGCMA) $(NEWRINGCMA) $(FIELDCMA) \ - $(FOURIERCMA) $(RTAUTOCMA) + $(FOURIERCMA) $(RTAUTOCMA) $(GBCMA) INITPLUGINSOPT:=$(INITPLUGINS:.cma=.cmxs) PLUGINSOPT:=$(PLUGINS:.cma=.cmxs) ifeq ($(BEST),opt) @@ -383,21 +359,23 @@ OBJSCMO:=$(CONFIG) $(LIBREP) $(KERNEL) $(LIBRARY) $(PRETYPING) $(INTERP) \ $(PROOFS) $(PARSING) $(TACTICS) $(TOPLEVEL) $(HIGHPARSING) \ $(HIGHTACTICS) $(CONTRIBSTATIC) -COQIDECMO:=ide/utils/okey.cmo ide/utils/config_file.cmo \ - ide/utils/configwin_keys.cmo ide/utils/configwin_types.cmo \ - ide/utils/configwin_messages.cmo ide/utils/configwin_ihm.cmo \ - ide/utils/configwin.cmo \ - ide/utils/editable_cells.cmo ide/config_parser.cmo \ - ide/typed_notebook.cmo \ - ide/config_lexer.cmo ide/utf8_convert.cmo ide/preferences.cmo \ - ide/ideutils.cmo ide/blaster_window.cmo ide/undo.cmo \ - ide/highlight.cmo ide/coq.cmo ide/coq_commands.cmo \ - ide/coq_tactics.cmo ide/command_windows.cmo ide/coqide.cmo +COQIDECMO:=\ + $(addprefix ide/utils/, \ + okey.cmo config_file.cmo configwin_keys.cmo \ + configwin_types.cmo configwin_messages.cmo configwin_ihm.cmo \ + configwin.cmo editable_cells.cmo ) \ + $(addprefix ide/, \ + config_parser.cmo typed_notebook.cmo config_lexer.cmo \ + utf8_convert.cmo preferences.cmo ideutils.cmo \ + blaster_window.cmo undo.cmo highlight.cmo \ + coq.cmo coq_commands.cmo coq_tactics.cmo \ + command_windows.cmo coqide.cmo ) COQIDECMX:=$(COQIDECMO:.cmo=.cmx) -COQENVCMO:=$(CONFIG) lib/pp_control.cmo lib/pp.cmo lib/compat.cmo lib/flags.cmo \ - lib/util.cmo lib/system.cmo lib/envars.cmo +COQENVCMO:=$(CONFIG) \ + lib/pp_control.cmo lib/pp.cmo lib/compat.cmo lib/flags.cmo \ + lib/util.cmo lib/system.cmo lib/envars.cmo COQMKTOPCMO:=$(COQENVCMO) scripts/tolink.cmo scripts/coqmktop.cmo COQMKTOPCMX:=$(COQMKTOPCMO:.cmo=.cmx) @@ -405,16 +383,10 @@ COQMKTOPCMX:=$(COQMKTOPCMO:.cmo=.cmx) COQCCMO:=$(COQENVCMO) $(REVISIONCMO) toplevel/usage.cmo scripts/coqc.cmo COQCCMX:=$(COQCCMO:.cmo=.cmx) -INTERFACE:=\ - contrib/interface/vtp.cmo contrib/interface/xlate.cmo \ - contrib/interface/paths.cmo contrib/interface/translate.cmo \ - contrib/interface/pbp.cmo \ - contrib/interface/dad.cmo \ - contrib/interface/history.cmo \ - contrib/interface/name_to_ast.cmo contrib/interface/debug_tac.cmo \ - contrib/interface/showproof_ct.cmo contrib/interface/showproof.cmo \ - contrib/interface/blast.cmo contrib/interface/depends.cmo \ - contrib/interface/centaur.cmo +INTERFACE:=$(addprefix contrib/interface/, \ + vtp.cmo xlate.cmo paths.cmo translate.cmo pbp.cmo dad.cmo \ + history.cmo name_to_ast.cmo debug_tac.cmo showproof_ct.cmo showproof.cmo \ + blast.cmo depends.cmo centaur.cmo ) INTERFACECMX:=$(INTERFACE:.cmo=.cmx) @@ -439,9 +411,9 @@ else INTERFACERC:= contrib/interface/vernacrc contrib/interface/CoqParser.v endif -CSDPCERTCMO:= contrib/micromega/mutils.cmo contrib/micromega/micromega.cmo \ - contrib/micromega/mfourier.cmo contrib/micromega/certificate.cmo \ - contrib/micromega/sos.cmo contrib/micromega/csdpcert.cmo +CSDPCERTCMO:=$(addprefix contrib/micromega/, \ + mutils.cmo micromega.cmo mfourier.cmo certificate.cmo \ + sos.cmo csdpcert.cmo ) CSDPCERTCMX:= $(CSDPCERTCMO:.cmo=.cmx) DEBUGPRINTERS:=dev/top_printers.cmo dev/vm_printers.cmo dev/printers.cma @@ -452,15 +424,14 @@ COQDEPCMX:=$(COQDEPCMO:.cmo=.cmx) GALLINACMO:=tools/gallina_lexer.cmo tools/gallina.cmo GALLINACMX:=$(GALLINACMO:.cmo=.cmx) -COQDOCCMO:=$(CONFIG) tools/coqdoc/cdglobals.cmo tools/coqdoc/alpha.cmo \ - tools/coqdoc/index.cmo tools/coqdoc/output.cmo \ - tools/coqdoc/pretty.cmo tools/coqdoc/main.cmo +COQDOCCMO:=$(CONFIG) $(addprefix tools/coqdoc/, \ + cdglobals.cmo alpha.cmo index.cmo output.cmo pretty.cmo main.cmo ) COQDOCCMX:=$(COQDOCCMO:.cmo=.cmx) # checker MCHECKER:=\ - config/coq_config.cmo \ + $(CONFIG) \ lib/pp_control.cmo lib/pp.cmo lib/compat.cmo \ lib/flags.cmo lib/util.cmo lib/option.cmo lib/hashcons.cmo \ lib/system.cmo lib/envars.cmo \ @@ -480,27 +451,22 @@ MCHECKER:=\ # grammar modules with camlp4 GRAMMARNEEDEDCMO:=\ - lib/pp_control.cmo lib/pp.cmo lib/compat.cmo lib/flags.cmo \ - lib/util.cmo lib/bigint.cmo lib/dyn.cmo lib/hashcons.cmo lib/predicate.cmo \ - lib/rtree.cmo lib/option.cmo \ - kernel/names.cmo kernel/univ.cmo \ - kernel/esubst.cmo kernel/term.cmo kernel/mod_subst.cmo kernel/sign.cmo \ - kernel/cbytecodes.cmo kernel/copcodes.cmo kernel/cemitcodes.cmo \ - kernel/declarations.cmo \ - kernel/retroknowledge.cmo kernel/pre_env.cmo \ - kernel/cbytegen.cmo kernel/conv_oracle.cmo kernel/environ.cmo \ - kernel/closure.cmo kernel/reduction.cmo kernel/type_errors.cmo\ - kernel/entries.cmo \ - kernel/modops.cmo \ - kernel/inductive.cmo kernel/typeops.cmo \ - kernel/indtypes.cmo kernel/cooking.cmo kernel/term_typing.cmo \ - kernel/subtyping.cmo kernel/mod_typing.cmo kernel/safe_typing.cmo \ - library/nameops.cmo library/libnames.cmo library/summary.cmo \ - library/nametab.cmo library/libobject.cmo library/lib.cmo \ - library/goptions.cmo library/decl_kinds.cmo library/global.cmo \ - pretyping/termops.cmo pretyping/evd.cmo pretyping/reductionops.cmo \ - pretyping/inductiveops.cmo pretyping/rawterm.cmo pretyping/detyping.cmo \ - pretyping/pattern.cmo \ + $(addprefix lib/, \ + pp_control.cmo pp.cmo compat.cmo flags.cmo util.cmo \ + bigint.cmo dyn.cmo hashcons.cmo predicate.cmo rtree.cmo option.cmo ) \ + $(addprefix kernel/, \ + names.cmo univ.cmo esubst.cmo term.cmo mod_subst.cmo sign.cmo \ + cbytecodes.cmo copcodes.cmo cemitcodes.cmo declarations.cmo \ + retroknowledge.cmo pre_env.cmo cbytegen.cmo conv_oracle.cmo \ + environ.cmo closure.cmo reduction.cmo type_errors.cmo entries.cmo \ + modops.cmo inductive.cmo typeops.cmo indtypes.cmo cooking.cmo \ + term_typing.cmo subtyping.cmo mod_typing.cmo safe_typing.cmo ) \ + $(addprefix library/, \ + nameops.cmo libnames.cmo summary.cmo nametab.cmo libobject.cmo \ + lib.cmo goptions.cmo decl_kinds.cmo global.cmo ) \ + $(addprefix pretyping/, \ + termops.cmo evd.cmo reductionops.cmo inductiveops.cmo rawterm.cmo \ + detyping.cmo pattern.cmo ) \ interp/topconstr.cmo interp/genarg.cmo interp/ppextend.cmo \ proofs/tacexpr.cmo \ parsing/lexer.cmo parsing/extend.cmo \ @@ -510,9 +476,8 @@ GRAMMARNEEDEDCMO:=\ CAMLP4EXTENSIONSCMO:=\ parsing/argextend.cmo parsing/tacextend.cmo parsing/vernacextend.cmo -GRAMMARSCMO:=\ - parsing/g_prim.cmo parsing/g_tactic.cmo \ - parsing/g_ltac.cmo parsing/g_constr.cmo +GRAMMARSCMO:=$(addprefix parsing/, \ + g_prim.cmo g_tactic.cmo g_ltac.cmo g_constr.cmo ) GRAMMARCMO:=config/coq_config.cmo $(GRAMMARNEEDEDCMO) $(CAMLP4EXTENSIONSCMO) $(GRAMMARSCMO) @@ -530,30 +495,23 @@ STAGE1:=parsing/grammar.cma parsing/q_constr.cmo PRINTERSCMO:=\ config/coq_config.cmo lib/lib.cma \ - kernel/names.cmo kernel/univ.cmo kernel/esubst.cmo kernel/term.cmo \ - kernel/mod_subst.cmo kernel/copcodes.cmo kernel/cemitcodes.cmo \ - kernel/sign.cmo kernel/declarations.cmo kernel/retroknowledge.cmo \ - kernel/pre_env.cmo \ - kernel/retroknowledge.cmo kernel/pre_env.cmo \ - kernel/cbytecodes.cmo kernel/cbytegen.cmo kernel/environ.cmo \ - kernel/conv_oracle.cmo kernel/closure.cmo kernel/reduction.cmo \ - kernel/modops.cmo kernel/type_errors.cmo kernel/inductive.cmo \ - kernel/typeops.cmo kernel/subtyping.cmo kernel/indtypes.cmo \ - kernel/cooking.cmo \ - kernel/term_typing.cmo kernel/mod_typing.cmo kernel/safe_typing.cmo \ - library/summary.cmo library/global.cmo library/nameops.cmo \ - library/libnames.cmo library/nametab.cmo library/libobject.cmo \ - library/lib.cmo library/goptions.cmo library/decls.cmo library/heads.cmo \ - pretyping/termops.cmo pretyping/evd.cmo pretyping/rawterm.cmo \ - pretyping/reductionops.cmo pretyping/inductiveops.cmo \ - pretyping/retyping.cmo pretyping/cbv.cmo \ - pretyping/pretype_errors.cmo pretyping/typing.cmo pretyping/evarutil.cmo \ - pretyping/term_dnet.cmo pretyping/recordops.cmo \ - pretyping/evarconv.cmo pretyping/tacred.cmo \ - pretyping/classops.cmo pretyping/typeclasses_errors.cmo pretyping/typeclasses.cmo \ - pretyping/detyping.cmo pretyping/indrec.cmo pretyping/coercion.cmo \ - pretyping/unification.cmo pretyping/cases.cmo \ - pretyping/pretyping.cmo pretyping/clenv.cmo pretyping/pattern.cmo \ + $(addprefix kernel/, \ + names.cmo univ.cmo esubst.cmo term.cmo mod_subst.cmo sign.cmo \ + cbytecodes.cmo copcodes.cmo cemitcodes.cmo declarations.cmo \ + retroknowledge.cmo pre_env.cmo cbytegen.cmo conv_oracle.cmo \ + environ.cmo closure.cmo reduction.cmo type_errors.cmo entries.cmo \ + modops.cmo inductive.cmo typeops.cmo indtypes.cmo cooking.cmo \ + term_typing.cmo subtyping.cmo mod_typing.cmo safe_typing.cmo ) \ + $(addprefix library/, \ + summary.cmo global.cmo nameops.cmo libnames.cmo nametab.cmo \ + libobject.cmo lib.cmo goptions.cmo decls.cmo heads.cmo ) \ + $(addprefix pretyping/, \ + termops.cmo evd.cmo rawterm.cmo reductionops.cmo inductiveops.cmo \ + retyping.cmo cbv.cmo pretype_errors.cmo typing.cmo evarutil.cmo \ + term_dnet.cmo recordops.cmo evarconv.cmo tacred.cmo \ + classops.cmo typeclasses_errors.cmo typeclasses.cmo \ + detyping.cmo indrec.cmo coercion.cmo unification.cmo cases.cmo \ + pretyping.cmo clenv.cmo pattern.cmo ) \ parsing/lexer.cmo interp/ppextend.cmo interp/genarg.cmo \ interp/topconstr.cmo interp/notation.cmo interp/dumpglob.cmo interp/reserve.cmo \ library/impargs.cmo interp/constrextern.cmo \ @@ -840,6 +798,9 @@ NEWRINGVO:=$(addprefix contrib/setoid_ring/, \ Field_theory.vo Field_tac.vo \ Field.vo RealField.vo ) +GBVO:=$(addprefix contrib/groebner/, \ + GroebnerR.vo GroebnerZ.vo ) + XMLVO:= FOURIERVO:=$(addprefix contrib/fourier/, \ @@ -867,7 +828,7 @@ endif CONTRIBVO:= $(OMEGAVO) $(ROMEGAVO) $(MICROMEGAVO) $(RINGVO) $(FIELDVO) \ $(XMLVO) $(FOURIERVO) $(CCVO) $(FUNINDVO) \ $(RTAUTOVO) $(RECDEFVO) $(NEWRINGVO) $(DPVO) $(QUOTEVO) \ - $(INTERFACEVO) + $(INTERFACEVO) $(GBVO) ALLVO:= $(INITVO) $(THEORIESVO) $(CONTRIBVO) VFILES:= $(ALLVO:.vo=.v) diff --git a/contrib/groebner/GroebnerR.v b/contrib/groebner/GroebnerR.v new file mode 100644 index 0000000000..c64bc88159 --- /dev/null +++ b/contrib/groebner/GroebnerR.v @@ -0,0 +1,465 @@ +(************************************************************************) +(* 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 Reals. (* beware that Reals export Rlist *) +Require Import List. +Require Import ZArith. +Require Import Znumtheory. +Require Import List. +Require Import Ring_polynom. +Require Import Ring_tac. +Require Import Setoid. +Require Import BinPos. +Require Import BinList. +Require Import InitialRing. + +Declare ML Module "groebner". + +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. + +Ltac equalities_to_goal := + match goal with +| H: (@eq R ?x 0) |- _ => + try generalize H; clear H +| H: (@eq R 0 ?x) |- _ => + try generalize (sym_equal H); clear H +| H: (@eq R ?x ?y) |- _ => + try generalize (psos_r1 _ _ H); clear H + end. + +Lemma groebnerR_not1: forall x y:R, x<>y -> exists z:R, z*(x-y)-1=0. +intros. +exists (1/(x-y)). +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. + +Lemma groebnerR_not2: 1<>0. +auto with *. +Qed. + +Lemma groebnerR_diff: forall x y:R, x<>y -> x-y<>0. +Admitted. + +(* 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 => 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%R Rplus Rmult (gen_phiZ 0%R 1%R Rplus Rmult Ropp)). + +Definition PEevalR : list R -> PEZ -> R := + PEeval 0%R Rplus Rmult Rminus Ropp (gen_phiZ 0%R 1%R Rplus Rmult Ropp) + Nnat.nat_of_N pow. + +Lemma P0Z_correct : forall l, PhiR l P0Z = 0%R. +Proof. trivial. Qed. + + +Lemma PolZadd_correct : forall P' P l, + PhiR l (PolZadd P P') = (PhiR l P + PhiR l P')%R. +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')%R. +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%R /\ 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%R. +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%R. +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=0 -> p=0. +intros. +elim (Rmult_integral _ _ H0);intros. +absurd (c=0);auto. + 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. + +Ltac AddFvTail2 a l := + match l with + | (@nil _) => constr:(cons a l) + | (cons a _) => l + | (cons ?x ?l) => let l' := AddFvTail2 a l in constr:(cons x l') + 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 := AddFvTail2 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. + +(* Pompe sur Ring *) + +Ltac groebnerR_gen2 Cst_tac CstPow_tac lemma1 req n lH := + let Main lhs rhs R radd rmul rsub ropp rpow C := + let mkFV := FV Cst_tac CstPow_tac radd rmul rsub ropp rpow in + let mkPol := mkPolexpr C Cst_tac CstPow_tac radd rmul rsub ropp rpow in + match goal with + | Hparam: (@eq (@prod (@List.list ?R) (@List.list ?R)) + (?lparam, ?lvar) + ?lparam2) |- _ => + clear Hparam; + generalise_eq_hyps; + match goal with + | |- ?t => + let lpol:=lpol_goal t in + let rec fv_rec l lv:= + match l with + |?a::?l1 => + let lv1 := mkFV a lv in + fv_rec l1 lv1 + | nil => lv + end in + let fv := + match lvar with + | (@nil _) => let fv1 := FV_hypo_tac mkFV req lH in + let fv1 := fv_rec lpol fv1 in + rev fv1 + (* heuristique: les dernieres variables auront le poid le plus fort *) + | _ => lvar + end in + (*let fv1 := FV_hypo_tac mkFV req lH in + let fv := fv_rec lpol fv1 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 rec mkPol_list lp:= + match lp with + |?p::?lp1 => + let r := mkPol p fv in + let lr := mkPol_list lp1 in + constr:(r::lr) + | nil => constr:(@nil (@PExpr Z)) + end + in + let lpol := mkPol_list lpol in + let lpol := eval compute in (List.rev lpol) in + let lpol := constr:((@PEc Z nparam)::lpol) in + (*idtac lpol;*) + groebner_compute lpol; + match goal with + | |- ?res=?res1 -> _ => + let a1:= eval compute in lpol in + match a1 with + | ?np::?p2::?lp2 => + match res with + | (?p ::?lp0)::(?c::?lq)::?lci => + set (lci1:=lci); + set (lq1:=lq); + set (p21:=p2) ; + let q:= constr:(PEmul c p21) in + set (q1:=q); + set (lp21:=lp2); + let Hg := fresh "Hg" in + assert (Hg:check lp21 q1 (lci1,lq1) = true); + [vm_compute;trivial + |intros; + let Hg2 := fresh "Hg" in + assert (Hg2:(interpret3 q1 fv)=0); + [simpl; apply (@check_correct fv lp21 q1 (lci1,lq1) Hg); simpl; + repeat (match goal with h:_=0%R |- _ => rewrite h end); + repeat split; auto + |simpl in Hg2; simpl; + apply groebnerR_l3 with (interpret3 c fv);simpl; + [discrR; let Hc := fresh "Hg" in + unfold not; intro Hc; ring_simplify in Hc; + generalize Hc; clear Hc + |auto] + ] + ] + end + end + end + end + end + in + ParseRingComponents lemma1 ltac:(OnEquation req Main) + . + +Ltac groebnerR_gen + req sth ext morph arth cst_tac pow_tac lemma1 lemma2 pre post lH rl := + pre();groebnerR_gen2 cst_tac pow_tac lemma1 req ring_subst_niter lH . + +Ltac groebnerRpv lparam lvar:= + groebnerR_begin; + intros; + generalize (refl_equal (lparam,lvar)); intro; + let G := Get_goal in + ring_lookup groebnerR_gen [] G. + +Ltac groebnerR := groebnerRpv (@nil R) (@nil R). + +Ltac groebnerRp lparam := groebnerRpv lparam (@nil R). + + +(*************** Examples *) +(* +Open Scope R_scope. +Goal forall x y, + x+y=0 -> + x*y=0 -> + x^2=0. +Time groebnerR. +Qed. + +Goal forall x y z u v, + x+y+z+u+v=0 -> + x*y+x*z+x*u+x*v+y*z+y*u+y*v+z*u+z*v+u*v=0-> + 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. +Time groebnerR. +Qed. +*) + + + + + + + + + + + + + + + diff --git a/contrib/groebner/GroebnerZ.v b/contrib/groebner/GroebnerZ.v new file mode 100644 index 0000000000..ca58dbdb38 --- /dev/null +++ b/contrib/groebner/GroebnerZ.v @@ -0,0 +1,67 @@ +(************************************************************************) +(* 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 ZArith. +Require Import Zpow_facts. +Require Import Wf_Z. +Require Import Znumtheory. +Require Import List. +Require Export GroebnerR. + + +Open Scope Z_scope. + + +Lemma groebnerZhypR: forall x y:Z, x=y -> IZR x = IZR y. +Admitted. + +Lemma groebnerZconclR: forall x y:Z, IZR x = IZR y -> x = y. +Admitted. + +Lemma groebnerZhypnotR: forall x y:Z, x<>y -> IZR x <> IZR y. +Admitted. + +Lemma groebnerZconclnotR: forall x y:Z, IZR x <> IZR y -> x <> y. +Admitted. + +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. +Admitted. +Lemma groebnerZR2: forall x y:Z, IZR(x*y) = (IZR x * IZR y)%R. +Admitted. +Lemma groebnerZR3: forall x y:Z, IZR(x-y) = (IZR x - IZR y)%R. +Admitted. +Lemma groebnerZR4: forall (x:Z) p, IZR(x ^ Zpos p) = (IZR x ^ nat_of_P p)%R. +Admitted. +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/contrib/groebner/groebner.ml4 b/contrib/groebner/groebner.ml4 new file mode 100644 index 0000000000..895725f83d --- /dev/null +++ b/contrib/groebner/groebner.ml4 @@ -0,0 +1,450 @@ +(************************************************************************) +(* 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 +open Entiers +open Ideal + +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') + +(* ------------------------------------------------------------------------- *) +(* 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 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 + +(* 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 (Polynom.Pint (Polynom.coef_of_num r)) + | Var v -> + let v = int_of_string v in + if v <= np + then polconst d (Polynom.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 + |Polynom.Pint n -> Const (Polynom.num_of_coef n) + |Polynom.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 rec aux p = + match p with + |[] -> Const (num_of_string "0") + |(a,m)::p1 -> + let n = (Array.length m)-1 in + let (i0,e0) = + (List.fold_left + (fun (r,d) (a,m) -> + let i0= ref 0 in + for k=1 to n do + if m.(k)>0 + then i0:=k + done; + if !i0 = 0 + then (r,d) + else if !i0 > r + then (!i0, m.(!i0)) + else if !i0 = r && m.(!i0)<d + then (!i0, m.(!i0)) + else (r,d)) + (0,0) + p) in + if i0=0 + then + let mp = ref (polrec_to_term a) in + if p1=[] + then !mp + else Add(!mp,aux p1) + else ( + let p1=ref [] in + let p2=ref [] in + List.iter + (fun (a,m) -> + if m.(i0)>=e0 + then (m.(i0)<-m.(i0)-e0; + p1:=(a,m)::(!p1)) + else p2:=(a,m)::(!p2)) + p; + let vm = + if e0=1 + then Var (string_of_int (i0)) + else Pow (Var (string_of_int (i0)),e0) in + Add(Mul(vm, aux (List.rev (!p1))), aux (List.rev (!p2)))) + in let pp = aux p in + + pp + + +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 (c,lp0,p,lci,lc) = in_ideal m lpol p in + let lpc = List.rev !poldepcontent in + info ("temps: "^Format.sprintf "@[%10.3f@]s\n" (Unix.gettimeofday ()-.t1)); + if lc<>[] + then (c,lp0,p,lci,lc,1,lpc) + else (c,[],zeroP,[],[],0,[]) + + +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 (c,lp0,p,lci,lq,d,lct)= theoremedeszeros lpol p in + + let lci1 = List.rev lci in + match remove_zeros (fun x -> x=zeroP) (lq::lci1) with + | [] -> assert false + | (lq::lci) -> + (* lci commence par les nouveaux polynomes *) + let m= !nvars in + let lci = List.rev lci in + let c = pol_sparse_to_term m [c, Array.create (m+1) 0] in + let lp0 = List.map (pol_sparse_to_term m) lp0 in + let p = (Pow ((pol_sparse_to_term m p),d)) in + let lci = List.map (fun lc -> List.map (pol_sparse_to_term m) lc) lci in + let lq = List.map (pol_sparse_to_term m) lq in + info ("nombre de parametres: "^string_of_int nparam^"\n"); + info "terme calcule\n"; + (c,lp0,p,lci,lq) + ) + |_ -> assert false + + + + +let groebner lpol = + let lp= parse_request lpol in + let (c,lp0,p,lci,lq) = theoremedeszeros_termes lp in + let res = [p::lp0;c::lq]@lci in + let 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 + +(* +open Util +open Term +open Tactics +open Coqlib +open Utile +open Groebner +*) + +let groebner_compute t gl = + let lpol = groebner t in + let a = + mkApp(gen_constant "CC" ["Init";"Logic"] "refl_equal",[|tllp ();lpol|]) in + generalize [a] gl + +TACTIC EXTEND groebner_compute +| [ "groebner_compute" constr(lt) ] -> + [ groebner_compute lt ] +END + + diff --git a/contrib/groebner/ideal.ml b/contrib/groebner/ideal.ml new file mode 100644 index 0000000000..a87a9972bd --- /dev/null +++ b/contrib/groebner/ideal.ml @@ -0,0 +1,792 @@ +(************************************************************************) +(* 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 + +(*********************************************************************** + Coefficients: polynomes recursifs + on n'utilise pas les big_int car leur egalite n'est pas generique: on ne peut pas utiliser de tables de hash simplement. + *) + +type coef = Polynom.poly +let coef_of_int = Polynom.cf +let eq_coef = Polynom.eqP +let plus_coef =Polynom.plusP +let mult_coef = Polynom.multP +let sub_coef = Polynom.moinsP +let opp_coef = Polynom.oppP +let div_coef = Polynom.divP (* division exacte *) +let coef0 = Polynom.cf0 +let coef1 = Polynom.cf1 +let string_of_coef c = "["^(Polynom.string_of_P c)^"]" + + +let rec pgcd a b = Polynom.pgcdP a b + +(* +let htblpgcd = Hashtbl.create 1000 + +let pgcd a b = + try + (let c = Hashtbl.find htblpgcd (a,b) in + info "*"; + c) + with _ -> + (let c = Polynom.pgcdP a b in + info "-"; + Hashtbl.add htblpgcd (a,b) c; + c) +*) + +(*********************************************************************** + Monomes + tableau d'entiers + le premier coefficient est le degre + *) + +let hmon = Hashtbl.create 1000 + +type mon = int array + +(* 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) + + +(* Comparaison de monomes avec ordre du degre lexicographique = on commence par regarder la 1ere variable*) +(*let compare_mon d m m' = + let res=ref 0 in + let i=ref 1 in (* 1 si lexico pur 0 si degre*) + while (!res=0) && (!i<=d) do + res:= (compare m.(!i) m'.(!i)); + i:=!i+1; + done; + !res +*) + +(* degre lexicographique inverse *) +let compare_mon d m m' = + match compare m.(0) m'.(0) with + | 0 -> (* meme degre total *) + let res=ref 0 in + let i=ref d in + while (!res=0) && (!i>=1) do + res:= - (compare m.(!i) m'.(!i)); + i:=!i-1; + done; + !res + | x -> x + + +(* 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) -> (div_coef 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 + + *) + +type poly = (coef * mon) list + +let eq_poly = + Util.list_for_all2eq + (fun (c1,m1) (c2,m2) -> eq_coef c1 c2 && m1=m2) + +(* + 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=plus_coef (fst t) (fst t') in + match eq_coef 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 -> ((mult_coef 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 + match (div_mon_test d m m') with + true -> q + |false -> 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 + [((coef_of_int 1),m)] + + + +(* oppose d'un polynome *) +let oppP p = + let rec oppP p = + match p with + [] -> [] + |(b,m')::p -> ((opp_coef 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 -> ((mult_coef 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 rec contentP p = + match p with + |[] -> coef1 + |[a,m] -> a + |(a,m)::p1 -> pgcd a (contentP p1) + +let contentPlist lp = + match lp with + |[] -> coef1 + |p::l1 -> List.fold_left (fun r q -> pgcd r (contentP q)) (contentP p) l1 + +(*********************************************************************** + Division de polynomes + *) + +let pgcdpos a b = pgcd a b + + +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,[]) + |t::p' -> let (a,m)=t in + let q = (try Hashtbl.find hmon m + with Not_found -> + let q = selectdiv d m l in + match q with + t'::q' -> (Hashtbl.add hmon m q;q) + |[] -> q) in + match q with + [] -> if reduire_les_queues + then + let (c,r)=(reduce p') in + (c,((mult_coef a c,m)::r)) + else (coef1,p) + |(b,m')::q' -> + let c=(pgcdpos a b) in + let a'= (div_coef b c) in + let b'=(opp_coef (div_coef a c)) in + let (e,r)=reduce (div_pol d p' q' a' b' + (div_mon d m m')) in + (mult_coef a' e,r) + in let (c,r) = reduce p in + (c,r) + +(* trace des divisions *) + +(* liste des polynomes de depart *) +let poldep = ref [] +let poldepcontent = ref [] + +(* table des coefficients des polynomes en fonction des polynomes de depart *) +let coefpoldep = Hashtbl.create 51 + +(* coefficient de q dans l expression de p = sum_i c_i*q_i *) +let coefpoldep_find p q = + try (Hashtbl.find coefpoldep (p,q)) + with _ -> [] + +let coefpoldep_set p q c = + Hashtbl.add coefpoldep (p,q) c + +let initcoefpoldep d lp = + poldep:=lp; + poldepcontent:= List.map contentP (!poldep); + List.iter + (fun p -> coefpoldep_set p p (polconst d (coef_of_int 1))) + 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'=(opp_coef (div_coef 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 eq_poly s q + then + plusP d x (mult_t_pol d a m (polconst d (coef_of_int 1))) + else x) + c0 + lq in + c) + lcp + !poldep, + r) + +(*********************************************************************** + Completion + *) +let homogeneous = ref false + +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 fsp p' q' = + plusP d + (mult_t_pol d + (div_coef b c) + (div_mon d m'' m) p') + (mult_t_pol d + (opp_coef (div_coef a c)) + (div_mon d m'' m') q') in + let sp = fsp p' q' in + coefpoldep_set sp p (fsp (polconst d (coef_of_int 1)) []); + coefpoldep_set sp q (fsp [] (polconst d (coef_of_int 1))); + 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 list_rec f0 f1 = + let rec f2 = function + [] -> f0 + | a0::l0 -> f1 a0 l0 (f2 l0) + in f2 + +let addRes i = function + Keep h'0 -> Keep (i::h'0) + | DontKeep h'0 -> DontKeep (i::h'0) + +let slice d i a q = + list_rec + (match etrangers d i a with + true -> DontKeep [] + | false -> Keep []) + (fun b q1 rec_ren -> + match div_ppcm d i a b with + true -> DontKeep (b::q1) + | false -> + (match div_ppcm d i b a with + true -> rec_ren + | false -> addRes b rec_ren)) q + +let rec addS x l = l @[x] + +(* ajoute les spolynomes de i avec la liste de polynomes aP, + a la liste q *) + +let genPcPf d i aP q = + (let rec genPc aP0 = + match aP0 with + [] -> (fun r -> r) + | a::l1 -> + (fun l -> + (match slice d i a l1 with + Keep l2 -> addS (spol d i a) (genPc l2 l) + | DontKeep l2 -> genPc l2 l)) + in genPc aP) q + +let genOCPf d h' = + list_rec [] (fun a l rec_ren -> + genPcPf d a l rec_ren) h' + +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 +let pol_courant = ref [] + +let test_dans_ideal d p lp lp0 = + let (c,r) = reduce2 d !pol_courant lp in + info ("reste: "^(stringPcut r)^"\n"); + coef_courant:= mult_coef !coef_courant c; + 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"); + + let rec aux lp = + match lp with + |[] -> [] + |p::lp -> + (List.map + (fun q -> coefpoldep_find p q) + lp)::(aux lp) + in + let coefficient_multiplicateur = c in + let liste_polynomes_de_depart = List.rev lp0 in + let polynome_a_tester = p in + let liste_des_coefficients_intermediaires = + (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 (coef_of_int (-1)) cq) + (List.rev lcq) in + (coefficient_multiplicateur, + liste_polynomes_de_depart, + polynome_a_tester, + liste_des_coefficients_intermediaires, + liste_des_coefficients)) + else ((*info "polynome non reduit a 0\n"; + info ("\nreste: "^(stringPcut r)^"\n");*) + (c,[],[],[],[])) + +let choix_spol p l = l + +let deg_hom p = + match p with + | [] -> -1 + | (a,m)::_ -> m.(0) + +let pbuchf d pq p lp0= + info "calcul de la base de Groebner\n"; + step:=0; + Hashtbl.clear hmon; + let rec pbuchf x = + let (lp, lpc) = x in + infobuch lp lpc; +(* step:=(!step+1)mod 10;*) + match lpc with + [] -> test_dans_ideal d p lp lp0 + | _ -> + (match choix_spol !pol_courant lpc with + |[] -> assert false + |(a::lpc2) -> + if !homogeneous && a<>[] && deg_hom a > deg_hom !pol_courant + then (info "h";pbuchf (lp, lpc2)) + else + let (ca,a0)= reduce2 d a lp in + match a0 with + [] -> info "0";pbuchf (lp, lpc2) + | _ -> + let a1 = emultP ca a in + List.iter + (fun q -> + coefpoldep_set a1 q (emultP ca (coefpoldep_find a q))) + !poldep; + let (lca,a0) = reduce2_trace d a1 lp + (List.map (fun q -> coefpoldep_find a1 q) !poldep) in + List.iter2 (fun c q -> coefpoldep_set a0 q c) lca !poldep; + info ("\nnouveau polynome: "^(stringPcut a0)^"\n"); + let ct = coef1 (* contentP a0 *) in + (*info ("contenu: "^(string_of_coef ct)^"\n");*) + poldep:=addS a0 lp; + poldepcontent:=addS ct (!poldepcontent); + + let (c1,lp1,p1,lci1,lc1) + = test_dans_ideal d p (addS a0 lp) lp0 in + if lc1<>[] + then (c1,lp1,p1,lci1,lc1) + else + pbuchf + (((addS a0 lp), + (genPcPf d a0 lp lpc2)))) + in pbuchf pq + +let 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 = + info ("p: "^(stringPcut p)^"\n"); + info ("lp:\n"^(List.fold_left (fun r p -> r^(stringPcut p)^"\n") "" lp)); + initcoefpoldep d lp; + coef_courant:=coef1; + pol_courant:=p; + homogeneous := List.for_all is_homogeneous (p::lp); + if !homogeneous then info "polynomes homogenes\n"; + let (c1,lp1,p1,lci1,lc1) = test_dans_ideal d p lp lp in + if lc1<>[] + then (c1,lp1,p1,lci1,lc1) + else pbuchf d (lp, (genOCPf d lp)) p lp + + diff --git a/contrib/groebner/ideal.mli b/contrib/groebner/ideal.mli new file mode 100644 index 0000000000..48b81102e5 --- /dev/null +++ b/contrib/groebner/ideal.mli @@ -0,0 +1,82 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(************************************************************************) + +type coef = Polynom.poly +val coef_of_int : int -> coef +val eq_coef : coef -> coef -> bool +val plus_coef : coef -> coef -> coef +val mult_coef : coef -> coef -> coef +val sub_coef : coef -> coef -> coef +val opp_coef : coef -> coef +val div_coef : coef -> coef -> coef +val coef0 : coef +val coef1 : coef +val string_of_coef : coef -> string +val pgcd : coef -> coef -> coef + +type mon = int array +type poly = (coef * mon) list + +val eq_poly : poly -> poly -> bool +val mult_mon : int -> mon -> mon -> mon +val deg : mon -> int +val compare_mon : int -> mon -> mon -> int +val div_mon : int -> mon -> mon -> mon +val div_pol_coef : poly -> coef -> poly +val div_mon_test : int -> mon -> mon -> bool +val set_deg : int -> mon -> mon +val ppcm_mon : int -> mon -> mon -> mon + +val name_var : string list ref +val getvar : string list -> int -> string +val stringP : poly -> string +val stringPcut : poly -> string +val lstringP : poly list -> string +val printP : poly -> unit +val lprintP : poly list -> unit + +val polconst : int -> coef -> poly +val zeroP : poly +val plusP : int -> poly -> poly -> poly +val mult_t_pol : int -> coef -> mon -> poly -> poly +val selectdiv : int -> mon -> poly list -> poly +val gen : int -> int -> poly +val oppP : poly -> poly +val emultP : coef -> poly -> poly +val multP : int -> poly -> poly -> poly +val puisP : int -> poly -> int -> poly +val contentP : poly -> coef +val contentPlist : poly list -> coef +val pgcdpos : coef -> coef -> coef +val div_pol : int -> poly -> poly -> coef -> coef -> mon -> poly +val reduce2 : int -> poly -> poly list -> coef * poly + +val poldepcontent : coef list ref +val coefpoldep_find : poly -> poly -> poly +val coefpoldep_set : poly -> poly -> poly -> unit +val initcoefpoldep : int -> poly list -> unit +val reduce2_trace : + int -> poly -> poly list -> poly list -> poly list * poly +val homogeneous : bool ref +val spol : int -> poly -> poly -> poly +val etrangers : int -> ('a * mon) list -> ('b * mon) list -> bool +val div_ppcm : int -> poly -> poly -> poly -> bool + +type 'a cpRes = Keep of 'a list | DontKeep of 'a list +val list_rec : 'a -> ('b -> 'b list -> 'a -> 'a) -> 'b list -> 'a +val addRes : 'a -> 'a cpRes -> 'a cpRes +val slice : int -> poly -> poly -> poly list -> poly cpRes +val addS : 'a -> 'a list -> 'a list +val genPcPf : int -> poly -> poly list -> poly list -> poly list +val genOCPf : int -> poly list -> poly list + +val is_homogeneous : poly -> bool + +val in_ideal : + int -> poly list -> poly -> + coef * poly list * poly * poly list list * poly list diff --git a/contrib/groebner/polynom.ml b/contrib/groebner/polynom.ml new file mode 100644 index 0000000000..13acaa240d --- /dev/null +++ b/contrib/groebner/polynom.ml @@ -0,0 +1,1055 @@ +(************************************************************************) + +(* 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 + +(*********************************************************************** + 1. Opérations sur les coefficients. +*) + +open Big_int + +type coef = big_int +let coef_of_int = big_int_of_int +let coef_of_num = Num.big_int_of_num +let num_of_coef = Num.num_of_big_int +let eq_coef = eq_big_int +let lt_coef = lt_big_int +let le_coef = le_big_int +let abs_coef = abs_big_int +let plus_coef =add_big_int +let mult_coef = mult_big_int +let sub_coef = sub_big_int +let opp_coef = minus_big_int +let div_coef = div_big_int +let mod_coef = mod_big_int +let coef0 = coef_of_int 0 +let coef1 = coef_of_int 1 +let string_of_coef = string_of_big_int +let int_of_coef x = int_of_big_int x +let hash_coef x = + try (int_of_big_int x) + with _-> 1 +let puis_coef = power_big_int_positive_int + +(* +type coef = Entiers.entiers +let coef_of_int = Entiers.ent_of_int +let coef_of_num x = (* error if x is a ratio *) + Entiers.ent_of_string + (Big_int.string_of_big_int (Num.big_int_of_num x)) +let num_of_coef x = Num.num_of_string (Entiers.string_of_ent x) +let eq_coef = Entiers.eq_ent +let lt_coef = Entiers.lt_ent +let le_coef = Entiers.le_ent +let abs_coef = Entiers.abs_ent +let plus_coef =Entiers.add_ent +let mult_coef = Entiers.mult_ent +let sub_coef = Entiers.moins_ent +let opp_coef = Entiers.opp_ent +let div_coef = Entiers.div_ent +let mod_coef = Entiers.mod_ent +let coef0 = Entiers.ent0 +let coef1 = Entiers.ent1 +let string_of_coef = Entiers.string_of_ent +let int_of_coef x = Entiers.int_of_ent x +let hash_coef x =Entiers.hash_ent x +let signe_coef = Entiers.signe_ent + +let rec puis_coef p n = match n with + 0 -> coef1 + |_ -> (mult_coef p (puis_coef p (n-1))) +*) + +(* a et b positifs, résultat positif *) +let rec pgcd a b = + if eq_coef b coef0 + then a + else if lt_coef a b then pgcd b a else pgcd b (mod_coef a b) + + +(* signe du pgcd = signe(a)*signe(b) si non nuls. *) +let pgcd2 a b = + if eq_coef a coef0 then b + else if eq_coef b coef0 then a + else let c = pgcd (abs_coef a) (abs_coef b) in + if ((lt_coef coef0 a)&&(lt_coef b coef0)) + ||((lt_coef coef0 b)&&(lt_coef a coef0)) + then opp_coef c else c + +(*********************************************************************** + 2. Le type des polynômes, operations. +*) + +type variable = int + +type poly = + Pint of coef (* polynome constant *) + | Prec of variable * (poly array) (* coefficients par degre croissant *) + +(* 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 cf x = Pint (coef_of_int x) +let cf0 = (cf 0) +let cf1 = (cf 1) + +(* la n-ième variable *) +let x n = Prec (n,[|cf0;cf1|]) + +(* 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 eqP p q = + match (p,q) with + (Pint a,Pint b) -> eq_coef a b + |(Prec(x,p1),Prec(y,q1)) -> + if x<>y then false + else if (Array.length p1)<>(Array.length q1) then false + else (try (Array.iteri (fun i a -> if not (eqP a q1.(i)) + 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 && (eqP 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 (plus_coef 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 -> abs_coef a + | Prec (x ,p1) -> + Array.fold_left pgcd coef0 (Array.map content p1) + + +(* divise tous les coefficients de p par l'entier a*) +let rec div_int p a= + match p with + Pint b -> Pint (div_coef b a) + | Prec(x,p1) -> Prec(x,Array.map (fun x -> div_int x a) p1) + +(* divise p par son contenu entier positif. *) +let vire_contenu p = + let c = content p in + if eq_coef c coef0 then p else div_int p c + + +(* 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)))) + + +(* conversion d'un poly cst en entier*) +let int_of_Pint=function + (Pint x)->x + |_->failwith "non" + + +(* multiplie p par v^n, v >= max_var p *) +let rec multx n v p = + match p with + 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 (mult_coef a b) + |(Pint a, Prec (y,q1)) -> + if eq_coef a coef0 then Pint coef0 + else let q2 = Array.map (fun z-> multP p z) q1 in + Prec (y,q2) + + |(Prec (x,p1), Pint b) -> + if eq_coef b coef0 then Pint coef0 + 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 (opp_coef 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 le_coef coef0 a then p else oppP p + + +(* teste si un polynome est constant *) +let is_constantP p = match p with + Pint _ -> true + |Prec (_,_) -> false + + +(* teste si un poly est identiquement nul *) +let is_zero p = + match p with Pint n -> if eq_coef n coef0 then true else false |_-> false + +(* crée rapidement v^n *) +let monome v n = + match n with + 0->Pint coef1; + |_->let tmp = Array.create (n+1) (Pint coef0) in + tmp.(n)<-(Pint coef1); + Prec (v, tmp) + + +(*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 le_coef coef0 a + then string_of_coef a + else "("^(string_of_coef 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 string_of_P p = + nsP:=20; + string_of_Pcut p + +let printP p = Format.printf "@[%s@]" (string_of_P p) + + +let print_tpoly lp = + let s = ref "\n{ " in + Array.iter (fun p -> s:=(!s)^(string_of_P 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 eq_coef (mod_coef a b) coef0 + then (Pint (div_coef a b), cf 0) + 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 (eqP !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 eqP r cf0 + then s + else failwith ("div_pol:\n" + ^"p:"^(string_of_P p)^"\n" + ^"q:"^(string_of_P q)^"\n" + ^"r:"^(string_of_P r)^"\n" + ^"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 (pgcd (abs_coef a) (abs_coef b)) + |_ -> + if eqP p cf0 + then q + else if eqP q cf0 + then p + else if (deg x q) = 0 + then pgcd_coef_pol q p x + else if (deg x p) = 0 + then pgcd_coef_pol p q x + else ( + let a = content_pol p x in + let b = content_pol q x in + let c = pgcd_pol_rec a b (x-1) in + pr (string_of_int x); + let p1 = div_pol p c x in + let q1 = div_pol q c x in + let r = gcd_sub_res p1 q1 x in + let cr = content_pol r x in + let res = c @@ (div_pol r cr x) in + res + ) + +(* version avec memo*) +(* +and htblpgcd = Hashtbl.create 1000 + +and pgcd_pol_rec p q x = + try + (let c = Hashtbl.find htblpgcd (p,q,x) in + (*info "*";*) + c) + with _ -> + (let c = + (match (p,q) with + (Pint a,Pint b) -> Pint (pgcd (abs_coef a) (abs_coef b)) + |_ -> + if eqP p cf0 + then q + else if eqP q cf0 + then p + else if (deg x q) = 0 + then pgcd_coef_pol q p x + else if (deg x p) = 0 + then pgcd_coef_pol p q x + else ( + let a = content_pol p x in + let b = content_pol q x in + let c = pgcd_pol_rec a b (x-1) in + pr (string_of_int x); + let p1 = div_pol p c x in + let q1 = div_pol q c x in + let r = gcd_sub_res p1 q1 x in + let cr = content_pol r x in + let res = c @@ (div_pol r cr x) in + res + ) + ) + in + (*info "-";*) + Hashtbl.add htblpgcd (p,q,x) c; + c) +*) +(* Sous-résultants: + + ai*Ai = Qi*Ai+1 + bi*Ai+2 + + 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 eqP 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 eqP 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 + +let prcontentP () = + Hashtbl.iter (fun s c -> + prn (s^" -> "^(string_of_P c))) + hcontentP + + +let contentP = + memos "c" hcontentP (fun (p,x) -> ((string_of_P p)^":"^(string_of_int x))) + (fun (p,x) -> content_pol p x) + + +(* Tables de hash et polynômes, mémo *) + +(* dans le mli on écrit: + module Hashpol:Hashtbl.S with type t = poly +*) + +(* fonction de hachage des polynômes *) +let hashP p = + let rec hP p = + match p with + Pint a -> (hash_coef a) + |Prec (v,p) -> + (Array.fold_right (fun q h -> h+(hP q)) p 0) + in (* Hashtbl.hash *) (hP p) + + +module Hashpol = Hashtbl.Make( + struct + type t = poly + let equal = eqP + let hash = hashP + 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 ((string_of_P 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_select (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 -> eqP 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 -> le_coef coef0 a + |Prec(x,p1) -> + (array_test is_positif p1) + && (try (Array.iteri (fun i c -> if (i mod 2)<>0 && not (eqP c cf0) + then failwith "pas pair") + p1; + true) + with _ -> 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 (eqP verif cf0) + then (prn ("p:"^(string_of_P p)); + prn ("q:"^(string_of_P q)); + prn ("c:"^(string_of_P c)); + prn ("r:"^(string_of_P r)); + 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 eqP 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 (eqP 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: " + ^(string_of_P 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 eqP 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 eqP (lp@lf@[p]) in + r) + else lp + diff --git a/contrib/groebner/polynom.mli b/contrib/groebner/polynom.mli new file mode 100644 index 0000000000..99c1701b8f --- /dev/null +++ b/contrib/groebner/polynom.mli @@ -0,0 +1,111 @@ +type coef (*= Entiers.entiers*) +val coef_of_num : Num.num -> coef +val num_of_coef : coef -> Num.num +val eq_coef : coef -> coef -> bool +val lt_coef : coef -> coef -> bool +val le_coef : coef -> coef -> bool +val abs_coef : coef -> coef +val plus_coef : coef -> coef -> coef +val mult_coef : coef -> coef -> coef +val sub_coef : coef -> coef -> coef +val opp_coef : coef -> coef +val div_coef : coef -> coef -> coef +val mod_coef : coef -> coef -> coef +val coef0 : coef +val coef1 : coef +val string_of_coef : coef -> string +val int_of_coef : coef -> int +val hash_coef : coef -> int +val puis_coef : coef -> int -> coef +val pgcd : coef -> coef -> coef +val pgcd2 : coef -> coef -> coef + +type variable = int +type poly = Pint of coef | Prec of variable * poly array + +val cf : int -> poly +val cf0 : poly +val cf1 : poly +val x : variable -> poly +val max_var_pol : poly -> variable +val max_var_pol2 : poly -> variable +val max_var : poly array -> variable +val eqP : poly -> poly -> bool +val norm : poly -> poly +val deg : variable -> poly -> int +val deg_total : poly -> int +val copyP : poly -> poly +val coef : variable -> int -> poly -> poly +val plusP : poly -> poly -> poly +val content : poly -> coef +val div_int : poly -> coef -> poly +val vire_contenu : poly -> poly +val vars : poly -> variable list +val int_of_Pint : poly -> coef +val multx : int -> variable -> poly -> poly +val multP : poly -> poly -> poly +val deriv : variable -> poly -> poly +val oppP : poly -> poly +val moinsP : poly -> poly -> poly +val puisP : poly -> int -> poly +val ( @@ ) : poly -> poly -> poly +val ( -- ) : poly -> poly -> poly +val ( ^^ ) : poly -> int -> poly +val coefDom : variable -> poly -> poly +val coefConst : variable -> poly -> poly +val remP : variable -> poly -> poly +val coef_int_tete : poly -> coef +val normc : poly -> poly +val is_constantP : poly -> bool +val is_zero : poly -> bool +val monome : variable -> int -> poly +val coef_constant : poly -> coef +val univ : bool ref +val string_of_var : int -> string +val nsP : int ref +val string_of_Pcut : poly -> string +val string_of_P : poly -> string +val printP : poly -> unit +val print_tpoly : poly array -> unit +val print_lpoly : poly list -> unit +val quo_rem_pol : poly -> poly -> variable -> poly * poly +val div_pol : poly -> poly -> variable -> poly +val divP : poly -> poly -> poly +val div_pol_rat : poly -> poly -> bool +val pseudo_div : poly -> poly -> variable -> poly * poly * int * poly +val pgcdP : poly -> poly -> poly +val pgcd_pol : poly -> poly -> variable -> poly +val content_pol : poly -> variable -> poly +val pgcd_coef_pol : poly -> poly -> variable -> poly +val pgcd_pol_rec : poly -> poly -> variable -> poly +val gcd_sub_res : poly -> poly -> variable -> poly +val gcd_sub_res_rec : poly -> poly -> poly -> poly -> int -> variable -> poly +val lazard_power : poly -> poly -> int -> variable -> poly +val sans_carre : poly -> variable -> poly +val facteurs : poly -> variable -> poly list +val facteurs_impairs : poly -> variable -> poly list +val hcontentP : (string, poly) Hashtbl.t +val prcontentP : unit -> unit +val contentP : poly * variable -> poly +val hashP : poly -> int +module Hashpol : Hashtbl.S with type key=poly +val memoP : + string -> 'a Hashpol.t -> (poly -> 'a) -> poly -> 'a +val hfactorise : poly list list Hashpol.t +val prfactorise : unit -> unit +val factorise : poly -> poly list list +val facteurs2 : poly -> poly list +val pol_de_factorisation : poly list list -> poly +val set_of_array_facteurs : poly list array -> poly list +val factorise_tableauP2 : + poly array -> + poly list array -> poly array * (poly * int list) array +val factorise_tableauP : + poly array -> poly array * (poly * int list) array +val is_positif : poly -> bool +val is_negatif : poly -> bool +val pseudo_euclide : + poly list -> poly -> poly -> variable -> + poly * poly * int * poly * poly * (poly * int) list * (poly * int) list +val implique_non_nul : poly list -> poly -> bool +val ajoute_non_nul : poly -> poly list -> poly list diff --git a/contrib/groebner/utile.ml b/contrib/groebner/utile.ml new file mode 100644 index 0000000000..cc4d31978a --- /dev/null +++ b/contrib/groebner/utile.ml @@ -0,0 +1,228 @@ +(*********************************************************************** + 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 = + output_string stderr s; flush stderr + +(********************************************************************** + Tableaux +*) +(* compte les occurrences de i dans le tableau t *) +let compte t i = + let r = ref 0 in + Array.iter (fun x -> if x=i then r:= !r + 1) t; + !r ;; + +(* maximum d'un tableau d'entiers *) +let maximum t = + Array.fold_left max 0 t +;; + + +(* appliquer une fonction a tous les elements d'une matrice *) +let matrix_map f m= + let s=Array.length m in + let res=Array.create s [| |] in + for i=0 to s-1 do + res.(i)<-Array.map f m.(i) + done; +res +;; +(* selectionne dans un tableau *) +let array_select f l = + let res = ref [] in + Array.iter (fun x -> if (f x) then res:=(!res)@[x]) l; + Array.of_list !res +;; +(* teste si tous les elements d'un tableau verifient f *) +let array_test f t = + try (Array.iter (fun x -> if not (f x) then failwith "raté") t; + true) + with _ -> false +;; + +(* cherche a dans t, rend l'indice ou il se trouve ou declenche une exception*) +let array_find a t = + let n = Array.length t in + let ok = ref true in + let res = ref 0 in + while !ok do + if (!res = n) then raise Not_found +else (if t.(!res) = a then ok := false + else res := !res +1;) + done; +!res +;; + +(********************************************************************** + Listes +*) + +let set_of_list l = + let r = Hashtbl.create 51 in + List.iter (fun x -> + try (Hashtbl.find r x;()) + with _ -> Hashtbl.add r x true) l; + let res = ref [] in + Hashtbl.iter (fun x _ -> res:=x::(!res)) r; + !res + +(* 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 +;; + +(* selectionne dans une liste *) +let list_select f l = + let res = ref [] in + List.iter (fun x -> if (f x) then res:=(!res)@[x]) l; + !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_select (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/contrib/groebner/utile.mli b/contrib/groebner/utile.mli new file mode 100644 index 0000000000..927e6f4332 --- /dev/null +++ b/contrib/groebner/utile.mli @@ -0,0 +1,34 @@ + +(* Printing *) +val pr : string -> unit +val prn : string -> unit +val prt0 : 'a -> unit +val prt : string -> unit +val info : string -> unit + +(* Array *) +val compte : 'a array -> 'a -> int +val maximum : int array -> int +val matrix_map : ('a -> 'b) -> 'a array array -> 'b array array +val array_select : ('a -> bool) -> 'a array -> 'a array +val array_test : ('a -> bool) -> 'a array -> bool +val array_find : 'a -> 'a array -> int + +(* Listes *) +val set_of_list : 'a list -> 'a list +val list_mem_eq : ('a -> 'b -> bool) -> 'a -> 'b list -> bool +val set_of_list_eq : ('a -> 'a -> bool) -> 'a list -> 'a list +val list_select : ('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 |