1 files changed, 679 insertions, 0 deletions

diff --git a/plugins/nsatz/polynom.ml b/plugins/nsatz/polynom.ml
new file mode 100644
index 00000000..14e279b5
--- /dev/null
+++ b/plugins/nsatz/polynom.ml
@@ -0,0 +1,679 @@
+(************************************************************************)
+(*  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        *)
+(************************************************************************)
+
+(* Recursive polynomials: R[x1]...[xn]. *)
+open Utile
+open Util
+
+(* 1. Coefficients: R *)
+
+module type Coef = sig
+  type t
+  val equal : t -> t -> bool
+  val lt : t -> t -> bool
+  val le : t -> t -> bool
+  val abs : t -> t
+  val plus : t -> t -> t
+  val mult : t -> t -> t
+  val sub : t -> t -> t
+  val opp : t -> t
+  val div : t -> t -> t
+  val modulo : t -> t -> t
+  val puis : t -> int -> t
+  val pgcd : t -> t -> t
+
+  val hash : t -> int
+  val of_num : Num.num -> t
+  val to_string : t -> string
+end
+
+module type S = sig
+  type coef
+  type variable = int
+  type t = Pint of coef | Prec of variable * t array
+
+  val of_num : Num.num -> t
+  val x : variable -> t
+  val monome : variable -> int -> t
+  val is_constantP : t -> bool
+  val is_zero : t -> bool
+
+  val max_var_pol : t -> variable
+  val max_var_pol2 : t -> variable
+  val max_var : t array -> variable
+  val equal : t -> t -> bool
+  val norm : t -> t
+  val deg : variable -> t -> int
+  val deg_total : t -> int
+  val copyP : t -> t
+  val coef : variable -> int -> t -> t
+
+  val plusP : t -> t -> t
+  val content : t -> coef
+  val div_int : t -> coef -> t
+  val vire_contenu : t -> t
+  val vars : t -> variable list
+  val int_of_Pint : t -> coef
+  val multx : int -> variable -> t -> t
+  val multP : t -> t -> t
+  val deriv : variable -> t -> t
+  val oppP : t -> t
+  val moinsP : t -> t -> t
+  val puisP : t -> int -> t
+  val ( @@ ) : t -> t -> t
+  val ( -- ) : t -> t -> t
+  val ( ^^ ) : t -> int -> t
+  val coefDom : variable -> t -> t
+  val coefConst : variable -> t -> t
+  val remP : variable -> t -> t
+  val coef_int_tete : t -> coef
+  val normc : t -> t
+  val coef_constant : t -> coef
+  val univ : bool ref
+  val string_of_var : int -> string
+  val nsP : int ref
+  val to_string : t -> string
+  val printP : t -> unit
+  val print_tpoly : t array -> unit
+  val print_lpoly : t list -> unit
+  val quo_rem_pol : t -> t -> variable -> t * t
+  val div_pol : t -> t -> variable -> t
+  val divP : t -> t -> t
+  val div_pol_rat : t -> t -> bool
+  val pseudo_div : t -> t -> variable -> t * t * int * t
+  val pgcdP : t -> t -> t
+  val pgcd_pol : t -> t -> variable -> t
+  val content_pol : t -> variable -> t
+  val pgcd_coef_pol : t -> t -> variable -> t
+  val pgcd_pol_rec : t -> t -> variable -> t
+  val gcd_sub_res : t -> t -> variable -> t
+  val gcd_sub_res_rec : t -> t -> t -> t -> int -> variable -> t
+  val lazard_power : t -> t -> int -> variable -> t
+  val hash : t -> int
+  module Hashpol : Hashtbl.S with type key=t
+end
+
+(***********************************************************************
+  2. Type of polynomials, operations.
+*)
+module Make (C:Coef) = struct
+
+type coef = C.t
+let coef_of_int i = C.of_num (Num.Int i)
+let coef0 = coef_of_int 0
+let coef1 = coef_of_int 1
+
+type variable = int
+
+type t =
+    Pint of coef                    (* constant polynomial *)
+  | Prec of variable * (t array)    (* coefficients, increasing degree *)
+
+(* by default, operations work with normalized polynomials:
+- variables are positive integers
+- coefficients of a polynomial in x only use variables < x
+- no zero coefficient at beginning
+- no Prec(x,a) where a is constant in x
+*)
+
+(* constant polynomials *)
+let of_num x = Pint (C.of_num x)
+let cf0 = of_num (Num.Int 0)
+let cf1 = of_num (Num.Int 1)
+
+(* nth variable *)
+let x n = Prec (n,[|cf0;cf1|])
+
+(* create 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)
+
+let is_constantP = function
+    Pint _ -> true
+  | Prec _ -> false
+
+let int_of_Pint = function
+    Pint x -> x
+  | _ -> failwith "non"
+
+let is_zero p =
+  match p with Pint n -> if C.equal n coef0 then true else false |_-> false
+
+let max_var_pol p =
+  match p with
+      Pint _ -> 0
+    |Prec(x,_) -> x
+
+(* p not normalized *)
+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
+
+let rec max_var l = Array.fold_right (fun p m -> max (max_var_pol2 p) m) l 0
+
+(* equality between polynomials *)
+
+let rec equal p q =
+  match (p,q) with
+      (Pint a,Pint b) -> C.equal a b
+    |(Prec(x,p1),Prec(y,q1)) ->
+       if x<>y then false
+       else if (Array.length p1)<>(Array.length q1) then false
+       else (try (Array.iteri (fun i a -> if not (equal a q1.(i))
+			       then failwith "raté")
+		    p1;
+		  true)
+	     with _ -> false)
+    | (_,_) -> false
+
+(* normalize polynomial: remove head zeros, coefficients are normalized
+   if constant, returns the coefficient
+*)
+
+let rec norm p = match p with
+    Pint _ -> p
+  |Prec (x,a)->
+     let d = (Array.length a -1) in
+     let n = ref d in
+       while !n>0 && (equal a.(!n) (Pint coef0)) do
+	 n:=!n-1;
+       done;
+       if !n<0 then Pint coef0
+       else if !n=0 then a.(0)
+       else if !n=d then p
+       else (let b=Array.create (!n+1) (Pint coef0) in
+               for i=0 to !n do b.(i)<-a.(i);done;
+               Prec(x,b))
+
+
+(* degree in v, v >= max var of p *)
+let rec deg v p =
+  match p with
+      Prec(x,p1) when x=v -> Array.length p1 -1
+    |_ -> 0
+
+
+(* total degree *)
+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
+
+let rec copyP p =
+  match p with
+      Pint i -> Pint i
+    |Prec(x,q) -> Prec(x,Array.map copyP q)
+
+(* coefficient of degree i in v, v >= max var of 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
+
+(* addition *)
+
+let rec plusP p q =
+  let res =
+    (match (p,q) with
+	 (Pint a,Pint b) -> Pint (C.plus a b)
+       |(Pint a, Prec (y,q1)) -> let q2=Array.map copyP q1 in
+           q2.(0)<- plusP p q1.(0);
+           Prec (y,q2)
+       |(Prec (x,p1),Pint b) -> let p2=Array.map copyP p1 in
+           p2.(0)<- plusP p1.(0) q;
+           Prec (x,p2)
+       |(Prec (x,p1),Prec (y,q1)) ->
+          if x<y then (let q2=Array.map copyP q1 in
+                         q2.(0)<- plusP p q1.(0);
+                         Prec (y,q2))
+          else if x>y then (let p2=Array.map copyP p1 in
+                              p2.(0)<- plusP p1.(0) q;
+                              Prec (x,p2))
+          else
+            (let n=max (deg x p) (deg x q) in
+             let r=Array.create (n+1) (Pint coef0) in
+               for i=0 to n do
+                 r.(i)<- plusP (coef x i p) (coef x i q);
+               done;
+               Prec(x,r)))
+  in norm res
+
+
+(* content, positive integer *)
+let rec content p =
+  match p with
+      Pint a -> C.abs a
+    | Prec (x ,p1) ->
+       Array.fold_left C.pgcd coef0 (Array.map content p1)
+
+let rec div_int p a=
+  match p with
+      Pint b -> Pint (C.div b a)
+    | Prec(x,p1) -> Prec(x,Array.map (fun x -> div_int x a) p1)
+
+let vire_contenu p =
+  let c = content p in
+    if C.equal c coef0 then p else div_int p c
+
+(* sorted list of variables of a polynomial *)
+
+let rec vars=function
+    Pint _->[]
+  | Prec (x,l)->(List.flatten ([x]::(List.map vars (Array.to_list l))))
+
+
+(* multiply p by 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))
+
+
+(* product *)
+let rec multP p q =
+  match (p,q) with
+      (Pint a,Pint b) -> Pint (C.mult a b)
+    |(Pint a, Prec (y,q1)) ->
+       if C.equal a coef0 then Pint coef0
+       else let q2 = Array.map (fun z-> multP p z) q1 in
+         Prec (y,q2)
+
+    |(Prec (x,p1), Pint b) ->
+       if C.equal b coef0 then Pint coef0
+       else let p2 = Array.map (fun z-> multP z q) p1 in
+         Prec (x,p2)
+    |(Prec (x,p1), Prec(y,q1)) ->
+       if x<y
+       then (let q2 = Array.map (fun z-> multP p z) q1 in
+               Prec (y,q2))
+       else if x>y
+       then (let p2 = Array.map (fun z-> multP z q) p1 in
+               Prec (x,p2))
+       else Array.fold_left plusP (Pint coef0)
+         (Array.mapi (fun i z-> (multx i x (multP z q))) p1)
+
+
+
+(* derive p with 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
+
+
+(* opposite *)
+let rec oppP p =
+  match p with
+      Pint a -> Pint (C.opp a)
+    |Prec(x,p1) -> Prec(x,Array.map oppP p1)
+
+let moinsP p q=plusP p (oppP q)
+
+let rec puisP p n = match n with
+    0 -> cf1
+  |_ -> (multP p (puisP p (n-1)))
+
+
+(* infix notations *)
+(*let (++) a b = plusP a b
+*)
+let (@@) a b = multP a b
+
+let (--) a b = moinsP a b
+
+let (^^) a b = puisP a b
+
+
+(* leading coefficient in v,  v>= max_var p *)
+
+let coefDom v p= coef v (deg v p) p
+
+let coefConst v p = coef v 0 p
+
+(* tail of a polynomial *)
+let remP v p =
+  moinsP p (multP (coefDom v p) (puisP (x v) (deg v p)))
+
+
+(* first interger coefficient of 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)
+
+
+(* divide by the content and make the head int coef positive *)
+let normc p =
+  let p = vire_contenu p in
+  let a = coef_int_tete p in
+    if C.le coef0 a then p else oppP p
+
+
+(* constant coef of normalized polynomial *)
+let rec coef_constant p =
+  match p with
+      Pint a->a
+    |Prec(_,q)->coef_constant q.(0)
+
+
+(***********************************************************************
+  3. Printing polynomials.
+*)
+
+(* if univ = false, we use x,y,z,a,b,c,d... as variables, else x1,x2,...
+*)
+let univ=ref true
+
+let string_of_var x=
+  if !univ then
+    "u"^(string_of_int x)
+  else
+    if x<=3 then String.make 1 (Char.chr(x+(Char.code 'w')))
+    else String.make 1 (Char.chr(x-4+(Char.code 'a')))
+
+let nsP = ref 0
+
+let rec string_of_Pcut p =
+  if (!nsP)<=0
+  then "..."
+  else
+  match p with
+  |Pint a-> nsP:=(!nsP)-1;
+      if C.le coef0 a
+      then C.to_string a
+      else "("^(C.to_string a)^")"
+  |Prec (x,t)->
+      let v=string_of_var x
+      and s=ref ""
+      and sp=ref "" in
+    let st0 = string_of_Pcut t.(0) in
+      if st0<>"0"
+      then s:=st0;
+    let fin = ref false in
+    for i=(Array.length t)-1 downto 1 do
+      if (!nsP)<0
+      then (sp:="...";
+	    if not (!fin) then s:=(!s)^"+"^(!sp);
+	    fin:=true)
+      else (
+	let si=string_of_Pcut  t.(i) in
+	sp:="";
+	if i=1
+	then (
+	  if si<>"0"
+	  then (nsP:=(!nsP)-1;
+		if si="1"
+		then sp:=v
+		else
+		  (if (String.contains si '+')
+		  then sp:="("^si^")*"^v
+		  else sp:=si^"*"^v)))
+	else (
+	  if si<>"0"
+	  then (nsP:=(!nsP)-1;
+		if si="1"
+		then sp:=v^"^"^(string_of_int i)
+		else (if (String.contains si '+')
+		then sp:="("^si^")*"^v^"^"^(string_of_int i)
+		else  sp:=si^"*"^v^"^"^(string_of_int i))));
+	if !sp<>"" && not (!fin)
+	then (nsP:=(!nsP)-1;
+	      if !s=""
+	      then s:=!sp
+	      else s:=(!s)^"+"^(!sp)));
+    done;
+    if !s="" then (nsP:=(!nsP)-1;
+		   (s:="0"));
+    !s
+
+let to_string p =
+  nsP:=20;
+  string_of_Pcut p
+
+let printP p = Format.printf "@[%s@]" (to_string p)
+
+let print_tpoly lp =
+  let s = ref "\n{ " in
+    Array.iter (fun p -> s:=(!s)^(to_string p)^"\n") lp;
+    prt0 ((!s)^"}")
+
+let print_lpoly lp = print_tpoly (Array.of_list lp)
+
+(***********************************************************************
+  4. Exact division of polynomials.
+*)
+
+(* return (s,r) s.t. p = s*q+r *)
+let rec quo_rem_pol p q x =
+  if x=0
+  then (match (p,q) with
+          |(Pint a, Pint b) ->
+	     if C.equal (C.modulo a b) coef0
+             then (Pint (C.div a b), cf0)
+             else failwith "div_pol1"
+	  |_ -> assert false)
+  else
+    let m = deg x q in
+    let b = coefDom x q in
+    let q1 = remP x q in (* q = b*x^m+q1 *)
+    let r = ref p in
+    let s = ref cf0 in
+    let continue =ref true in
+      while (!continue) && (not (equal !r cf0)) do
+	let n = deg x !r in
+	  if n<m
+	  then continue:=false
+	  else (
+            let a = coefDom x !r in
+            let p1 = remP x !r in  (* r = a*x^n+p1 *)
+            let c = div_pol a b (x-1) in  (* a = c*b *)
+	    let s1 = c @@ ((monome x (n-m))) in
+              s:= plusP (!s) s1;
+              r:= p1 -- (s1 @@ q1);
+          )
+      done;
+      (!s,!r)
+
+(* returns quotient p/q if q divides p, else fails *)
+and div_pol p q x =
+  let (s,r) = quo_rem_pol p q x in
+    if equal r cf0
+    then s
+    else  failwith ("div_pol:\n"
+		   ^"p:"^(to_string p)^"\n"
+		   ^"q:"^(to_string q)^"\n"
+		   ^"r:"^(to_string r)^"\n"
+		   ^"x:"^(string_of_int x)^"\n"
+		   )
+let divP p q=
+  let x = max (max_var_pol p) (max_var_pol q) in
+  div_pol p q x
+
+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 and gcd with subresultants.
+*)
+
+(* pseudo division :
+   q = c*x^m+q1
+   retruns (r,c,d,s) s.t. 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)
+	)
+
+(* gcd with subresultants *)
+
+let rec pgcdP p q =
+  let x = max (max_var_pol p) (max_var_pol q) in
+    pgcd_pol p q x
+
+and pgcd_pol p q x =
+  pgcd_pol_rec p q x
+
+and content_pol p x =
+  match p with
+      Prec(v,p1) when v=x ->
+        Array.fold_left (fun a b -> pgcd_pol_rec a b (x-1)) cf0 p1
+    | _ -> p
+
+and pgcd_coef_pol c p x =
+  match p with
+      Prec(v,p1) when x=v ->
+        Array.fold_left (fun a b -> pgcd_pol_rec a b (x-1)) c  p1
+    |_ -> pgcd_pol_rec c p (x-1)
+
+and pgcd_pol_rec p q x =
+ match (p,q) with
+	(Pint a,Pint b) -> Pint (C.pgcd (C.abs a) (C.abs b))
+      |_ ->
+	  if equal p cf0
+	  then q
+	  else if equal q cf0
+	  then p
+	  else if (deg x q) = 0
+	  then pgcd_coef_pol q p x
+	  else if (deg x p) = 0
+	  then pgcd_coef_pol p q x
+	  else (
+	    let a = content_pol p x in
+	    let b = content_pol q x in
+	    let c = pgcd_pol_rec a b (x-1) in
+	    pr (string_of_int x);
+	    let p1 = div_pol p c x in
+	    let q1 = div_pol q c x in
+	    let r = gcd_sub_res p1 q1 x in
+	    let cr = content_pol r x in
+	    let res = c @@ (div_pol r cr x) in
+	    res
+	   )
+
+(* Sub-résultants:
+
+   ai*Ai = Qi*Ai+1 + bi*Ai+2
+
+   deg Ai+2 < deg Ai+1
+
+   Ai = ci*X^ni + ...
+   di = ni - ni+1
+
+   ai = (- ci+1)^(di + 1)
+   b1 = 1
+   bi = ci*si^di  si i>1
+
+   s1 = 1
+   si+1 = ((ci+1)^di*si)/si^di
+
+*)
+and gcd_sub_res p q x =
+  if equal q cf0
+  then p
+  else
+    let d = deg x p in
+    let d' = deg x q in
+      if d<d'
+      then gcd_sub_res q p x
+      else
+	let delta = d-d' in
+	let c' = coefDom x q in
+	let r = snd (quo_rem_pol (((oppP c')^^(delta+1))@@p) (oppP q) x) in
+	  gcd_sub_res_rec q r (c'^^delta) c' d' x
+
+and gcd_sub_res_rec p q s c d x =
+  if equal q cf0
+  then p
+  else (
+    let d' = deg x q in
+    let c' = coefDom x q in
+    let delta = d-d' in
+    let r = snd (quo_rem_pol (((oppP c')^^(delta+1))@@p) (oppP q) x) in
+    let s'= lazard_power c' s delta x in
+      gcd_sub_res_rec q (div_pol r (c @@ (s^^delta)) x) s' c' d' x
+  )
+
+and lazard_power c s d x =
+  let res = ref c in
+    for i=1 to d-1 do
+      res:= div_pol ((!res)@@c) s x;
+    done;
+    !res
+
+(* memoizations *)
+
+let rec hash = function
+    Pint a -> (C.hash a)
+  | Prec (v,p) ->
+      Array.fold_right (fun q h -> h + hash q) p 0
+
+module Hashpol = Hashtbl.Make(
+  struct
+    type poly = t
+    type t = poly
+    let equal = equal
+    let hash = hash
+  end)
+
+end