diff options
Diffstat (limited to 'plugins/nsatz/ideal.ml')
-rw-r--r-- | plugins/nsatz/ideal.ml | 818 |
1 files changed, 293 insertions, 525 deletions
diff --git a/plugins/nsatz/ideal.ml b/plugins/nsatz/ideal.ml index 48bdad82..f8fc9437 100644 --- a/plugins/nsatz/ideal.ml +++ b/plugins/nsatz/ideal.ml @@ -1,9 +1,11 @@ (************************************************************************) -(* v * The Coq Proof Assistant / The Coq Development Team *) -(* <O___,, * INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2016 *) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* <O___,, * (see CREDITS file for the list of authors) *) (* \VV/ **************************************************************) -(* // * This file is distributed under the terms of the *) -(* * GNU Lesser General Public License Version 2.1 *) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(* * (see LICENSE file for the text of the license) *) (************************************************************************) (* Nullstellensatz with Groebner basis computation @@ -23,7 +25,6 @@ exception NotInIdeal Global options *) let lexico = ref false -let use_hmon = ref false (* division of tail monomials *) @@ -33,31 +34,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 +104,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 +132,45 @@ 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 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 @@ -163,8 +198,6 @@ module Hashpol = Hashtbl.Make( (* A pretty printer for polynomials, with Maple-like syntax. *) -open Format - let getvar lv i = try (List.nth lv i) with Failure _ -> (List.fold_left (fun r x -> r^" "^x) "lv= " lv) @@ -179,8 +212,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 @@ -218,62 +251,7 @@ let string_of_pol zeroP hdP tlP coefterm monterm string_of_coef in (stringP p true) - - -let print_pol zeroP hdP tlP coefterm monterm string_of_coef - dimmon string_of_exp lvar p = - - let rec print_mon m coefone = - let s=ref [] in - for i=1 to (dimmon m) do - (match (string_of_exp m i) with - "0" -> () - | "1" -> s:= (!s) @ [(getvar !lvar (i-1))] - | e -> s:= (!s) @ [((getvar !lvar (i-1)) ^ "^" ^ e)]); - done; - (match !s with - [] -> if coefone - then print_string "1" - else () - | l -> if coefone - then print_string (String.concat "*" l) - else (print_string "*"; - print_string (String.concat "*" l))) - and print_term t start = let a = coefterm t and m = monterm t in - match (string_of_coef a) with - "0" -> () - | "1" ->(match start with - true -> print_mon m true - |false -> (print_string "+ "; - print_mon m true)) - | "-1" ->(print_string "-";print_space();print_mon m true) - | c -> if (String.get c 0)='-' - then (print_string "- "; - print_string (String.sub c 1 - ((String.length c)-1)); - print_mon m false) - else (match start with - true -> (print_string c;print_mon m false) - |false -> (print_string "+ "; - print_string c;print_mon m false)) - and printP p start = - if (zeroP p) - then (if start - then print_string("0") - else ()) - else (print_term (hdP p) start; - if start then open_hovbox 0; - print_space(); - print_cut(); - printP (tlP p) false) - in open_hovbox 3; - printP p true; - print_flush() - - -let name_var= ref [] - -let stringP p = +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 +259,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 +309,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 +340,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 +350,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 +391,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 +442,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 +462,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 +477,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 +490,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 +519,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,301 +530,183 @@ 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 + +module CPair = +struct +type t = (int * int) * Monomial.t +let compare ((i1, j1), m1) ((i2, j2), m2) = compare_mon m2 m1 +end + +module Heap : +sig + type elt = (int * int) * Monomial.t + type t + val length : t -> int + val empty : t + val add : elt -> t -> t + val pop : t -> (elt * t) option +end = +struct + include Heap.Functional(CPair) + let length h = fold (fun _ accu -> accu + 1) h 0 + let pop h = try Some (maximum h, remove h) with Heap.EmptyHeap -> None +end let ord i j = if i<j then (i,j) else (j,i) -let cpair p q = - if etrangers (ppol p) (ppol q) - then [] - else [(ord p.num q.num, - ppcm_mon (lm p) (lm q))] - -let cpairs1 p lq = - sortcpairs (List.fold_left (fun r q -> r @ (cpair p q)) [] lq) - -let cpairs lp = - let rec aux l = - match l with - []|[_] -> [] - |p::l1 -> mergecpairs (cpairs1 p l1) (aux l1) - in aux lp - - -let critere2 ((i,j),m) lp lcp = - 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 cpair p q accu = + if etrangers (ppol p) (ppol q) then accu + else Heap.add (ord p.num q.num, ppcm_mon (lm p) (lm q)) accu + +let cpairs1 p lq accu = + List.fold_left (fun r q -> cpair p q r) accu lq -let critere3 ((i,j),m) lp lcp = +let rec cpairs l accu = match l with +| [] | [_] -> accu +| p :: l -> + cpairs l (cpairs1 p l accu) + +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) (Heap.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 - | ((i,j),m) :: lpc2 -> -(* info "choosen pair\n";*) - if critere3 ((i,j),m) lp lpc2 - then (info "c"; pbuchf (lp, lpc2)) + match Heap.pop lpc with + | None -> + test_dans_ideal cur_pb table metadata p lp len0 + | Some (((i, j), m), 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 -> - let newlpc = add_cpairs a0 lp lpc2 in - pbuchf (((addS a0 lp), newlpc)) - in pbuchf pq + 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 = cpairs1 a0 lp lpc2 in + 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 +723,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, Heap.empty) p + with NotInIdealUpdate cur_pb -> + try pbuchf table metadata cur_pb homogeneous (lp, cpairs lp Heap.empty) p + with NotInIdealUpdate _ -> raise NotInIdeal + in + sinfo "computed"; + + cert +end |