(* 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 

text {*

  Below we reference 

  F. Winkler, Polynomial algorithyms in computer algebra. Springer 1996.

  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>d'. (d' dvd p1 \<and> d' dvd p2) \<longrightarrow> d' \<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, 'a poly

*}

subsection {* The connectives required *}

text {*

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

  arithmetic < \<le> + - * ^ 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 {* auxiliary functions *}

ML {* 

  (* a variant for div: 

    5 div 2 = 2; ~5 div 2 = ~3; BUT WEE NEED ~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

  (* why has "last" been removed since 2002, although "last" is in List.thy ? *)

  fun last ls = (hd o rev) ls

  (* drop tail elements equal 0 *)

  fun drop_hd_zeros (0 :: ps) = drop_hd_zeros ps

    | drop_hd_zeros p = p

  fun drop_tl_zeros p = (rev o drop_hd_zeros o rev) p;

  (* swaps the arguments of a function *)

  fun swapf f a b = f b a

*}

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> natO

     mod_div a b m = c

     assumes True

     yields a * b ^(-1) mod m = c   <\<Longrightarrow>  a mod m = (b * c) 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 * nat \<Rightarrow> int * nat \<Rightarrow> nat

     chinese_remainder (r1, m1) (r2, m2) = r

     assumes True

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

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

    (r1 mod m1) + 

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

*}

subsection {* creation of lists of primes for efficiency *}

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 <= last ps 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 {* basic notions for univariate polynomials *}

ML {*

  (* leading coefficient, includes drop_lc0_up ?FIXME130507 *)

  fun lcoeff_up (p: unipoly) = (last o drop_tl_zeros) p

  (* drop leading coefficients equal 0 *)

  fun drop_lc0_up (p: unipoly) = drop_tl_zeros p : unipoly;

  (* degree, includes drop_lc0_up ?FIXME130507 *)

  fun deg_up (p: unipoly) = ((Integer.add ~1) o length o drop_lc0_up) p;

  (* normalise a unipoly such that lcoeff_up 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_up p

      val lcp = lcoeff_up 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;

*}

subsection {* addition, multiplication, division *}

ML {*

  (* scalar multiplication *)

  infix %*

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

  (* scalar divison *)

  infix %/ 

  fun (p: unipoly) %/ n = map (swapf (curry op div2) n) p : unipoly

  (* 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] = (p mod d = 0)

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

      let 

        val ds = drop_lc0_up ds;

        val ps = drop_lc0_up ps;

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

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

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

      in

        if rest = [] then true else

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

          then ds %|% rest

          else false

      end

*}

subsection {* normalisation and Landau-Mignotte bound *}

ML {* 

  (* normalise :: centr_up \<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 centr_up (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 :: centr_up \<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 :: centr_up \<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_up p1) (deg_up p2)) 2) * (abs (Integer.gcd (lcoeff_up p1) (lcoeff_up p2))) * 

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

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

*}

subsection {* modulo calculations on polynomials *}

ML{*

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

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

     assume m1, m2 relatively prime

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

  fun chinese_remainder_up (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_up :: unipoly \<Rightarrow> int \<Rightarrow> unipoly

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

     assume m > 0

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

  infix mod_up 

  fun (p : unipoly) mod_up m = 

  let 

    fun mod' m p = (curry op mod) p m 

  in map (mod' m) p : unipoly end

  (* euclidean algoritm in Z_p[x].

     mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly

     mod_up_gcd p1 p2 m = g

       assumes True

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

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

    let 

      val p1m = p1 mod_up m;

      val p2m = drop_lc0_up (p2 mod_up m);

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

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

      val rest = drop_lc0_up ((p1m %-% (p2n %* quot)) mod_up m)

    in

      if rest = [] then p2 else

        if length rest < List.length p2

        then mod_up_gcd p2 rest m 

        else mod_up_gcd rest p2 m

    end;

  (* primitive_up :: unipoly \<Rightarrow> unipoly

     primitive_up 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 primitive_up ([c] : unipoly) = 

      if c = 0 then ([0] : unipoly) else ([1] : unipoly)

    | primitive_up p =

        let

          val d = abs (Integer.gcds p)

        in

          if d = 1 then p else p %/ d

        end

*}

subsection {* gcd_up, 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 (lcoeff_up p1, lcoeff_up p2)  \<and> 

               M = LANDAU_MIGNOTTE_bound  \<and>  p = prime \<and>  p ~| d  \<and>  P \<ge> p  \<and>

               p1 is primitive  \<and>  p2 is primitive

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

    argumentns "a b d M P g p" named according to [1] p.93 *)

  fun try_new_prime_up a b d M P g p =

(writeln ("try_new_prime_up:  a = " ^ PolyML.makestring a ^

", b = " ^ PolyML.makestring b ^

", d = " ^ PolyML.makestring d ^

", M = " ^ PolyML.makestring M ^

", P = " ^ PolyML.makestring P ^

", g = " ^ PolyML.makestring g ^

", p = " ^ PolyML.makestring p);

    if P > M then

(writeln ("try_new_prime_up 1 -----> " ^ PolyML.makestring g);

                  g 

)

                    else

      let

        val p = p next_prime_not_dvd d

        val g_p = centr_up (         (    (normalise (mod_up_gcd a b p) p) 

                                       %* (d mod p))                  

                              mod_up p) 

                            p

      in

        if deg_up g_p < deg_up g

        then

          if (deg_up g_p) = 0 then

(writeln ("try_new_prime_up 2 -----> " ^ PolyML.makestring [1]);

                                   [1]: unipoly

) 

                                                else try_new_prime_up a b d M p g_p p

        else

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

            let 

              val P = P * p

              val g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up 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 (lcoeff_up p1, lcoeff_up p2) \<and> 

               M = LANDAU_MIGNOTTE_bound \<and> p = prime  \<and>  p ~| d  \<and>

               p1 is primitive  \<and>  p2 is primitive

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

    argumentns "a b d M p" named according to [1] p.93 *)

  fun HENSEL_lifting_up a b d M p =

(writeln ("HENSEL_lifting_up: a = " ^ PolyML.makestring a ^

", b = " ^ PolyML.makestring b ^

", d = " ^ PolyML.makestring d ^

", M = " ^ PolyML.makestring M ^

", p = " ^ PolyML.makestring p);

    let

      val p = p next_prime_not_dvd d 

      val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p 

        (*.. nesting of functions see above*)

    in

      if deg_up g_p = 0 then 

(writeln ("HENSEL_lifting_up 1 ----> " ^ PolyML.makestring [1]);

                             [1] 

)

                                else

        let 

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

        in

          if (g %|% a) andalso (g %|% b) then 

(writeln ("HENSEL_lifting_up 2 -----> " ^ PolyML.makestring g);

                                              g:unipoly

) 

                                                        else HENSEL_lifting_up a b d M p

        end

    end

)

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

     gcd_up 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_up a b =

    let val d = abs (Integer.gcd (lcoeff_up a) (lcoeff_up b))

    in 

      if a = b then a else

        HENSEL_lifting_up a b d (2 * d * LANDAU_MIGNOTTE_bound a b) 1

    end;

*}

section {* gcd for multivariate polynomials *}

ML {*

type monom = (int * int list);

type poly = monom list;

*}

subsection {* auxiliary functions *}

ML {*

  fun land a b = a andalso b

  fun all_geq a b = fold land (map2 (curry op >=) a b) true

  fun all_geq' i es = fold land (map ((curry op <=) i) es) true

  fun maxs is = fold Integer.max is (hd is);

  fun div2_swapped d i = i div2 d

  fun nth_swapped i ls = nth ls i;

  fun drop_last ls = (rev o tl o rev) ls

  fun nth_ins n i xs = (*analogous to nth_drop: why not in library.ML?*)

    List.take (xs, n) @ [i] @ List.drop (xs, n);

*}

subsection {* basic notions for multivariate polynomials *}

ML {*

  fun lcoeff (p: poly) = (fst o hd o rev) p

  fun lexp (p: poly) = (snd o hd o rev) p

  fun lmonom (p: poly) = (hd o rev) p

  fun zero_poly (nro_vars: int) = [(0, replicate nro_vars 0)] : poly;

  (* drop leading coefficients equal 0 TODO: signature?*)

  fun drop_lc0 (p: poly as [(c, es)]) = if c = 0 then zero_poly (length es) else p

    | drop_lc0 p = if lcoeff p = 0 then drop_lc0 (drop_last p) else p

  (* gcd of coefficients WN find right location in file *)

  fun gcd_coeff (up: unipoly) = abs (Integer.gcds up)

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

    add_variable p 2  =>  [(3,[1,0,0]),(4,[0,0,1])] with x<y<z *)

  fun add_variable pos (p: poly) = map (apsnd (nth_ins pos 0)) p : poly

*}

subsection {* ordering and other auxiliary funcions *}

ML {*

  (* monomial order: lexicographic order on exponents *)

  fun  lex_ord ((_, a): monom) ((_, b): monom) =

    let val d = drop_tl_zeros (a %-% b)

    in

      if d = [] then false else last d > 0

    end

  fun  lex_ord' ((_, a): monom, (_, b): monom) =

    let val d = drop_tl_zeros (a %-% b)

    in

      if  d = [] then EQUAL 

      else if last d > 0 then GREATER else LESS

    end

  (* add monomials in ordered polynomal *)

  fun add_monoms (p: poly) = 

    let 

      fun add [(0, es)] [] = zero_poly (length es)

        | add [(0, _)] (new: poly) = new

        | add [m] (new: poly) = new @ [m]

        | add ((c1, es1) :: (c2, es2) :: ms) (new: poly) =

          if es1 = es2

          then add ([(c1 + c2, es1)] @ ms) new 

          else

            if c1 = 0 

            then add ((c2, es2) :: ms) new 

            else add ((c2, es2) :: ms) (new @ [(c1, es1)])

    in add p [] : poly end

  (* brings the polynom in correct order and adds monomials *)

  fun order (p: poly) = (add_monoms o (sort lex_ord')) p;

  (* largest exponent of variable at location var *)

  fun max_deg_var (p: poly) var = (maxs o (map ((nth_swapped var) o snd))) p;

*}

subsection {* evaluation, translation uni--multivariate, etc *}

ML {*

  (* evaluate a polynomial in a certain variable and remove this var from the exp-list *)         

  fun eval_monom var value ((c, es): monom) =

    (c * (Integer.pow (nth es var) value), nth_drop var es) : monom

  fun eval (p: poly) var value = order (map (eval_monom var value) p)

  (* transform a multivariate to a univariate polynomial  TODO map?*)         

  fun multi_to_uni (p as (_, es)::_ : poly) = 

    if length es = 1 

    then 

      let fun transfer ([]: poly) (up: unipoly) = up 

            | transfer (mp as (c, es)::ps) up =

                if length up = hd es

                then transfer ps (up @ [c])

                else transfer mp (up @ [0])    

         in transfer (order p) [] end

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

  (* transform a univariate to a multivariate polynomial TODO map *)         

  fun uni_to_multi (up: unipoly) =

    let fun transfer ([]: unipoly) (mp: poly) (_: int) = add_monoms mp

          | transfer up mp var = 

            transfer (tl up) (mp @ [(hd up, [var])]) (var + 1)

    in  transfer up [] 0 end

  (* multiply a poly with a variable represented by a unipoly (at position var) 

    e.g new z: (3x + 2y)*(z - 4)  TODO mp %%*%% (uni_to_multi up) !TODO equal es ! *)

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

    | mult_with_new_var mp up var = 

    let 

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

        | mult mp up var new _ =

          let

            (* the inner loop *)

            fun mult' (_: poly) ([]: unipoly) (_) (new': poly) (_: int) = order new'

              | mult' (mp' as (c, es)::ps) up' var new' c2' =

                let val (first, last) = chop var es

                in mult' mp' (tl up') var

                  (new' @ [(c * hd up', first @ [c2'] @ last)]) (c2' + 1)

                end

          in mult (tl mp) up var (new @ (mult' mp up var [] 0)) (0) end

    in mult mp up var [] 0 : poly end

  fun greater_deg (es1: int list) (es2: int list) =

    if es1 = [] andalso es2 =[] then 0

    else if (nth es1 (length es1 -1)) = (nth es2 (length es1 -1) ) 

         then greater_deg (nth_drop (length es1 -1) es1) (nth_drop (length es1 -1) es2)

         else if (nth es1 (length es1 -1)) > (nth es2 (length es1 -1)) 

              then 1 else 2

  (* Winkler p.95 *) 

  fun find_new_r (p1 as (_, es1)::_ : poly) (p2 as (_, es2)::_ : poly) old =

    let val p1' as (_, es1')::_ = eval p1 (length es1 - 2) (old + 1)

        val p2' as (_, es2')::_ = eval p2 (length es2 - 2) (old + 1);

    in

      if max_deg_var p1' (length es1' - 1) = max_deg_var p1 (length es1 - 1) orelse

         max_deg_var p2' (length es2' - 1) = max_deg_var p2 (length es1 - 1) 

      then old + 1

      else find_new_r p1 p2 (old + 1)

    end

*}

subsection {* addition, multiplication, division *}

ML {*

  (* resulting 0 coefficients are removed by add_monoms *)

  infix %%/

  fun (p: poly) %%/ d = (add_monoms o (map (apfst (div2_swapped d)))) p;

  infix %%*

  fun (p: poly) %%* i = (add_monoms o (map (apfst (Integer.mult i)))) p;

  infix %%+%%

  (* the same algorithms as done by hand

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

    | p1 %%+%% [] = p1

    | p1 %%+%% p2 =

      let 

        fun plus [] p2 new = order (new @ p2)

          | plus p1 [] new = order (new @ p1)

          | plus (p1 as (c1, es1)::ps1) (p2 as (c2, es2)::ps2) new = 

            if greater_deg es1 es2 = 0

            then plus ps1 ps2 (new @ [((c1 + c2), es1)])

            else 

              if greater_deg es1 es2 = 1

              then plus p1 ps2 (new @ [(c2, es2)])

              else plus ps1 p2 (new @ [(c1, es1)])

      in plus p1 p2 [] end

  *)

  fun (p1 : poly) %%+%% (p2 : poly) = order (p1 @ p2)

  infix %%-%%

  fun p1 %%-%% (p2: poly) = p1 %%+%% (map (apfst (Integer.mult ~1)) p2)

  infix %%*%%

  fun mult_monom ((c1, es1): monom) ((c2, es2): monom) =

    (c1 * c2, map2 Integer.add es1 es2): monom;

  fun mult_monoms (ms: poly) m = map (mult_monom m) ms : poly

  fun p1 %%*%% (p2 : poly) = order (fold (curry op @) (map (mult_monoms p1) p2) [])

  infix %%|%% 

  (* the same algorithms as done by hand; TODO fold / map ? *)

  fun ([(c1, es1)]: poly) %%|%%  ([(c2, es2)]: poly) = all_geq es2 es1 andalso c2 mod c1 = 0

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

    | p1 %%|%% p2 =

      if [lmonom p1] %%|%% [lmonom p2]

      then 

        p1 %%|%% 

          (p2 %%-%%

            ([((lcoeff p2) div2 (lcoeff p1), (lexp p2) %-% (lexp p1))]

              %%*%% p1)) 

      else false

  (* %%/%% :: poly \<Rightarrow> poly \<Rightarrow> (poly \<times> poly)

     a %%/%% b = (c, r)

     assumes b \<noteq> [] \<and> 

     yields THE c r. c *\<^isub>p b  +\<^isub>p  r = a *)

  infix %%/%% 

  fun (p1: poly) %%/%% (p2: poly) =

    let 

      fun div_up p1 p2 quot_old = 

      let 

        val p1 as (c, es)::_ = drop_lc0 (order p1);

        val p2 = drop_lc0 (order p2);

        val [c] = zero_poly (length es);

        val d = (replicate (List.length p1 - List.length p2) c) @ p2 ;

        val quot = [(lcoeff p1 div2 lcoeff d, lexp p1 %-% lexp p2)];

        val rest = drop_lc0 (p1 %%-%% (d %%*%% quot));

      in 

        if rest <> [c] andalso all_geq' 0 (lexp rest %-% lexp p2)

        then div_up rest p2 (quot @ quot_old)

        else (quot @ quot_old : poly, rest)

      end

    in div_up p1 p2 [] end 

*}

subsection {* Newton interpolation *}

ML{* (* first step *)

  fun newton_first (x: int list) (f: poly list) (pos: int) =

    let 

      val new_vals = [(hd x) * ~1, 1];

      val p = (add_variable pos (hd f)) %%+%% 

        ((mult_with_new_var 

          (((nth f 1) %%-%% (hd f)) %%/ (nth x 1) - (hd x)))   new_vals pos);

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

    in (p, new_vals, steps) end

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

     | next_step steps new_steps x =  

       next_step (tl steps)

         (new_steps @ 

           [((last new_steps) %%-%% (hd steps)) %%/ ((last x) - (nth x (length x - 3)))])

         (nth_drop (length x -2) x)

  fun NEWTON x f (steps: poly list) (up: unipoly) (p: poly) pos = 

    if length x = 2

    then 

      let val (p', new_vals, steps) = newton_first x [(hd f), (nth f 1)] pos

      in (p', new_vals, steps) end

    else 

      let 

        val new_vals = 

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

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

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

        val steps = next_step steps new_steps x;

        val p' = p %%+%% (mult_with_new_var (last steps) new_vals pos);

      in (order p', new_vals, steps) end 

*}

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

ML {*

  (* naming of n, M, m, r, ...  follows Winkler p. 95 

    n: amount of variables 

    r: *)

  fun gcd_poly' (a as (_, es1)::_ : poly) (b as (_, es2)::_ : poly) n r =

    if lex_ord (lmonom b) (lmonom a) then gcd_poly' b a n r else

      if n > 1

      then

        let val M = 1 + Int.min (max_deg_var a (length es1 - 2), max_deg_var b (length es2 - 2)); 

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

      else 

        let val (a', b') = (multi_to_uni a, multi_to_uni b)

        in

          (* a b are made primitive and this is compensated by multiplying with gcd of divisors *)

          uni_to_multi ((gcd_up (primitive_up a') (primitive_up b')) %* 

                        (abs (Integer.gcd (gcd_coeff a') (gcd_coeff b'))))

        end

  (* fn: poly -> poly -> int -> 

                      int -> int -> int list -> poly list -> poly*)

  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 (eval a (n - 2) r)) 

                            (order (eval 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 ( zero_poly n ) [1] end

  (* fn: poly -> poly -> int -> 

                           int -> int -> int -> int list -> poly list -> 

                           |                  poly -> poly -> int list -> poly*)

  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 (eval a (n - 2) r)) 

                            (order (eval b (n - 2) r)) (n - 1) 0

      in  

        if lex_ord (lmonom g) (lmonom g_r)

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

        else if (lexp g) = (lexp g_r) 

          then

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

            in 

              if (nth steps (length steps - 1)) = (zero_poly (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 (* \<not> lex_ord (lmonom g) (lmonom g_r) *)

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

      end

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

     gcd_poly a b = ((a', b'), c)

     assumes a \<noteq> [] \<and> b \<noteq> [] \<and> 

     yields a = a' *\<^isub>p c  \<and>  b = b' *\<^isub>p c \<and> 

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

  fun gcd_poly (a as (_, es)::_ : poly) b = 

    let val c = gcd_poly' a b (length es) 0

        val (a': poly, _) = a %%/%% c

        val (b': poly, _) = b %%/%% c

    in ((a', b'), c) end

*}

section {* example from the last mail *}

ML {* (* test, see text below *)

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

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

  val ((a', b'), c) = gcd_poly a b;

  a' = [(~1, [1, 0]), (~1, [0, 1])];

  b' = [(~1, [1, 0])];

  c = [(~1, [1, 0]), (1, [0, 1])];

  (* checking the postcondition: *)

  a = (a' %%*%% c) andalso b = (b' %%*%% c); (* where "%%*%%" is "*\<^isub>p" *)

  (* \<forall>c'. (c' dvd\<^isub>p  a \<and> c' dvd\<^isub>p b) \<longrightarrow> c' \<le>\<^isub>p c) can NOT be checked this way, of course *)

*}

text {* last mail to Tobias Nipkow:

Eine erste Idee, wie die Integration der Diplomarbeit für einen Benutzer 

von Isabelle aussehen könnte, wäre zum Beispiel im

   lemma cancel:

     assumes asm3: "x^2 - x*y \<noteq> 0" and asm4: "x \<noteq> 0"

     shows "(x^2 - y^2) / (x^2 - x*y) = (x + y) / (x::real)"

     apply (insert asm3 asm4)

     apply simp

   sorry

die Assumptions

   asm1: "(x^2 - y^2) = (x + y) * (x - y)" and asm2: "x^2 - x*y = x * (x - y)"

im Hintergrund automatisch zu erzeugen (mit der Garantie, dass "(x - y)" der GCD ist) 

und sie dem Simplifier (für die Rule nonzero_mult_divide_mult_cancel_right) zur 

Verfügung zu stellen, sodass anstelle von "sorry" einfach "done" stehen kann.

Und weiters wäre eventuell asm3 zu "x - y \<noteq> 0" zu vereinfachen, eine 

Rewriteorder zum Herstellen einer Normalform festzulegen, etc.

*}

end