aboutsummaryrefslogtreecommitdiffhomepage
path: root/plugins/nsatz
diff options
context:
space:
mode:
authorGravatar Pierre-Marie Pédrot <pierre-marie.pedrot@inria.fr>2017-02-24 12:23:43 +0100
committerGravatar Pierre-Marie Pédrot <pierre-marie.pedrot@inria.fr>2017-04-09 22:00:20 +0200
commitba636891121821b4943a0464b49e14de54de8253 (patch)
treed70caf013b7e428f376d48ea018d880b9d36a505 /plugins/nsatz
parent2c0287fe8445bd4b599bf8498bcb71b2a7df0d51 (diff)
Academic prescriptivism strikes back: down with baroque programming in Nsatz.
Several cleanups were performed. 1. Removal of dead code lurking around. 2. Removal of global variables used to pass arguments to functions, as well as unnecessary mutable state here and there. We rely on state-passing and encapsulated mutable state. 3. Removal of crazy reference manipulation and its replacement with proper list handling, as well as cleaning up the source and taking advantage of invariants. This should solve algorithmic limitations of the previous code. 4. Opacification of some structures to have a clearer idea of the code requirements. 5. Cleaning of debug printing functions. We thunk the computation of the debugging data, whose computation can be costly for no reason, and we rely on Feedback-based interaction instead of Printf-debugging.
Diffstat (limited to 'plugins/nsatz')
-rw-r--r--plugins/nsatz/ideal.ml702
-rw-r--r--plugins/nsatz/ideal.mli21
-rw-r--r--plugins/nsatz/nsatz.ml207
-rw-r--r--plugins/nsatz/utile.ml4
-rw-r--r--plugins/nsatz/utile.mli3
5 files changed, 356 insertions, 581 deletions
diff --git a/plugins/nsatz/ideal.ml b/plugins/nsatz/ideal.ml
index 48bdad826..f90e76386 100644
--- a/plugins/nsatz/ideal.ml
+++ b/plugins/nsatz/ideal.ml
@@ -23,7 +23,6 @@ exception NotInIdeal
Global options
*)
let lexico = ref false
-let use_hmon = ref false
(* division of tail monomials *)
@@ -33,31 +32,30 @@ let reduire_les_queues = false
let nouveaux_pol_en_tete = false
-(***********************************************************************
- Functor
-*)
-
-module Make (P:Polynom.S) = struct
-
- type coef = P.t
- let coef0 = P.of_num (Num.Int 0)
- let coef1 = P.of_num (Num.Int 1)
- let coefm1 = P.of_num (Num.Int (-1))
- let string_of_coef c = "["^(P.to_string c)^"]"
-
-(***********************************************************************
- Monomials
- array of integers, first is the degree
-*)
-
-type mon = int array
-type deg = int
-type poly = (coef * mon) list
-type polynom =
- {pol : poly ref;
- num : int;
- sugar : int}
-
+type metadata = {
+ name_var : string list;
+}
+
+module Monomial :
+sig
+type t
+val repr : t -> int array
+val make : int array -> t
+val deg : t -> int
+val nvar : t -> int
+val var_mon : int -> int -> t
+val mult_mon : t -> t -> t
+val compare_mon : t -> t -> int
+val div_mon : t -> t -> t
+val div_mon_test : t -> t -> bool
+val ppcm_mon : t -> t -> t
+val const_mon : int -> t
+end =
+struct
+type t = int array
+type mon = t
+let repr m = m
+let make m = m
let nvar (m : mon) = Array.length m - 1
let deg (m : mon) = m.(0)
@@ -104,9 +102,6 @@ let div_mon m m' =
done;
m''
-let div_pol_coef p c =
- List.map (fun (a,m) -> (P.divP a c,m)) p
-
(* m' divides m *)
let div_mon_test m m' =
let d = nvar m in
@@ -135,7 +130,46 @@ let ppcm_mon m m' =
done;
set_deg m''
+(* returns a constant polynom ial with d variables *)
+let const_mon d =
+ let m = Array.make (d+1) 0 in
+ let m = set_deg m in
+ m
+
+let var_mon d i =
+ let m = Array.make (d+1) 0 in
+ m.(i) <- 1;
+ let m = set_deg m in
+ m
+
+end
+
+(***********************************************************************
+ Functor
+*)
+
+module Make (P:Polynom.S) = struct
+ type coef = P.t
+ let coef0 = P.of_num (Num.Int 0)
+ let coef1 = P.of_num (Num.Int 1)
+ let coefm1 = P.of_num (Num.Int (-1))
+ let string_of_coef c = "["^(P.to_string c)^"]"
+
+(***********************************************************************
+ Monomials
+ array of integers, first is the degree
+*)
+
+open Monomial
+
+type mon = Monomial.t
+type deg = int
+type poly = (coef * mon) list
+type polynom = {
+ pol : poly;
+ num : int;
+}
(**********************************************************************
Polynomials
@@ -179,8 +213,8 @@ let string_of_pol zeroP hdP tlP coefterm monterm string_of_coef
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)]);
+ | "1" -> s:= (!s) @ [(getvar lvar (i-1))]
+ | e -> s:= (!s) @ [((getvar lvar (i-1)) ^ "^" ^ e)]);
done;
(match !s with
[] -> if coefone
@@ -228,8 +262,8 @@ let print_pol zeroP hdP tlP coefterm monterm string_of_coef
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)]);
+ | "1" -> s:= (!s) @ [(getvar lvar (i-1))]
+ | e -> s:= (!s) @ [((getvar lvar (i-1)) ^ "^" ^ e)]);
done;
(match !s with
[] -> if coefone
@@ -271,9 +305,7 @@ let print_pol zeroP hdP tlP coefterm monterm string_of_coef
print_flush()
-let name_var= ref []
-
-let stringP p =
+let stringP metadata (p : poly) =
string_of_pol
(fun p -> match p with [] -> true | _ -> false)
(fun p -> match p with (t::p) -> t |_ -> failwith "print_pol dans dansideal")
@@ -281,55 +313,29 @@ let stringP p =
(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
+ (fun m -> (Array.length (Monomial.repr m))-1)
+ (fun m i -> (string_of_int ((Monomial.repr m).(i))))
+ metadata.name_var
p
-let nsP2 = ref max_int
+let nsP2 = 10
-let stringPcut p =
+let stringPcut metadata (p : poly) =
(*Polynomesrec.nsP1:=20;*)
- nsP2:=10;
let res =
- if (List.length p)> !nsP2
- then (stringP [List.hd p])^" + "^(string_of_int (List.length p))^" terms"
- else stringP p in
+ if (List.length p)> nsP2
+ then (stringP metadata [List.hd p])^" + "^(string_of_int (List.length p))^" terms"
+ else stringP metadata p in
(*Polynomesrec.nsP1:= max_int;*)
- nsP2:= max_int;
res
-let rec lstringP l =
- match l with
- [] -> ""
- |p::l -> (stringP p)^("\n")^(lstringP l)
-
-let printP = print_pol
- (fun p -> match p with [] -> true | _ -> false)
- (fun p -> match p with (t::p) -> t |_ -> failwith "print_pol dans dansideal")
- (fun p -> match p with (t::p) -> p |_ -> failwith "print_pol dans dansideal")
- (fun (a,m) -> a)
- (fun (a,m) -> m)
- string_of_coef
- (fun m -> (Array.length m)-1)
- (fun m i -> (string_of_int (m.(i))))
- name_var
-
-
-let rec lprintP l =
- match l with
- [] -> ()
- |p::l -> printP p;print_string "\n"; lprintP l
-
-
(* Operations *)
let zeroP = []
(* returns a constant polynom ial with d variables *)
let polconst d c =
- let m = Array.make (d+1) 0 in
- let m = set_deg m in
+ let m = const_mon d in
[(c,m)]
let plusP p q =
@@ -357,9 +363,7 @@ let coef_of_int x = P.of_num (Num.Int x)
(* variable i *)
let gen d i =
- let m = Array.make (d+1) 0 in
- m.(i) <- 1;
- let m = set_deg m in
+ let m = var_mon d i in
[((coef_of_int 1),m)]
let oppP p =
@@ -390,7 +394,7 @@ let puisP p n=
|_ ->
if n = 0 then
let d = nvar (snd (List.hd p)) in
- [coef1, Array.make (d+1) 0]
+ [coef1, const_mon d]
else
let rec puisP p n =
if n = 1 then p
@@ -400,49 +404,34 @@ let puisP p n=
if n mod 2 = 0 then q else multP p q
in puisP p n
-let rec contentP p =
- match p with
- |[] -> coef1
- |[a,m] -> a
- |(a,m)::p1 ->
- if P.equal a coef1 || P.equal a coefm1
- then a
- else P.pgcdP a (contentP p1)
-
-let contentPlist lp =
- match lp with
- |[] -> coef1
- |p::l1 ->
- List.fold_left
- (fun r q ->
- if P.equal r coef1 || P.equal r coefm1
- then r
- else P.pgcdP r (contentP q))
- (contentP p) l1
-
(***********************************************************************
Division of polynomials
*)
+type table = {
+ hmon : (mon, poly) Hashtbl.t option;
+ (* coefficients of polynomials when written with initial polynomials *)
+ coefpoldep : ((int * int), poly) Hashtbl.t;
+ mutable nallpol : int;
+ mutable allpol : polynom array;
+ (* list of initial polynomials *)
+}
+
let pgcdpos a b = P.pgcdP a b
-let polynom0 = {pol = ref []; num = 0; sugar = 0}
+let polynom0 = { pol = []; num = 0 }
-let ppol p = !(p.pol)
+let ppol p = p.pol
let lm p = snd (List.hd (ppol p))
-let nallpol = ref 0
-
-let allpol = ref (Array.make 1000 polynom0)
-
-let new_allpol p s =
- nallpol := !nallpol + 1;
- if !nallpol >= Array.length !allpol
+let new_allpol table p =
+ table.nallpol <- table.nallpol + 1;
+ if table.nallpol >= Array.length table.allpol
then
- allpol := Array.append !allpol (Array.make !nallpol polynom0);
- let p = {pol = ref p; num = !nallpol; sugar = s} in
- !allpol.(!nallpol)<- p;
+ table.allpol <- Array.append table.allpol (Array.make table.nallpol polynom0);
+ let p = { pol = p; num = table.nallpol } in
+ table.allpol.(table.nallpol) <- p;
p
(* returns a polynomial of l whose head monomial divides m, else [] *)
@@ -456,43 +445,42 @@ let rec selectdiv m l =
|false -> selectdiv m r
let div_pol p q a b m =
-(* info ".";*)
plusP (emultP a p) (mult_t_pol b m q)
-let hmon = Hashtbl.create 1000
-
-let use_hmon = ref false
-
-let find_hmon m =
- if !use_hmon
- then Hashtbl.find hmon m
- else raise Not_found
-
-let add_hmon m q =
- if !use_hmon
- then Hashtbl.add hmon m q
- else ()
+let find_hmon table m = match table.hmon with
+| None -> raise Not_found
+| Some hmon -> Hashtbl.find hmon m
+
+let add_hmon table m q =
+match table.hmon with
+| None -> ()
+| Some hmon -> Hashtbl.add hmon m q
+
+let selectdiv table m l =
+ try find_hmon table m
+ with Not_found ->
+ let q = selectdiv m l in
+ let q = ppol q in
+ match q with
+ | [] -> q
+ | _ :: _ ->
+ let () = add_hmon table m q in
+ q
let div_coef a b = P.divP a b
(* remainder r of the division of p by polynomials of l, returns (c,r) where c is the coefficient for pseudo-division : c p = sum_i q_i p_i + r *)
-let reduce2 p l =
+let reduce2 table p l =
let l = if nouveaux_pol_en_tete then List.rev l else l in
let rec reduce p =
match p with
[] -> (coef1,[])
|t::p' ->
let (a,m)=t in
- let q = (try find_hmon m
- with Not_found ->
- let q = selectdiv m l in
- match (ppol q) with
- t'::q' -> (add_hmon m q;
- q)
- |[] -> q) in
- match (ppol q) with
+ let q = selectdiv table m l in
+ match q with
[] -> if reduire_les_queues
then
let (c,r)=(reduce p') in
@@ -508,37 +496,19 @@ let reduce2 p l =
in let (c,r) = reduce p in
(c,r)
-(* trace of divisions *)
-
-(* list of initial polynomials *)
-let poldep = ref []
-let poldepcontent = ref []
-
-(* coefficients of polynomials when written with initial polynomials *)
-let coefpoldep = Hashtbl.create 51
-
(* coef of q in p = sum_i c_i*q_i *)
-let coefpoldep_find p q =
- try (Hashtbl.find coefpoldep (p.num,q.num))
+let coefpoldep_find table p q =
+ try (Hashtbl.find table.coefpoldep (p.num,q.num))
with Not_found -> []
-let coefpoldep_remove p q =
- Hashtbl.remove coefpoldep (p.num,q.num)
-
-let coefpoldep_set p q c =
- Hashtbl.add coefpoldep (p.num,q.num) c
-
-let initcoefpoldep d lp =
- poldep:=lp;
- poldepcontent:= List.map (fun p -> contentP (ppol p)) lp;
- List.iter
- (fun p -> coefpoldep_set p p (polconst d (coef_of_int 1)))
- lp
+let coefpoldep_set table p q c =
+ Hashtbl.add table.coefpoldep (p.num,q.num) c
(* keeps trace in coefpoldep
divides without pseudodivisions *)
-let reduce2_trace p l lcp =
+let reduce2_trace table p l lcp =
+ let lp = l in
let l = if nouveaux_pol_en_tete then List.rev l else l in
(* rend (lq,r), ou r = p + sum(lq) *)
let rec reduce p =
@@ -546,15 +516,8 @@ let reduce2_trace p l lcp =
[] -> ([],[])
|t::p' ->
let (a,m)=t in
- let q =
- (try find_hmon m
- with Not_found ->
- let q = selectdiv m l in
- match (ppol q) with
- t'::q' -> (add_hmon m q;
- q)
- |[] -> q) in
- match (ppol q) with
+ let q = selectdiv table m l in
+ match q with
[] ->
if reduire_les_queues
then
@@ -568,19 +531,12 @@ let reduce2_trace p l lcp =
let (lq,r)=reduce p1 in
((b',m'',q)::lq, r)
in let (lq,r) = reduce p in
- (*info "reduce2_trace:\n";
- iter
- (fun (a,m,s) ->
- let x = mult_t_pol a m s in
- info ((stringP x)^"\n"))
- lq;
- info "ok\n";*)
(List.map2
(fun c0 q ->
let c =
List.fold_left
(fun x (a,m,s) ->
- if equal (ppol s) (ppol q)
+ if equal s (ppol q)
then
plusP x (mult_t_pol a m (polconst (nvar m) (coef_of_int 1)))
else x)
@@ -588,25 +544,14 @@ let reduce2_trace p l lcp =
lq in
c)
lcp
- !poldep,
+ lp,
r)
-let homogeneous = ref false
-let pol_courant = ref polynom0
-
(***********************************************************************
Completion
*)
-let sugar_flag = ref true
-
-let compute_sugar p =
- List.fold_left (fun s (a,m) -> max s m.(0)) 0 p
-
-let mk_polynom p =
- new_allpol p (compute_sugar p)
-
-let spol ps qs=
+let spol0 ps qs=
let p = ppol ps in
let q = ppol qs in
let m = snd (List.hd p) in
@@ -628,14 +573,9 @@ let spol ps qs=
(P.oppP (div_coef a c))
m2 q') in
let sp = fsp p' q' in
- let sps =
- new_allpol
- sp
- (max (m1.(0) + ps.sugar) (m2.(0) + qs.sugar)) in
- coefpoldep_set sps ps (fsp (polconst (nvar m) (coef_of_int 1)) []);
- coefpoldep_set sps qs (fsp [] (polconst (nvar m) (coef_of_int 1)));
- sps
-
+ let p0 = fsp (polconst (nvar m) (coef_of_int 1)) [] in
+ let q0 = fsp [] (polconst (nvar m) (coef_of_int 1)) in
+ (sp, p0, q0)
let etrangers p p'=
let m = snd (List.hd p) in
@@ -644,118 +584,34 @@ let etrangers p p'=
let res=ref true in
let i=ref 1 in
while (!res) && (!i<=d) do
- res:= (m.(!i) = 0) || (m'.(!i)=0);
+ res:= ((Monomial.repr m).(!i) = 0) || ((Monomial.repr m').(!i)=0);
i:=!i+1;
done;
!res
-(* teste if head monomial of p'' divides lcm of lhead monomials of p and p' *)
-
-let div_ppcm p p' p'' =
- let m = snd (List.hd p) in
- let m'= snd (List.hd p') in
- let m''= snd (List.hd p'') in
- let d = nvar m in
- let res=ref true in
- let i=ref 1 in
- while (!res) && (!i<=d) do
- res:= ((max m.(!i) m'.(!i)) >= m''.(!i));
- i:=!i+1;
- done;
- !res
-
-(* code from extraction of Laurent Théry Coq program *)
-
-type 'poly cpRes =
- Keep of ('poly list)
- | DontKeep of ('poly list)
-
-let list_rec f0 f1 =
- let rec f2 = function
- [] -> f0
- | a0::l0 -> f1 a0 l0 (f2 l0)
- in f2
-
-let addRes i = function
- Keep h'0 -> Keep (i::h'0)
- | DontKeep h'0 -> DontKeep (i::h'0)
-
-let slice i a q =
- list_rec
- (match etrangers (ppol i) (ppol a) with
- true -> DontKeep []
- | false -> Keep [])
- (fun b q1 rec_ren ->
- match div_ppcm (ppol i) (ppol a) (ppol b) with
- true -> DontKeep (b::q1)
- | false ->
- (match div_ppcm (ppol i) (ppol b) (ppol a) with
- true -> rec_ren
- | false -> addRes b rec_ren)) q
-
-(* sugar strategy *)
-
let addS x l = l @ [x] (* oblige de mettre en queue sinon le certificat deconne *)
-
-let addSsugar x l =
- if !sugar_flag
- then
- let sx = x.sugar in
- let rec insere l =
- match l with
- | [] -> [x]
- | y::l1 ->
- if sx <= y.sugar
- then x::l
- else y::(insere l1)
- in insere l
- else addS x l
-
-(* ajoute les spolynomes de i avec la liste de polynomes aP,
- a la liste q *)
-
-let genPcPf i aP q =
- (let rec genPc aP0 =
- match aP0 with
- [] -> (fun r -> r)
- | a::l1 ->
- (fun l ->
- (match slice i a l1 with
- Keep l2 -> addSsugar (spol i a) (genPc l2 l)
- | DontKeep l2 -> genPc l2 l))
- in genPc aP) q
-
-let genOCPf h' =
- list_rec [] (fun a l rec_ren ->
- genPcPf a l rec_ren) h'
-
+
(***********************************************************************
critical pairs/s-polynomials
*)
let ordcpair ((i1,j1),m1) ((i2,j2),m2) =
-(* let s1 = (max
- (!allpol.(i1).sugar + m1.(0)
- - (snd (hd (ppol !allpol.(i1)))).(0))
- (!allpol.(j1).sugar + m1.(0)
- - (snd (hd (ppol !allpol.(j1)))).(0))) in
- let s2 = (max
- (!allpol.(i2).sugar + m2.(0)
- - (snd (hd (ppol !allpol.(i2)))).(0))
- (!allpol.(j2).sugar + m2.(0)
- - (snd (hd (ppol !allpol.(j2)))).(0))) in
- match compare s1 s2 with
- | 1 -> 1
- |(-1) -> -1
- |0 -> compare_mon m1 m2*)
-
compare_mon m1 m2
let sortcpairs lcp =
List.sort ordcpair lcp
let mergecpairs l1 l2 =
- List.merge ordcpair l1 l2
+ let rec merge l1 l2 accu = match l1, l2 with
+ | [], [] -> List.rev accu
+ | l, [] | [], l -> List.rev_append accu l
+ | x1 :: t1, x2 :: t2 ->
+ if ordcpair x1 x2 <= 0 then
+ merge t1 l2 (x1 :: accu)
+ else
+ merge l1 t2 (x2 :: accu)
+ in
+ merge l1 l2 []
let ord i j =
if i<j then (i,j) else (j,i)
@@ -775,170 +631,140 @@ let cpairs lp =
[]|[_] -> []
|p::l1 -> mergecpairs (cpairs1 p l1) (aux l1)
in aux lp
-
-
-let critere2 ((i,j),m) lp lcp =
- List.exists
- (fun h ->
- h.num <> i && h.num <> j
- && (div_mon_test m (lm h))
- && (let c1 = ord i h.num in
- not (List.exists (fun (c,_) -> c1 = c) lcp))
- && (let c1 = ord j h.num in
- not (List.exists (fun (c,_) -> c1 = c) lcp)))
- lp
-let critere3 ((i,j),m) lp lcp =
+let critere3 table ((i,j),m) lp lcp =
List.exists
(fun h ->
h.num <> i && h.num <> j
&& (div_mon_test m (lm h))
&& (h.num < j
|| not (m = ppcm_mon
- (lm (!allpol.(i)))
+ (lm (table.allpol.(i)))
(lm h)))
&& (h.num < i
|| not (m = ppcm_mon
- (lm (!allpol.(j)))
+ (lm (table.allpol.(j)))
(lm h))))
lp
let add_cpairs p lp lcp =
mergecpairs (cpairs1 p lp) lcp
-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))
- ^ "]"))
+ (info (fun () -> Printf.sprintf "[%i,%i]" (List.length p) (List.length q)))
(* in lp new polynomials are at the end *)
-let coef_courant = ref coef1
-
type certificate =
{ coef : coef; power : int;
gb_comb : poly list list; last_comb : poly list }
-let test_dans_ideal p lp lp0 =
- let (c,r) = reduce2 (ppol !pol_courant) lp in
- info ("remainder: "^(stringPcut r)^"\n");
- coef_courant:= P.multP !coef_courant c;
- pol_courant:= mk_polynom r;
- if r=[]
- then (info "polynomial reduced to 0\n";
- let lcp = List.map (fun q -> []) !poldep in
- let c = !coef_courant in
- let (lcq,r) = reduce2_trace (emultP c p) lp lcp in
- info "r ok\n";
- info ("r: "^(stringP r)^"\n");
- let res=ref (emultP c p) in
- List.iter2
- (fun cq q -> res:=plusP (!res) (multP cq (ppol q));
- )
- lcq !poldep;
- info ("verif sum: "^(stringP (!res))^"\n");
- info ("coefficient: "^(stringP (polconst 1 c))^"\n");
- let rec aux lp =
- match lp with
- |[] -> []
- |p::lp ->
- (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 (* (map 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
- (liste_polynomes_de_depart,
- polynome_a_tester,
- {coef = coefficient_multiplicateur;
- power = 1;
- gb_comb = liste_des_coefficients_intermediaires;
- last_comb = liste_des_coefficients})
- )
- else ((*info "polynomial not reduced to 0\n";
- info ("\nremainder: "^(stringPcut r)^"\n");*)
- raise NotInIdeal)
-
-let divide_rem_with_critical_pair = ref false
-
-let list_diff l x =
- List.filter (fun y -> y <> x) l
+type current_problem = {
+ cur_poly : poly;
+ cur_coef : coef;
+}
+
+exception NotInIdealUpdate of current_problem
+
+let test_dans_ideal cur_pb table metadata p lp len0 =
+ (** Invariant: [lp] is [List.tl (Array.to_list table.allpol)] *)
+ let (c,r) = reduce2 table cur_pb.cur_poly lp in
+ info (fun () -> "remainder: "^(stringPcut metadata r));
+ let cur_pb = {
+ cur_coef = P.multP cur_pb.cur_coef c;
+ cur_poly = r;
+ } in
+ match r with
+ | [] ->
+ sinfo "polynomial reduced to 0";
+ let lcp = List.map (fun q -> []) lp in
+ let c = cur_pb.cur_coef in
+ let (lcq,r) = reduce2_trace table (emultP c p) lp lcp in
+ sinfo "r ok";
+ info (fun () -> "r: "^(stringP metadata r));
+ info (fun () ->
+ let fold res cq q = plusP res (multP cq (ppol q)) in
+ let res = List.fold_left2 fold (emultP c p) lcq lp in
+ "verif sum: "^(stringP metadata res)
+ );
+ info (fun () -> "coefficient: "^(stringP metadata (polconst 1 c)));
+ let coefficient_multiplicateur = c in
+ let liste_des_coefficients_intermediaires =
+ let rec aux accu lp =
+ match lp with
+ | [] -> accu
+ | p :: lp ->
+ let elt = List.map (fun q -> coefpoldep_find table p q) lp in
+ aux (elt :: accu) lp
+ in
+ let lci = aux [] (List.rev lp) in
+ CList.skipn len0 lci
+ in
+ let liste_des_coefficients =
+ List.rev_map (fun cq -> emultP (coef_of_int (-1)) cq) lcq
+ in
+ {coef = coefficient_multiplicateur;
+ power = 1;
+ gb_comb = liste_des_coefficients_intermediaires;
+ last_comb = liste_des_coefficients}
+ | _ -> raise (NotInIdealUpdate cur_pb)
let deg_hom p =
match p with
| [] -> -1
- | (a,m)::_ -> m.(0)
-
-let pbuchf pq p lp0=
- info "computation of the Groebner basis\n";
- step:=0;
- Hashtbl.clear hmon;
- let rec pbuchf (lp, lpc) =
+ | (a,m)::_ -> Monomial.deg m
+
+let pbuchf table metadata cur_pb homogeneous (lp, lpc) p =
+ (** Invariant: [lp] is [List.tl (Array.to_list table.allpol)] *)
+ sinfo "computation of the Groebner basis";
+ let () = match table.hmon with
+ | None -> ()
+ | Some hmon -> Hashtbl.clear hmon
+ in
+ let len0 = List.length lp in
+ let rec pbuchf cur_pb (lp, lpc) =
infobuch lp lpc;
-(* step:=(!step+1)mod 10;*)
match lpc with
[] ->
-
- (* info ("List of polynomials:\n"^(fold_left (fun r p -> r^(stringP p)^"\n") "" lp));
- info "--------------------\n";*)
- test_dans_ideal (ppol p) lp lp0
+ test_dans_ideal cur_pb table metadata p lp len0
| ((i,j),m) :: lpc2 ->
-(* info "choosen pair\n";*)
- if critere3 ((i,j),m) lp lpc2
- then (info "c"; pbuchf (lp, lpc2))
+ if critere3 table ((i,j),m) lp lpc2
+ then (sinfo "c"; pbuchf cur_pb (lp, lpc2))
else
- let a = spol !allpol.(i) !allpol.(j) in
- if !homogeneous && (ppol a)<>[] && deg_hom (ppol a)
- > deg_hom (ppol !pol_courant)
- then (info "h"; pbuchf (lp, lpc2))
+ let (a0, p0, q0) = spol0 table.allpol.(i) table.allpol.(j) in
+ if homogeneous && a0 <>[] && deg_hom a0 > deg_hom cur_pb.cur_poly
+ then (sinfo "h"; pbuchf cur_pb (lp, lpc2))
else
(* let sa = a.sugar in*)
- let (ca,a0)= reduce2 (ppol a) lp in
- match a0 with
- [] -> info "0";pbuchf (lp, lpc2)
- | _ ->
+ match reduce2 table a0 lp with
+ _, [] -> sinfo "0";pbuchf cur_pb (lp, lpc2)
+ | ca, _ :: _ ->
(* info "pair reduced\n";*)
- a.pol := emultP ca (ppol a);
- let (lca,a0) = reduce2_trace (ppol a) lp
- (List.map (fun q -> emultP ca (coefpoldep_find a q))
- !poldep) in
+ let map q =
+ let r =
+ if q.num == i then p0 else if q.num == j then q0 else []
+ in
+ emultP ca r
+ in
+ let lcp = List.map map lp in
+ let (lca, a0) = reduce2_trace table (emultP ca a0) lp lcp in
(* info "paire re-reduced";*)
- a.pol := a0;
-(* let a0 = new_allpol a0 sa in*)
- List.iter2 (fun c q ->
- coefpoldep_remove a q;
- coefpoldep_set a q c) lca !poldep;
+ let a = new_allpol table a0 in
+ List.iter2 (fun c q -> coefpoldep_set table a q c) lca lp;
let a0 = a in
- info ("\nnew polynomial: "^(stringPcut (ppol a0))^"\n");
- let ct = coef1 (* contentP a0 *) in
- (*info ("content: "^(string_of_coef ct)^"\n");*)
- poldep:=addS a0 lp;
- poldepcontent:=addS ct (!poldepcontent);
-
- try test_dans_ideal (ppol p) (addS a0 lp) lp0
- with NotInIdeal ->
+ info (fun () -> "new polynomial: "^(stringPcut metadata (ppol a0)));
+ let nlp = addS a0 lp in
+ try test_dans_ideal cur_pb table metadata p nlp len0
+ with NotInIdealUpdate cur_pb ->
let newlpc = add_cpairs a0 lp lpc2 in
- pbuchf (((addS a0 lp), newlpc))
- in pbuchf pq
+ pbuchf cur_pb (nlp, newlpc)
+ in pbuchf cur_pb (lp, lpc)
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
+ | (a,m)::p1 -> let d = deg m in
+ List.for_all (fun (b,m') -> deg m' =d) p1
(* returns
c
@@ -955,33 +781,33 @@ let is_homogeneous p =
where pn+k = a(n+k,n+k-1)*pn+k-1 + ... + a(n+k,1)* p1
*)
-let in_ideal d lp p =
- Hashtbl.clear hmon;
- Hashtbl.clear coefpoldep;
- nallpol := 0;
- allpol := Array.make 1000 polynom0;
- homogeneous := List.for_all is_homogeneous (p::lp);
- if !homogeneous then info "homogeneous polynomials\n";
- info ("p: "^(stringPcut p)^"\n");
- info ("lp:\n"^(List.fold_left (fun r p -> r^(stringPcut p)^"\n") "" lp));
- (*info ("p: "^(stringP p)^"\n");
- info ("lp:\n"^(fold_left (fun r p -> r^(stringP p)^"\n") "" lp));*)
-
- let lp = List.map mk_polynom lp in
- let p = mk_polynom p in
- initcoefpoldep d lp;
- coef_courant:=coef1;
- pol_courant:=p;
-
- let (lp1,p1,cert) =
- try test_dans_ideal (ppol p) lp lp
- with NotInIdeal -> pbuchf (lp, (cpairs lp)) p lp in
- info "computed\n";
-
- (List.map ppol lp1, p1, cert)
-
-(* *)
-end
-
-
+let in_ideal metadata d lp p =
+ let table = {
+ hmon = None;
+ coefpoldep = Hashtbl.create 51;
+ nallpol = 0;
+ allpol = Array.make 1000 polynom0;
+ } in
+ let homogeneous = List.for_all is_homogeneous (p::lp) in
+ if homogeneous then sinfo "homogeneous polynomials";
+ info (fun () -> "p: "^(stringPcut metadata p));
+ info (fun () -> "lp:\n"^(List.fold_left (fun r p -> r^(stringPcut metadata p)^"\n") "" lp));
+
+ let lp = List.map (fun c -> new_allpol table c) lp in
+ List.iter (fun p -> coefpoldep_set table p p (polconst d (coef_of_int 1))) lp;
+ let cur_pb = {
+ cur_poly = p;
+ cur_coef = coef1;
+ } in
+
+ let cert =
+ try pbuchf table metadata cur_pb homogeneous (lp, []) p
+ with NotInIdealUpdate cur_pb ->
+ try pbuchf table metadata cur_pb homogeneous (lp, (cpairs lp)) p
+ with NotInIdealUpdate _ -> raise NotInIdeal
+ in
+ sinfo "computed";
+
+ cert
+end
diff --git a/plugins/nsatz/ideal.mli b/plugins/nsatz/ideal.mli
index d1a2a0a7d..b7ec901af 100644
--- a/plugins/nsatz/ideal.mli
+++ b/plugins/nsatz/ideal.mli
@@ -6,6 +6,17 @@
(* * GNU Lesser General Public License Version 2.1 *)
(************************************************************************)
+type metadata = {
+ name_var : string list;
+}
+
+module Monomial :
+sig
+type t
+val repr : t -> int array
+val make : int array -> t
+end
+
module Make (P : Polynom.S) :
sig
(* Polynomials *)
@@ -14,32 +25,26 @@ type deg = int
type coef = P.t
type poly
-val repr : poly -> (coef * int array) list
+val repr : poly -> (coef * Monomial.t) list
val polconst : int -> coef -> poly
val zeroP : poly
val gen : int -> int -> poly
val equal : poly -> poly -> bool
-val name_var : string list ref
val plusP : poly -> poly -> poly
val oppP : poly -> poly
val multP : poly -> poly -> poly
val puisP : poly -> int -> poly
-val poldepcontent : coef list ref
-
type certificate =
{ coef : coef; power : int;
gb_comb : poly list list; last_comb : poly list }
-val in_ideal : deg -> poly list -> poly -> poly list * poly * certificate
+val in_ideal : metadata -> deg -> poly list -> poly -> certificate
module Hashpol : Hashtbl.S with type key = poly
-val sugar_flag : bool ref
-val divide_rem_with_critical_pair : bool ref
-
end
exception NotInIdeal
diff --git a/plugins/nsatz/nsatz.ml b/plugins/nsatz/nsatz.ml
index 36bce780b..8fba92e79 100644
--- a/plugins/nsatz/nsatz.ml
+++ b/plugins/nsatz/nsatz.ml
@@ -257,21 +257,19 @@ let rec parse_request lp =
(parse_term p)::(parse_request lp1)
|_-> assert false
-let nvars = ref 0
-
-let set_nvars_term t =
- let rec aux t =
+let set_nvars_term nvars t =
+ let rec aux t nvars =
match t with
- | Zero -> ()
- | Const r -> ()
+ | Zero -> nvars
+ | Const r -> nvars
| 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
+ max nvars n
+ | Opp t1 -> aux t1 nvars
+ | Add (t1,t2) -> aux t2 (aux t1 nvars)
+ | Sub (t1,t2) -> aux t2 (aux t1 nvars)
+ | Mul (t1,t2) -> aux t2 (aux t1 nvars)
+ | Pow (t1,n) -> aux t1 nvars
+ in aux t nvars
let string_of_term p =
let rec aux p =
@@ -301,8 +299,8 @@ open PIdeal
varaibles <=np are in the coefficients
*)
-let term_pol_sparse np t=
- let d = !nvars in
+let term_pol_sparse nvars np t=
+ let d = nvars in
let rec aux t =
(* info ("conversion de: "^(string_of_term t)^"\n");*)
let res =
@@ -336,14 +334,8 @@ let polrec_to_term p =
match p with
|Poly.Pint n -> const (Coef.to_num n)
|Poly.Prec (v,coefs) ->
- let res = ref Zero in
- Array.iteri
- (fun i c ->
- res:=add(!res, mul(aux c,
- pow (Var (string_of_int v),
- i))))
- coefs;
- !res
+ let fold i c res = add (res, mul (aux c, pow (Var (string_of_int v), i))) in
+ Array.fold_right_i fold coefs Zero
in aux p
(* approximation of the Horner form used in the tactic ring *)
@@ -355,9 +347,11 @@ let pol_sparse_to_term n2 p =
match p with
[] -> const (num_of_string "0")
| (a,m)::p1 ->
+ let m = Ideal.Monomial.repr m in
let n = (Array.length m)-1 in
let (i0,e0) =
List.fold_left (fun (r,d) (a,m) ->
+ let m = Ideal.Monomial.repr m in
let i0= ref 0 in
for k=1 to n do
if m.(k)>0
@@ -374,45 +368,28 @@ let pol_sparse_to_term n2 p =
p in
if Int.equal i0 0
then
- let mp = ref (polrec_to_term a) in
- if List.is_empty 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 mp = polrec_to_term a in
+ if List.is_empty p1 then mp else add (mp, aux p1)
+ else
+ let fold (p1, p2) (a, m) =
+ if (Ideal.Monomial.repr m).(i0) >= e0 then begin
+ let m0 = Array.copy (Ideal.Monomial.repr m) in
+ let () = m0.(i0) <- m0.(i0) - e0 in
+ let m0 = Ideal.Monomial.make m0 in
+ ((a, m0) :: p1, p2)
+ end else
+ (p1, (a, m) :: p2)
+ in
+ let (p1, p2) = List.fold_left fold ([], []) p in
let vm =
if Int.equal 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))))
+ add (mul(vm, aux (List.rev p1)), aux (List.rev p2))
in (*info "-> pol_sparse_to_term\n";*)
aux p
-let remove_list_tail l i =
- let rec aux l i =
- if List.is_empty l
- then []
- else if i<0
- then l
- else if Int.equal 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]
@@ -422,49 +399,35 @@ let remove_list_tail l i =
removes intermediate polynomials not useful to compute the last one.
*)
-let remove_zeros zero lci =
- let n = List.length (List.hd lci) in
- let m=List.length lci in
+let remove_zeros lci =
+ let m = List.length lci in
let u = Array.make m false in
let rec utiles k =
- if k>=m
- then ()
- else (
- u.(k)<-true;
+ (** TODO: Find a more reasonable implementation of this traversal. *)
+ if k >= m || u.(k) then ()
+ else
+ let () = u.(k) <- true in
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 ("useless spolynomials: "
- ^string_of_int (m-List.length lr)^"\n");
- info ("useful spolynomials: "
- ^string_of_int (List.length lr)^"\n");
+ let iter i c = if not (PIdeal.equal c zeroP) then utiles (i + k + 1) in
+ List.iteri iter lc
+ in
+ let () = utiles 0 in
+ let filter i l =
+ let f j l = if m <= i + j + 1 then true else u.(i + j + 1) in
+ if u.(i) then Some (List.filteri f l)
+ else None
+ in
+ let lr = CList.map_filter_i filter lci in
+ info (fun () -> Printf.sprintf "useless spolynomials: %i" (m-List.length lr));
+ info (fun () -> Printf.sprintf "useful spolynomials: %i " (List.length lr));
lr
-let theoremedeszeros lpol p =
+let theoremedeszeros metadata nvars lpol p =
let t1 = Unix.gettimeofday() in
- let m = !nvars in
- let (lp0,p,cert) = in_ideal m lpol p in
- let lpc = List.rev !poldepcontent in
- info ("time: "^Format.sprintf "@[%10.3f@]s\n" (Unix.gettimeofday ()-.t1));
- (cert,lp0,p,lpc)
+ let m = nvars in
+ let cert = in_ideal metadata m lpol p in
+ info (fun () -> Printf.sprintf "time: @[%10.3f@]s" (Unix.gettimeofday ()-.t1));
+ cert
open Ideal
@@ -507,51 +470,33 @@ let expand_pol lb lp =
in List.rev (aux lb (List.rev lp))
let theoremedeszeros_termes lp =
- nvars:=0;(* mise a jour par term_pol_sparse *)
- List.iter set_nvars_term lp;
+ let nvars = List.fold_left set_nvars_term 0 lp in
match lp with
| Const (Int sugarparam)::Const (Int nparam)::lp ->
((match sugarparam with
- |0 -> info "computation without sugar\n";
+ |0 -> sinfo "computation without sugar";
lexico:=false;
- sugar_flag := false;
- divide_rem_with_critical_pair := false
- |1 -> info "computation with sugar\n";
+ |1 -> sinfo "computation with sugar";
lexico:=false;
- sugar_flag := true;
- divide_rem_with_critical_pair := false
- |2 -> info "ordre lexico computation without sugar\n";
+ |2 -> sinfo "ordre lexico computation without sugar";
lexico:=true;
- sugar_flag := false;
- divide_rem_with_critical_pair := false
- |3 -> info "ordre lexico computation with sugar\n";
+ |3 -> sinfo "ordre lexico computation with sugar";
lexico:=true;
- sugar_flag := true;
- divide_rem_with_critical_pair := false
- |4 -> info "computation without sugar, division by pairs\n";
+ |4 -> sinfo "computation without sugar, division by pairs";
lexico:=false;
- sugar_flag := false;
- divide_rem_with_critical_pair := true
- |5 -> info "computation with sugar, division by pairs\n";
+ |5 -> sinfo "computation with sugar, division by pairs";
lexico:=false;
- sugar_flag := true;
- divide_rem_with_critical_pair := true
- |6 -> info "ordre lexico computation without sugar, division by pairs\n";
+ |6 -> sinfo "ordre lexico computation without sugar, division by pairs";
lexico:=true;
- sugar_flag := false;
- divide_rem_with_critical_pair := true
- |7 -> info "ordre lexico computation with sugar, division by pairs\n";
+ |7 -> sinfo "ordre lexico computation with sugar, division by pairs";
lexico:=true;
- sugar_flag := true;
- divide_rem_with_critical_pair := true
| _ -> error "nsatz: bad parameter"
);
- let m= !nvars in
- let lvar=ref [] in
- for i=m downto 1 do lvar:=["x"^(string_of_int i)^""]@(!lvar); done;
- lvar:=["a";"b";"c";"d";"e";"f";"g";"h";"i";"j";"k";"l";"m";"n";"o";"p";"q";"r";"s";"t";"u";"v";"w";"x";"y";"z"] @ (!lvar); (* pour macaulay *)
- name_var:=!lvar;
- let lp = List.map (term_pol_sparse nparam) lp in
+ let lvar = List.init nvars (fun i -> Printf.sprintf "x%i" (i + 1)) in
+ let lvar = ["a";"b";"c";"d";"e";"f";"g";"h";"i";"j";"k";"l";"m";"n";"o";"p";"q";"r";"s";"t";"u";"v";"w";"x";"y";"z"] @ lvar in
+ (* pour macaulay *)
+ let metadata = { name_var = lvar } in
+ let lp = List.map (term_pol_sparse nvars nparam) lp in
match lp with
| [] -> assert false
| p::lp1 ->
@@ -561,16 +506,16 @@ let theoremedeszeros_termes lp =
lb is kept in order to fix the certificate in the post-processing
*)
let lpol, lb = clean_pol lpol in
- let (cert,lp0,p,_lct) = theoremedeszeros lpol p in
- info "cert ok\n";
+ let cert = theoremedeszeros metadata nvars lpol p in
+ sinfo "cert ok";
let lc = cert.last_comb::List.rev cert.gb_comb in
- match remove_zeros (fun x -> equal x zeroP) lc with
+ match remove_zeros lc with
| [] -> assert false
| (lq::lci) ->
(* post-processing : we apply the correction for the last line *)
let lq = expand_pol lb lq in
(* lci commence par les nouveaux polynomes *)
- let m = !nvars in
+ let m = nvars in
let c = pol_sparse_to_term m (polconst m cert.coef) in
let r = Pow(Zero,cert.power) in
let lci = List.rev lci in
@@ -578,8 +523,8 @@ let theoremedeszeros_termes lp =
let lci = List.map (expand_pol lb) lci in
let lci = List.map (List.map (pol_sparse_to_term m)) lci in
let lq = List.map (pol_sparse_to_term m) lq in
- info ("number of parameters: "^string_of_int nparam^"\n");
- info "term computed\n";
+ info (fun () -> Printf.sprintf "number of parameters: %i" nparam);
+ sinfo "term computed";
(c,r,lci,lq)
)
|_ -> assert false
@@ -619,7 +564,7 @@ let nsatz lpol =
mkt_app lcons [tlp ();ltterm;r])
res
(mkt_app lnil [tlp ()]) in
- info "term computed\n";
+ sinfo "term computed";
res
let return_term t =
@@ -633,5 +578,3 @@ let nsatz_compute t =
with Ideal.NotInIdeal ->
error "nsatz cannot solve this problem" in
return_term lpol
-
-
diff --git a/plugins/nsatz/utile.ml b/plugins/nsatz/utile.ml
index 922432460..d3cfd75e5 100644
--- a/plugins/nsatz/utile.ml
+++ b/plugins/nsatz/utile.ml
@@ -11,8 +11,8 @@ let prt0 s = () (* print_string s;flush(stdout)*)
let prt s =
if !Flags.debug then (print_string (s^"\n");flush(stdout)) else ()
-let info s =
- Flags.if_verbose prerr_string s
+let sinfo s = if !Flags.debug then Feedback.msg_debug (Pp.str s)
+let info s = if !Flags.debug then Feedback.msg_debug (Pp.str (s ()))
(* Lists *)
diff --git a/plugins/nsatz/utile.mli b/plugins/nsatz/utile.mli
index 1f8415752..9308577e0 100644
--- a/plugins/nsatz/utile.mli
+++ b/plugins/nsatz/utile.mli
@@ -4,7 +4,8 @@ val pr : string -> unit
val prn : string -> unit
val prt0 : 'a -> unit
val prt : string -> unit
-val info : string -> unit
+val info : (unit -> string) -> unit
+val sinfo : string -> unit
(* Listes *)
val list_mem_eq : ('a -> 'b -> bool) -> 'a -> 'b list -> bool