(* rationals, i.e. fractions of multivariate polynomials over the real field

   author: Diana Meindl

   (c) Diana Meindl

   Use is subject to license terms. 

1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890

        10        20        30        40        50        60        70        80         90      100

*)

theory GCD_Poly imports Complex_Main begin 

chapter {* The Isabelle Theory *}

text {*

  Below we reference 

  F. Winkler, Polynomial Algorithms. ...

  by page 93 and 95. This is the final version of the Isabelle theory,

  which conforms with Isabelle's coding standards and with requirements

  in terms of logics appropriate for a theorem prover.

*}

section {* Predicates for the specifications *}

text {*

  This thesis provides program code to be integrated into the Isabelle 

  distribution.  Thus the code will be verified against specifications, and 

  consequently the code, as \emph{explicit} specification of how to calculate 

  results, is accompanied by \emph{implicit} specifications describing the 

  properties of the program code.

  For describing the properties of code the format of Haskell will be used.

  Haskell \cite{TODO} is front-of-the-wave in research on functional 

  programming and as such advances beyond SML. The advancement relevant for 

  this thesis is the concise combination of explicit and implicit specification

  of a function. The final result of the code developed within this thesis is 

  described as follows in Haskell format:

  \begin{verbatim} % TODO: find respective Isabelle/LaTeX features

    gcd_poly :: int poly \<Rightarrow> int poly \<Rightarrow> int poly

    gcd_poly p1 p2 = d

    assumes and p1 \<noteq> [:0:] and p2 \<noteq> [:0:]

    yields d dvd p1 \<and> d dvd p2 \<and> \<forall>dd. (dd dvd p1 \<and> dd dvd p2) \<longrightarrow> dd \<le> d

  \end{verbatim

  The above specification combines Haskell format with the term language of

  Isabelle\HOL. The latter provides a large and ever growing library of 

  mechanised mathematics knowledge, which this thesis re-uses in order to 

  prepare for integration of the code.

  In the sequel there is a brief account of the notions required for the 

  specification of the functions implemented within this thesis. *}

subsection {* The Isabelle types required *}

text {*

  bool, nat, int, option, THE, list, poly

*}

subsection {* The connectives required *}

text {*

  logic ~~/doc/tutorial.pdf p.203

  arithmetic < \<le> + - * ^ (*\<equiv> Power.power FIXME: make uniform*) div \<bar> \<bar>

*}

subsection {* The Isabelle functions required *}

text {*

  number theory

    dvd, mod, gcd

    Primes.prime

    Primes.next_prime_bound

  lists

    List.hd, List.tl

    List.length, List.member

    List.fold, List.map

    List.nth, List.take, List.drop

    List.replicate

  others

    Fact.fact

*}

section {* gcd for univariate polynomials *}

ML {*

type unipoly = int list;

val a = [~18, ~15, ~20, 12, 20, ~13, 2];

val b = [8, 28, 22, ~11, ~14, 1, 2];

*}

subsection {* a variant for div *}

ML {* (* 5 div 2 = 2; ~5 div 2 = ~3; but ~5 div2 2 = ~2; *)

  infix div2

  fun a div2 b = 

    if a div b < 0 

    then (abs a) div (abs b) * ~1

    else a div b;

*}

subsection {* modulo calculations for integers *}

ML {*

  (* mod_inv :: int \<Rightarrow> nat 

     mod_inv r m = s

     assumes True

     yields r * s mod m = 1 *)

  fun mod_inv r m =

    let

      fun modi (r, rold, m, rinv) = 

        if m <= rinv then 0 else

          if r mod m = 1

          then rinv

          else modi (rold * (rinv + 1), rold, m, rinv + 1)

     in modi (r, r, m, 1) end;

  (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat

     mod_div a b m = c

     assumes True

     yields a * b ^(-1) mod m = c   <\<Longrightarrow>  a mod m = *)

  fun mod_div a b m = a * (mod_inv b m) mod m;

  (* required just for one approximation:

     approx_root :: nat \<Rightarrow> nat

     approx_root a = r

     assumes True

     yields r * r \<ge> a *)

  fun approx_root a =

    let

      fun root1 a n = if n * n < a then root1 a (n + 1) else n

    in root1 a 1 end

  (* chinese_remainder :: int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int

     chinese_remainder r1 m1 r2 m2 = r

     assumes True

     yields r = r1 mod m1 \<and> r = r2 mod m2 *)

  fun chinese_remainder (r1: int, m1: int, r2: int, m2: int ) =

   (r1 mod m1) + ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1

*}

subsection {* operations with primes *}

ML{*

  (* is_prime :: int list \<Rightarrow> int \<Rightarrow> bool

     is_prime ps n = b

     assumes max ps < n  \<and>  n \<le> (max ps)^2  \<and>  

       (*     FIXME: really ^^^^^^^^^^^^^^^? *)

       (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and> (\<forall>p. p < n \<and> Primes.prime p \<longrightarrow> List.member ps p)

     yields Primes.prime n *)

  fun is_prime [] _ = true

    | is_prime prs n = 

        if n mod (List.hd prs) <> 0 then is_prime (List.tl prs) n else false

  (* make_primes is just local to primes_upto only:

     make_primes :: int list \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int list

     make_primes ps last_p n = pps

     assumes last_p = maxs ps  \<and>  (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p)  \<and>

       (\<forall>p. p < last_p \<and> Primes.prime p \<longrightarrow> List.member ps p)

     yields n \<le> maxs pps  \<and>  (\<forall>p. List.member pps p \<longrightarrow> Primes.prime p)  \<and>

       (\<forall>p. (p < maxs pps \<and> Primes.prime p) \<longrightarrow> List.member pps p)*)

    fun make_primes ps last_p n =

        if n <= (nth ps (List.length ps - 1)) then ps else

          if is_prime ps (last_p + 2) 

          then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n

          else make_primes  ps                   (last_p + 2) n

  (* primes_upto :: nat \<Rightarrow> nat list

     primes_upto n = ps

       assumes True

       yields (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>

         n \<le> maxs ps \<and> maxs ps \<le> Fact.fact n + 1 \<and>

         (\<forall>p. p \<le> maxs ps \<and> Primes.prime p \<longrightarrow> List.member ps p) *)

  fun primes_upto n = if n < 3 then [2] else make_primes [2, 3] 3 n

  (* find the next prime greater p not dividing the number n:

    next_prime_not_dvd :: int \<Rightarrow> int \<Rightarrow> int (infix)

    n1 next_prime_not_dvd n2 = p

      assumes True

      yields p is_prime  \<and>  n1 < p  \<and>  \<not> p dvd n2  \<and> (* smallest with these properties... *)

        (\<forall> p'. (p' is_prime  \<and>  n1 < p'  \<and>  \<not> p' dvd n2)  \<longrightarrow>  p \<le> p') *)

  infix next_prime_not_dvd;

  fun n1 next_prime_not_dvd n2 = 

    let

      val ps = primes_upto (n1 + 1)

      val nxt = nth ps (List.length ps - 1);

    in

      if n2 mod nxt <> 0 then nxt else nxt next_prime_not_dvd n2

  end

*}

subsection {* auxiliary functions for univariate polynomials *}

ML {*

  (* leading coefficient *)

  fun lcoeff (p: unipoly) =

    if nth p (List.length p - 1) <> 0

    then nth p (List.length p - 1)

    else lcoeff (nth_drop (length p - 1) p : unipoly);

  (* degree *)

  fun deg (p: unipoly) = 

    if nth p (length p - 1) <> 0

    then length p - 1

    else deg (nth_drop (length p - 1) p : unipoly);

  (* drop leading coefficients equal 0 *)

  fun drop_lc0 [] = []

    | drop_lc0 (p: unipoly) = 

      if nth p (length p - 1) = 0

      then drop_lc0 (nth_drop (length p - 1) p : unipoly)

      else p;

  (* normalise a unipoly such that lcoeff mod m = 1.

     normalise :: unipoly \<Rightarrow> nat \<Rightarrow> unipoly

     normalise [p_0, .., p_k] m = [q_0, .., q_k]

       assumes True

       yields \<exists> t. 0 \<le> i \<le> k \<Rightarrow> (p_i * t) mod m = q_i  \<and>  (p_k * t) mod m = 1 *)

  fun normalise (p: unipoly) (m: int) =

    let

      val p = drop_lc0 p

      val lcp = lcoeff p;

      fun norm nrm i =

        if i = 0

        then [mod_div (nth p i) lcp m] @ nrm

        else norm ([mod_div (nth p i) lcp m] @  nrm) (i - 1) ;

    in 

      if List.length p = 0 then p else norm [] (List.length p - 1)

  end;

  (* scalar multiplication *)

  infix %*

  fun (p: unipoly) %* n =  map (Integer.mult n) p : unipoly

  (* scalar divison *)

  infix %/ 

  fun (p: unipoly) %/ n =

    let fun swapf f a b = f b a (* swaps the arguments of a function *)

    in map (swapf (curry op div2) n) p : unipoly end;

  (* minus of polys *)

  infix %-%

  fun (p1: unipoly) %-% (p2: unipoly) = 

    map2 (curry op -) p1 (p2 @ replicate (List.length p1 - List.length p2) 0) : unipoly

  (* dvd for univariate polynomials *)

  infix %|%

  fun [d] %|% [p] =

      if abs d <= abs p andalso p mod d = 0 then true else false

    | (ds: unipoly) %|% (ps: unipoly) = 

      let 

        val ds = drop_lc0 ds;

        val ps = drop_lc0 ps;

        val d000 = (replicate (List.length ps - List.length ds) 0) @ ds ;

        val quot = (lcoeff ps) div2 (lcoeff d000);

        val rest = drop_lc0 (ps %-% (d000 %* quot));

      in 

        if rest = [] then true else

          if quot <> 0 andalso List.length ds <= List.length rest

          then ds %|% rest 

          else false

      end

  (* normalise :: poly_centr \<Rightarrow> unipoly => int \<Rightarrow> unipoly

     normalise [p_0, .., p_k] m = [q_0, .., q_k]

       assumes 

       yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^| 

        (where |^ x ^| means round x up to the next greater integer) *)

  fun poly_centr (p: unipoly) (m: int) =

    let

      val mi = m div2 2;

      val mid = if m mod mi = 0 then mi else mi + 1;

      fun centr m mid p_i = if mid < p_i then p_i - m else p_i;

    in map (centr m mid) p : unipoly end;

  (*** landau mignotte bound (Winkler p91)

    every coefficient of the gcd of a  and b in Z[x] is bounded in absolute value by

    2^{min(m,n)} * gcd(a_m,b_n) * min(1/|a_m| 

    * root(sum_{i=0}^{m}a_i^2) ,1/|b_n| * root( sum_{i=0}^{n}b_i^2 ) ***)

   (*  sum_lmb :: poly_centr \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int

     sum_lmb [p_0, .., p_k] e = s

       assumes exp >= 0

       yields. p_0^e + p_1^e + ... + p_k^e *)

  fun sum_lmb (p: unipoly) exp = fold ((curry op +) o (Integer.pow exp)) p 0;

   (* LANDAU_MIGNOTTE_bound :: poly_centr \<Rightarrow> unipoly \<Rightarrow> unipoly \<Rightarrow> int

      LANDAU_MIGNOTTE_bound [a_0, ..., a_m] [b_0, .., b_n] = lmb

        assumes True

        yields lmb = 2^min(m,n) * gcd(a_m,b_n) * 

          min( 1/|a_m| * root(sum_lmb [a_0,...a_m] 2 , 1/|b_n| * root(sum_lmb [b_0,...b_n] 2)*)

  fun LANDAU_MIGNOTTE_bound (p1: unipoly) (p2: unipoly) = 

    (Integer.pow (Integer.min (deg p1) (deg p2)) 2) * (abs (Integer.gcd (lcoeff p1) (lcoeff p2))) * 

    (Int.min (abs((approx_root (sum_lmb p1 2)) div2 ~(lcoeff p1)), 

              abs((approx_root (sum_lmb p2 2)) div2 ~(lcoeff p2))));

*}

subsection {* modulo calculations for polynomials *}

ML{*

  (* chinese_remainder_poly :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly

     chinese_remainder_poly (m1, m2) (p1, p2) = p

     assume m1, m2 relatively prime

     yields p1 = p mod m1 \<and> p2 = p mod m2

 *)

  fun chinese_remainder_poly (m1, m2) (p1: unipoly, p2: unipoly) = 

  let 

    fun chinese_rem (m1, m2) (p_1i, p_2i) = chinese_remainder (p_1i, m1, p_2i, m2)

  in map (chinese_rem (m1, m2)) (p1 ~~ p2): unipoly end

  (* mod_poly :: unipoly \<Rightarrow> int \<Rightarrow> unipoly

     mod_poly [p1, p2, ..., pk] m = up 

     assume m > 0

     yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)

  infix mod_poly 

    fun uvp mod_poly m = 

    let 

      fun modp m up = (curry op mod) up m 

    in map (modp m) uvp end

  (* euclidean algoritm in Z_p[x].

     mod_poly_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly

     mod_poly_gcd p1 p2 m = g

       assumes True

       yields gcd ((p1 mod m), (p2 mod m)) = g *)

  fun mod_poly_gcd (p1: unipoly) (p2: unipoly) (m: int) =

    let 

      val p1m = p1 mod_poly m;

      val p2m = drop_lc0 (p2 mod_poly m);

      val p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;

      val quot =  mod_div (lcoeff p1m) (lcoeff p2n) m;

      val rest = drop_lc0 ((p1m %-% (p2n %* quot)) mod_poly m)

    in

      if rest = [] then p2 else

        if length rest < List.length p2

        then mod_poly_gcd p2 rest m 

        else mod_poly_gcd rest p2 m

    end;

  (* prim_poly :: unipoly \<Rightarrow> unipoly

     prim_poly p = pp

     assumes True

     yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)

  fun prim_poly [c: int] = if c = 0 then [0] else [1]

    | prim_poly (p: unipoly) =

        let

          val d = abs (Integer.gcds p)

        in

          if d = 1 then p else p %/ d

        end

*}

subsection {* gcd_unipoly, code from [1] p.93 *}

ML {*

  (* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int  \<Rightarrow> unipoly

     try_new_prime_up p1 p2 d M P g p = new_g

       assumes d = gcd(lc(p1), lc(p2)) \<and> M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>

               p1 is primitive \<and> p2 is primitive  ----> primitive = prim_poly

       yields  new_g = [1] \<or> (new_g \<ge> g \<and> P > M)*)

  fun try_new_prime_up a b d M P g p =

    if P > M then g else

      let

        val p = p next_prime_not_dvd d

        val g_p = 

          poly_centr (((normalise (mod_poly_gcd a b p) p) %* (d mod p)) mod_poly p) p: unipoly

      in

        if deg g_p < deg g

        then

          if (deg g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p

        else

          if deg g_p <> deg g then try_new_prime_up a b d M P g p else

            let 

              val P = P * p

              val g = poly_centr ((chinese_remainder_poly (P, p) (g, g_p)) mod_poly P) P

            in try_new_prime_up a b d M P g p end

      end

  (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly

     HENSEL_lifting_up p1 p2 d M p = gcd

       assumes d = gcd(lc(p1), lc(p2)) \<and> M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d  \<and>

               p1 is primitive \<and> p2 is primitive  ----> primitive = prim_poly

       yields  gcd|p1  \<and> gcd|p2  \<and> gcd is primitive*)

  fun HENSEL_lifting_up a b d M p = 

    let 

      val p = p next_prime_not_dvd d 

      val g_p = poly_centr(((normalise ( mod_poly_gcd a b p) p) %* (d mod p)) mod_poly p) p

    in

      if deg g_p = 0 then [1] else

        let 

          val g = prim_poly (try_new_prime_up a b d M p g_p p)

        in

          if (g %|% a) andalso (g %|% b) then g else HENSEL_lifting_up a b d M p

        end

    end

  (* gcd_poly :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly

     gcd_poly a b = c

     assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and> 

             a is primitive \<and> b is primitive

     yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)

  fun gcd_unipoly (a: unipoly) (b: unipoly) =

    if a = b then a else

      HENSEL_lifting_up a b (abs (Integer.gcd (lcoeff a) (lcoeff b)))  

        (2 * (abs (Integer.gcd (lcoeff a) (lcoeff b))) * LANDAU_MIGNOTTE_bound a b) 1;

*}

(* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)

subsection {* gcd_poly Algorgithmus, auxiliary functions multivariate *}

(* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)

ML {*

type monom = (int * int list);

type multipoly = monom list;

*}

subsection {* calculations for multivariate polynomials *}

ML {*

(* assumes a and b have the same length*)

  fun listgreater a b = 

    let fun gr a b = a>=b;

    fun all [] = true 

      | all list = 

      if List.hd list then all (List.tl list)  else false 

    in all (map2 gr a b)

  end

  fun get_coef (mvp: multipoly) (place: int) =

    let fun coef ((c,_): monom) = c;

    in coef (List.nth (mvp, place)) end 

  fun get_varlist (mvp: multipoly) (place: int) =

    let fun varlist ((_,var): monom) = var;

    in varlist (nth mvp place) end

  fun get_coeflist (poly: multipoly) = 

    let 

      fun get_clist (poly: multipoly) list = 

        if poly = [] then list

        else get_clist (List.tl poly) (list @ [get_coef poly 0])

    in get_clist poly [] end

*}

ML{*

(* like 3x+4z = [(3,[1,0]),(4,[0,1])] with x<z

add_var y with order 2  =>  3x+4y = [(3,[1,0,0]),(4,[0,0,1])] with x<y<z*)

(* OFIZIELLE BESCHREIBUNG!!*)

  fun add_var_to_multipoly (mvp: multipoly) (order: int) =

    let fun add ([]: multipoly) (_: int) (new: multipoly) = new

          | add (mvp: multipoly) (order: int) (new: multipoly) =

          let val (first, last) = chop order (get_varlist mvp 0)

          in add (List.tl mvp) (order) (new @ [((get_coef mvp 0),(first @ [0] @ last))]) end 

    in add mvp order [] end

  fun cero_multipoly (different_var: int) = [(0, replicate different_var 0)];

*}

ML{*

  (* if a is greater than b --> true or false*)

  fun  greater_var (a: int list) (b: int list) =

   if drop_lc0 (a %-% b) = [] then false else lcoeff (drop_lc0 (a %-% b)) > 0 ;

*}

ML{*

   (* drop leading coefficients equal 0 *)

  fun drop_lc0_multipoly [m: monom] = 

    if get_coef [m] (length [m] - 1) = 0

    then cero_multipoly (length (get_varlist [m] 0))

    else [m]

    | drop_lc0_multipoly (mvp: multipoly) = 

      if get_coef mvp (length mvp - 1) = 0  

      then drop_lc0_multipoly (nth_drop (length mvp - 1) mvp)

      else mvp

  (* add same variables *)

  fun short_multipoly (mvp: multipoly) = 

    let fun short (mvp: multipoly) (new: multipoly) =

              if length mvp =1 

              then if (get_coef mvp 0) = 0 

                   then if new = [] then cero_multipoly (length (get_varlist mvp 0)) else new

                   else  new @ mvp 

              else if get_varlist  mvp 0 = get_varlist mvp 1

                then short ( [(get_coef mvp 0 + get_coef mvp 1,get_varlist  mvp 0)]

                             @ nth_drop 0 (nth_drop 0 mvp))   new 

                else if (get_coef mvp 0) = 0 then short (nth_drop 0 mvp) new 

                  else short (nth_drop 0 mvp) (new @ [nth mvp 0])

    in short mvp [] end

  (* brings the polynom in correct order and add same variables*)

  fun order_multipoly (a: multipoly)=

    let 

      val ordered = [nth a 0];

      fun order_mp [] [] = cero_multipoly (length (get_varlist a 0))

        | order_mp [] (ordered: multipoly) = short_multipoly ordered

        | order_mp a (ordered: multipoly) =

        if  greater_var (get_varlist a 0) (get_varlist ordered (length ordered -1))

          then order_mp (nth_drop 0 a)(ordered @ [nth a 0])

          else let 

            val rest = [nth ordered (length ordered -1)];

            fun order_mp' [] (new: multipoly) (rest: multipoly) = new @ rest

              | order_mp' (ordered': multipoly) (new: multipoly) (rest: multipoly) =

                if greater_var (get_varlist new 0) (get_varlist ordered' (length ordered' -1))

                  then ordered' @ new @ rest

                  else order_mp' (nth_drop (length ordered' -1) ordered') new 

                                 ([nth ordered' (length ordered' -1)]@ rest)

            in order_mp (nth_drop 0 a) 

                        (order_mp' (nth_drop (length ordered -1) ordered) [nth a 0] rest)

          end

    in order_mp (nth_drop 0 a) ordered

    end

 (* greatest variablegroup  *)

  fun deg_multipoly (mvp: multipoly) = 

    get_varlist (order_multipoly mvp) (length (order_multipoly mvp) -1)

  fun  max_deg_var [m: monom] (x: int) = nth (get_varlist [m] 0) x |

   max_deg_var (mvp: multipoly) (x: int) =

    let

      fun max_monom (m1: monom) (m2: monom) (x: int)=

        if nth (get_varlist [m1] 0) x < nth (get_varlist [m2] 0) x  then m2  else m1;

    in

    max_deg_var ([max_monom (nth(mvp) 0)( nth(mvp) 1) x] @ nth_drop 0 (nth_drop 0 mvp)) x

    end

  infix %%/

  fun (poly: multipoly) %%/ b = (* =quotient*)

    let fun division monom b = (((get_coef [monom] 0) div2 b),get_varlist [monom] 0);

    in order_multipoly(map2 division poly (replicate (length(poly)) b)) end

  infix %%*

    fun (poly: multipoly) %%* (b: int) = 

      let fun mult monom b = (((get_coef [monom] 0) * b),get_varlist [monom] 0);

      in order_multipoly(map2 mult poly (replicate (length(poly)) b)) end

  fun cont [a: int] = a

    | cont (poly1: unipoly) = abs (Integer.gcds poly1)

  fun evaluate_multipoly (mvp: multipoly) (var: int) (value: int) = 

    let fun eval ([]: multipoly) (_: int) (_: int) (new: multipoly) = order_multipoly new |

      eval (mvp: multipoly) (var: int) (value: int) (new: multipoly) =

      eval (nth_drop 0 mvp) var value

           (new @  [((Integer.pow  (nth (get_varlist mvp 0) var)value) * get_coef mvp 0, 

                    nth_drop var (get_varlist mvp 0))]);

    in eval mvp var value [] end

  fun multi_to_uni (mvp: multipoly) = 

    if length (get_varlist mvp 0) = 1 

    then let fun mtu ([]: multipoly) (uvp: unipoly) = uvp |

                 mtu (mvp: multipoly) (uvp: unipoly) =

               if length uvp = (nth(get_varlist mvp 0) 0)

               then mtu (nth_drop 0 mvp) (uvp @ [get_coef mvp 0])

               else mtu mvp (uvp @ [0])    

         in mtu (order_multipoly mvp) [] end

    else error "Polynom has more than one variable!";

  fun uni_to_multi (uvp: unipoly) =

    let fun utm ([]: unipoly) (mvp: multipoly) (_: int)=  short_multipoly mvp

          | utm (uvp: unipoly) (mvp: multipoly) (counter: int) = 

          utm (nth_drop 0 uvp) (mvp @ [((nth uvp 0),[counter])]) (counter+1)

    in  utm uvp [] 0 end

  fun mult_with_new_var ([]: multipoly) (_: unipoly) (_: int) = []

    | mult_with_new_var (mvp: multipoly) (uvp: unipoly) (order: int) = 

    let 

      fun mult ([]: multipoly) (_: unipoly) (_: int) (new: multipoly) (_: int) = 

            short_multipoly new

        | mult (mvp: multipoly) (uvp: unipoly) (order: int) (new: multipoly) (_: int)  =

        let 

          fun mult' (_: multipoly) ([]: unipoly) (_) (new': multipoly) (_) = order_multipoly new'

            | mult' (mvp': multipoly) (uvp': unipoly) (order: int) (new': multipoly) (c2': int) =

            let val (first, last) = chop order (get_varlist mvp' 0)

            in mult' mvp' (nth_drop 0 uvp') order

                    (new' @ [(((get_coef mvp' 0) * (nth uvp' 0)),(first @ [c2'] @ last))]) (c2'+1)

            end

        in mult (nth_drop 0 mvp) uvp order (new @ (mult' mvp uvp order [] 0)) (0) 

        end

    in mult mvp uvp order [] 0 end

  fun greater_deg_multipoly (var1: int list) (var2: int list) =

    if var1 = [] andalso var2 =[] then 0

    else if (nth var1 (length var1 -1)) = (nth var2 (length var1 -1) ) 

         then greater_deg_multipoly (nth_drop (length var1 -1) var1) (nth_drop (length var1 -1) var2)

         else if (nth var1 (length var1 -1)) > (nth var2 (length var1 -1)) 

              then 1 else 2 ;

  infix %%+%%

  fun ([]: multipoly) %%+%% (mvp2: multipoly) = mvp2

    | (mvp1: multipoly) %%+%% ([]: multipoly) = mvp1

    | (mvp1: multipoly) %%+%% (mvp2: multipoly) =

    let fun plus ([]: multipoly) (mvp2: multipoly) (new: multipoly) = order_multipoly (new @ mvp2)

          | plus (mvp1: multipoly) ([]: multipoly) (new: multipoly) = order_multipoly (new @ mvp1)

          | plus (mvp1: multipoly) (mvp2: multipoly) (new: multipoly) = 

          if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 0

          then plus (nth_drop 0 mvp1) (nth_drop 0 mvp2) 

                    (new @ [(((get_coef mvp1 0) + (get_coef mvp2 0)), get_varlist mvp1 0)])

          else if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 1

               then plus mvp1 (nth_drop 0 mvp2) (new @ [nth mvp2 0])

               else plus (nth_drop 0 mvp1) mvp2 (new @ [nth mvp1 0])

    in plus mvp1 mvp2 [] end

  infix %%-%%

  fun (mvp1: multipoly) %%-%% (mvp2: multipoly) = 

    let fun neg ([]: multipoly) (new: multipoly) = order_multipoly new

          | neg (mvp: multipoly) (new: multipoly) =

          neg (nth_drop 0 mvp) (new @ [(((get_coef mvp 0) * ~1), get_varlist mvp 0)])

        val neg_mvp2 = neg mvp2

    in mvp1 %%+%% (neg_mvp2 [])  end         

  infix %%*%%

  fun (mvp1: multipoly) %%*%% (mvp2: multipoly) =

    let fun mult ([]: multipoly) (_: multipoly) (_: multipoly) (new: multipoly) = 

              order_multipoly new

          | mult (mvp1: multipoly) ([]: multipoly) (regular_mvp2: multipoly) (new: multipoly) =

            mult (nth_drop 0 mvp1) regular_mvp2 regular_mvp2 new

          | mult (mvp1: multipoly) (mvp2: multipoly) (regular_mvp2: multipoly) (new: multipoly) = 

            mult mvp1 (nth_drop 0 mvp2) regular_mvp2 (new @ [(((get_coef mvp1 0) * (get_coef mvp2 0)),

            (map2 Integer.add (get_varlist mvp1 0) (get_varlist mvp2 0)))])

    in if (length mvp1) > (length mvp2) then mult mvp1 mvp2 mvp2 [] else mult mvp2 mvp1 mvp1 [] end

 (* if poly1|poly2 *)

  infix %%|%% 

  fun [(coef2, var2): monom] %%|%%  [(coef1, var1): monom]  = 

    ( (listgreater var1 var2) 

        andalso (coef1 mod coef2) = 0)

    | (_: multipoly) %%|%% [(0,_)]=  true

    | (poly2: multipoly)  %%|%% (poly1: multipoly) =

    if [nth poly2 (length poly2 -1)]   %%|%% [nth poly1 (length poly1 -1)]

    then poly2 %%|%% (poly1 %%-%%

      ([((get_coef poly1 (length poly1 -1)) div2 (get_coef poly2 (length poly2 -1)),

      (get_varlist poly1 (length poly1 -1)) %-%(get_varlist poly2 (length poly2 -1)))] %%*%%

      poly2)) 

    else false

  fun find_new_r (mvp1: multipoly) (mvp2: multipoly) (old: int) =

    let val poly1 = evaluate_multipoly mvp1 (length (get_varlist mvp1 0) - 2) (old + 1)

        val poly2 = evaluate_multipoly mvp2 (length (get_varlist mvp2 0) - 2) (old + 1);

    in

      if max_deg_var poly1 (length (get_varlist poly1 0) - 1)=

          max_deg_var mvp1 (length (get_varlist mvp1 0) - 1) orelse

         max_deg_var poly2 (length (get_varlist poly2 0) - 1)= 

          max_deg_var mvp2 (length (get_varlist mvp2 0) - 1) 

      then old + 1

      else find_new_r (mvp1) (mvp2) (old + 1)

    end

*}

subsection {* Newtoninterpolation *}

ML{* (* first step *)

  fun newton_first (x: int list) (f: multipoly list) (order: int) =

    let val polynom =(add_var_to_multipoly (nth f 0) order) %%+%% 

            ((mult_with_new_var (((nth f 1)%%-%% (nth f 0))%%/

               (nth x 1) - (nth x 0))) [(nth x 0) * ~1, 1] order);

        val new_value_poly =   [(nth x 0) * ~1, 1];

        val steps = [((nth f 1) %%-%%(nth f 0))%%/ ((nth x 1) - (nth x 0))]

    in (polynom, new_value_poly, steps) end

  fun newton_multipoly (x: int list) (f: multipoly list) (steps: multipoly list) 

                       (t: unipoly) (p: multipoly) order = 

    if length x = 2

    then let val (polynom, new_value_poly, steps) =  newton_first x [(nth f 0), (nth f 1)] order

         in (polynom, new_value_poly, steps) end

    else let val new_value_poly = multi_to_uni((uni_to_multi t)  %%*%% 

                                  (uni_to_multi [(nth x (length x -2) )* ~1, 1]));

             val new_steps = [((nth f (length f -1))  %%-%% (nth f (length f -2)))  %%/ 

                               ((nth x (length x - 1)) - (nth x (length x - 2)))];

             fun next_step ([]: multipoly list) (new_steps: multipoly list) (_: int list) = new_steps

               | next_step (steps: multipoly list) (new_steps: multipoly list) (x': int list) =  

                 next_step (nth_drop 0 steps)

                           (new_steps @ [(((nth new_steps (length new_steps -1)) %%-%%(nth steps 0)) %%/

                                    ((nth x' (length x' - 1)) - (nth x' (length x' - 3))))])

                           (nth_drop (length x' -2) x')

             val steps = next_step steps new_steps x;

             val polynom' = p %%+%% 

                            (mult_with_new_var (nth steps (length steps -1)) new_value_poly order);

       in (order_multipoly(polynom'), new_value_poly, steps) end 

*}

ML{*

fun lc_multipoly (mpoly: multipoly) = 

   get_coef (order_multipoly mpoly) (length (order_multipoly mpoly) - 1);

fun bool_list [] = true 

  | bool_list list = if nth list 0 then bool_list (List.tl list) else false;

fun list_smaler (a: int) (list: int list)= 

  let fun s a b = a <= b;

  in map (s a) list

 end

(* assume up2 | up1  and use the same variables *)

infix %%/%% 

fun (up1: multipoly) %%/%% (up2: multipoly) =

    let 

      fun div_up (up1: multipoly) (up2: multipoly) (quot_old: multipoly) = 

      let 

        val up1 = drop_lc0_multipoly (order_multipoly up1);

        val up2 = drop_lc0_multipoly (order_multipoly up2);

        val [c] = cero_multipoly (length (get_varlist up1 0));

        val d = (replicate (List.length up1 - List.length up2) c) @ up2 ;

        val quot = [(lc_multipoly up1 div2 lc_multipoly d,

                    get_varlist up1 (length up1 - 1) %-% get_varlist up2 (length up2 - 1))];

        val rest = drop_lc0_multipoly (up1 %%-%% (d %%*%% quot));

      in if  rest <>  [c] andalso

             bool_list ((list_smaler 0) (get_varlist rest (length rest - 1) %-% get_varlist up2 (length up2 - 1)))

        then div_up rest up2 (quot @ quot_old)

        else (quot @ quot_old, rest)

      end

    in div_up up1 up2 [] 

 end 

*}

subsection {* gcd_poly algorithm, code for [1] p.95: multivariate *}

ML {*

  fun gcd_poly' (a: multipoly) (b: multipoly) n r=

    if greater_var (deg_multipoly b) (deg_multipoly a) then gcd_poly' b a n r else

    if (n-1) = 0 

    then  uni_to_multi((gcd_unipoly (prim_poly(multi_to_uni a)) (prim_poly(multi_to_uni b))) %* 

                       (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))))

    else 

      let 

        val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)-2), 

                             max_deg_var b (length(get_varlist b 0)-2)); 

      in try_new_r a b n M r [] [] end

and try_new_r a b n M r r_list steps = 

     let 

        val r = find_new_r a b r;

        val r_list = r_list @ [r];

        val g_r = gcd_poly' (order_multipoly (evaluate_multipoly a (n-2) r)) 

                            ( order_multipoly (evaluate_multipoly b (n-2) r)) (n-1) 0

        val steps = steps @ [g_r];

      in Hensel_lifting a b n M 1 r r_list steps g_r ( cero_multipoly n ) [1] end

  and Hensel_lifting a b n M m r r_list steps g g_n mult =  

    if m > M then if g_n %%|%% a andalso g_n %%|%% b then  g_n else

      try_new_r a b n M r r_list steps 

    else

      let 

        val r = find_new_r a b r; 

        val r_list = r_list @ [r];

        val g_r = gcd_poly' (order_multipoly  (evaluate_multipoly a (n-2) r)) 

                            (order_multipoly  (evaluate_multipoly b (n-2) r)) (n-1) 0

      in  

        if greater_var  (deg_multipoly g) (deg_multipoly g_r)

        then Hensel_lifting a b n M 1  r r_list steps g g_n mult

        else if (deg_multipoly g)= (deg_multipoly g_r) 

          then let 

            val (g_n, new, steps) = newton_multipoly r_list [g, g_r] steps mult g_n (n-2)

          in 

            if (nth steps (length steps -1)) = (cero_multipoly (n-1))

            then Hensel_lifting a b n M (M+1) r r_list steps g_r g_n new

            else Hensel_lifting a b n M (m+1) r r_list steps g_r g_n new 

          end

          else Hensel_lifting a b n M (m+1) r r_list steps g  g_n mult

      end     

  fun gcd_poly (a: multipoly) (b: multipoly) = 

    let val gcd = gcd_poly' a b (length(get_varlist a 0)) 0;

        val (m1, _) = a %%/%% gcd;

        val (m2, _) = b %%/%% gcd;

    in ((m1, m2), gcd)  end

*}

ML {*

gcd_poly [(6,[0])] [(9,[0])];

(* gcd_poly [(6,[0])] [(0,[0])]; ERROR div with 0*)

(* (gcd_poly "x^2 - y^2"  "x^2 - x*y")-> (("x + y", "x"), "x - y") *)

val p1 = [(1,[2,0]),(~1,[0,2])]; val p2 =  [(1,[2,0]),(~1,[1,1])];

val (a,g) = gcd_poly p1 p2;

val t1 = p1 %%/%% g;

val t2 = p2 %%/%% g;

*}

end