(* 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 this been removed since 2002 ? *)

   fun last_elem 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;

 *}

 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 * 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_elem 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_elem 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 =

     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] = (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 =

     if P > M then 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 [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 = 

     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 = 0 then [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 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 cero_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 cero_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 (p: poly) var = map (apsnd (nth_ins var 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_elem 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_elem d > 0 then GREATER else LESS

     end

   (* add monomials in ordered polynomal *)

   fun add_monoms (p: poly) = 

     let 

       fun add [(0, es)] [] = cero_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 from the exp-list *)         

   fun eval (p: poly) var value = 

     let 

       fun eval [] _ _ new = order new 

         | eval ((c, es)::ps) var value new =

           eval ps var value (new @ [((Integer.pow  (nth es var) value) * c, nth_drop var es)]);

     in eval p var value [] end

   (* transform a multivariate to a univariate polynomial *)         

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

   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 fold / map ? *)

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

   (* the same algorithms as done by hand

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

     let 

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

         | mult p1 [] regular_p2 new = mult (tl p1) regular_p2 regular_p2 new

         | mult (p1 as (c1, es1)::_) ((c2, es2)::ps2) regular_p2 (new: poly) = 

             mult p1 ps2 regular_p2 (new @ [(c1 * c2, map2 Integer.add es1 es2)])

     in 

       if length p1 > length p2 

       then mult p1 p2 p2 [] 

       else mult p2 p1 p1 [] 

     end

   *)

   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] = cero_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) (ord: int) =

     let 

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

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

         ((mult_with_new_var 

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

       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_elem new_steps) %%-%% (hd steps)) %%/ ((last_elem x) - (nth x (length x - 3)))])

          (nth_drop (length x -2) x)

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

     if length x = 2

     then 

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

       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_elem x) - (nth x (length x - 2)))];

         val steps = next_step steps new_steps x;

         val p' = p %%+%% (mult_with_new_var (last_elem steps) new_vals ord);

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

   fun gcd_poly' (a as (c1, es1)::ps1 : poly) (b as (c2, es2)::ps2 : poly) n r =

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

       if n - 1 = 0 

       then

         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

       else 

         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

   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 ( cero_poly 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 (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)) = (cero_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 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 = (a' %%*%% c) andalso b = (b' %%*%% c); (* where "%%*%%" is "*\<^isub>p" *)

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

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

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

 *}

 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