aboutsummaryrefslogtreecommitdiffhomepage
path: root/plugins/groebner/polynom.ml
diff options
context:
space:
mode:
Diffstat (limited to 'plugins/groebner/polynom.ml')
-rw-r--r--plugins/groebner/polynom.ml1056
1 files changed, 1056 insertions, 0 deletions
diff --git a/plugins/groebner/polynom.ml b/plugins/groebner/polynom.ml
new file mode 100644
index 000000000..5c859f755
--- /dev/null
+++ b/plugins/groebner/polynom.ml
@@ -0,0 +1,1056 @@
+(************************************************************************)
+
+(* v * The Coq Proof Assistant / The Coq Development Team *)
+(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *)
+(* \VV/ **************************************************************)
+(* // * This file is distributed under the terms of the *)
+(* * GNU Lesser General Public License Version 2.1 *)
+(************************************************************************)
+
+(*
+ Polynomes récursifs: Z[x1]...[xn].
+*)
+open Utile
+open Util
+
+(***********************************************************************
+ 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.filter (fun q -> deg_total q >0)
+ (List.flatten (factorise (normc p))))
+
+
+(* produit des facteurs non constants d'une décomposition sans carré *)
+let pol_de_factorisation lf =
+ let p = ref cf1 in
+ List.iter (fun lq ->
+ Array.iteri (fun i q -> if (deg_total q)>0 then p:=(!p)@@(q^^(i+1)))
+ (Array.of_list lq))
+ lf;
+ !p
+
+
+let set_of_array_facteurs tf =
+ let h = Hashpol.create 51 in
+ Array.iter (fun lf ->
+ List.iter (fun p -> if not (Hashpol.mem h p)
+ then Hashpol.add h p true)
+ lf)
+ tf;
+ let res = ref [] in
+ Hashpol.iter (fun p _ -> res:=p::(!res)) h;
+ !res
+
+
+(* Factorise un tableau de polynômes f, et rend:
+ - un tableau p de facteurs (degré>0, contenu entier 1,
+ coefficient de tête >0) obtenu par décomposition sans carrés
+ puis par division mutuelle
+ - un tableau l de couples (constante, listes d'indices l)
+ tels que f.(i) = l.(i)_1*Produit(p.(j), j dans l.(i)_2)
+*)
+
+(* on donne le tableau des facteurs de chaque poly de f *)
+let factorise_tableauP2 f l1 =
+ let x = max_var f in
+ (* liste des facteurs sans carré des polynômes de f *)
+ pr"<";
+ let l1 = set_of_array_facteurs l1 in
+ (* on les divise entre eux pour éventuellement trouver
+ de nouveaux facteurs *)
+ pr "~";
+ let l1 = Sort.list (fun p q -> (deg_total p)<(deg_total q)) l1 in
+ let l1 = Array.of_list (facteurs_liste (fun a b -> div_pol a b x)
+ (fun p -> (deg_total p)<1)
+ l1) in
+ (* puis on décompose les polynômes de f avec ces facteurs *)
+ pr "-";
+ let res = factorise_tableau (fun a b -> div_pol a b x)
+ (fun p -> 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_for_all 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
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-10 05:59:33 +0000