(* 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_ML 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 remd = drop_lc0_up (ps %-% (d000 %* quot));

      in

        if remd = [] then true else

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

          then ds %|% remd

          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 remd = drop_lc0_up ((p1m %-% (p2n %* quot)) mod_up m)

    in

      if remd = [] then p2 else

        if length remd < List.length p2

        then mod_up_gcd p2 remd m 

        else mod_up_gcd remd 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

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

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

        (*.. nesting of functions see above*)

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

        iterate_CHINESE_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);

  fun coeffs (p : poly) = map fst p

*}

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)

     p1 %%/%% p2 = (q, r)

     assumes p2 \<noteq> []  \<and>  deg p1 \<ge> deg p2

     yields THE q r. p1 = q *\<^isub>p p2  +\<^isub>p  r *)

  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 remd = drop_lc0 (p1 %%-%% (d %%*%% quot));

      in 

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

        then div_up remd p2 (quot @ quot_old)

        else (quot @ quot_old : poly, remd)

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

  fun gcd_monom (c1, es1) (c2, es2) = (Integer.gcd c1 c2, map2 Integer.min es1 es2) 

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

    assumes: n = length es1 = length es2

    r: *)

  fun gcd_poly' a [monom as (c, es)] _ _ = [fold gcd_monom a monom] 

    | gcd_poly' [monom as (c, es)] b _ _ = [fold gcd_monom b monom] 

    | 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

  (* try_new_r :: poly -> poly -> int -> int -> int -> int list -> poly list -> poly

     try_new_r    a       b       n      M      r      r_list      steps = 

     assumes length a > 1  \<and>  length b > 1

     yields TODO

  *)

  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 iterate_CHINESE a b n M 1 r r_list steps g_r (zero_poly n) [1] end

  (* iterate_CHINESE :: poly -> poly -> int -> int -> int -> int -> int list -> 

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

     iterate_CHINESE    a    |  b       n      M      m      r      r_list 

                             steps        g       g_n     mult =

     assumes length a > 1  \<and>  length b > 1

     yields TODO

  *)

  and iterate_CHINESE 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 iterate_CHINESE 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 iterate_CHINESE a b n M (M + 1) r r_list steps g_r g_n new

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

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

      end

  fun majority_of_coefficients_is_negative_in a b c =

    let val cs = (coeffs a) @ (coeffs b) @ (coeffs c)

    in length (filter (curry op < 0) cs) < length cs div 2 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)  \<and>

       \<not> gcds a' < 0 \<and> gcds a' < 0                                                      *)

  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 

      if majority_of_coefficients_is_negative_in a' b' c

      then ((a' %%* ~1, b' %%* ~1), c %%* ~1) (*makes results look nicer*)

      else ((a', b'), c)

    end

*}

section {* example from the last mail to Tobias Nipkow *}

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

  val a = [(1,[2, 0]), (~1,[0, 2])]; (* x^2 - y^2 *)

  val b = [(1,[2, 0]), (~1,[1, 1])]; (* x^2 - x*y *)

  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.

*}

section {* Preparations for Lucas-Interpretation *}

text {*

  Lucas-Interpretation shall make "gcd_poly" transparent to the user,

  such that step-wise interpretation gives some intuition of the algorithm

  and opportunities for interactive learning by trial and error.

  The vision is to provide a learning environment similar to a good

  chess program, where the matter (chess) is selfcontained and

  learning happens by interaction with the system.

  GDC_Poly gives raise to a major case study on the above aims.

  Below pedagogical questions come first and respective technical

  requirements subsequently.

*}

subsection {* Questions of learners *}

text {*

  The earlier the phase in appropaching a new topic for a student, the more

  divergent the questions are. So here are just some which suggest themselves.

*}

subsubsection {* Why does naive transfer from Z to polynomials not work? *}

text {*

  Given the Eucliden Algorithm in Z ...

*}

ML {*

  fun euclid_int a 0 = a

    | euclid_int a b =

      if a < b then euclid_int b a else

      let

        val (q, r) = (a div b, a mod b)

        val _ = writeln (PolyML.makestring a ^ " = "  ^ PolyML.makestring q  ^ 

          " * "  ^ PolyML.makestring b ^ " + "  ^ PolyML.makestring r)

      in euclid_int b r : int end;

  euclid_int 192 162; (* ..stepwise execution is essential for comprehension *)

*}

text {* ... let it naively transfer to the polynomial ring Z[x]: *}

ML {*

  fun deg p = max_deg_var (p: poly) 0 (*TODO?*)

  fun EUCLID_naive a [(0, [0])] = a   (*TODO: for Z[x,y] we would need [(0, [0, 0])] here*)

    | EUCLID_naive a b =

      if deg a < deg b then EUCLID_naive b a else

      let

        val (q, r) = a %%/%% b

        val _ = writeln (PolyML.makestring a ^ "  ==  "  ^ PolyML.makestring q  ^ 

          "  **  "  ^ PolyML.makestring b ^ "  ++  "  ^ PolyML.makestring r)

      in EUCLID_naive b r : poly end;

*} 

ML {*  

  (*The example used in the inital mail to Tobias Nipkow, simplified to univariate*)

  "--------------------------------";

  val a = uni_to_multi [0,~1,1]; (* -x + x^2 =  x   *(-1 + x) *)

  val b = uni_to_multi [~1,0,1]; (* -1 + x^2 = (1+x)*(-1 + x) *)

  "--------------------------------";

  EUCLID_naive a b = [(1, [0]), (~1, [1])];          (* 1 - x  .. is acceptable*)

*}

ML {*

  "--------------------------------";

  val a = uni_to_multi [~1,0,0,1]; (* -1 + x^3 = (1+x+x^2)*(-1 + x) *)

  val b = uni_to_multi [0,~1,1];   (* -x + x^2 = x        *(-1 + x) *)

  "--------------------------------";

  EUCLID_naive a b = [(~1, [0]), (1, [1])];          (*-1 + x *) 

*}

text {*

  EUCLID_naive, however, does not work in general, because division is 

  not as general as in Z:

*}

ML {*

"--------------------------------";

val a = uni_to_multi [0,~1,0,0,1]; (* -x + x^4 = x*(1+x+x^2)*(-1 + x) *)

val b = uni_to_multi [~2,2];       (* -x + x^2 =  2         *(-1 + x) *)

"--------------------------------";

(*EUCLID_naive a b = [(~1, [0]), (1, [1])] happens to loop*)

*} 

text {*

(*  (x^4 - x) : (2*x - 2) = 1/2 * x^4 ... 

  This division, equal to above division, is impossible in Z[x], rather, it requires Q[x]. 

  This problem of division with integer coefficients proceeds down from the least 

  coefficient to the other coefficients.

  Nevertheless, below we are going to generalise division such that the

  Euclidean Algorithm can work on Z[x].

*}

subsubsection {* Can division of polynomials be adapted to Euclid? *}

text {*

  In principle it is simplest to go from Z[x] to Q[x], polynamials over rational numbers

  and thus, for instance, allow 

            (x^4 - x) : (2*x - 2) = 1/2 * x^4 ... 

  However, Q is computationally less efficient and the terms are more complicated.

  Fortunately, [1].p.83/84 tells, that division in Z[x] can be generalised for the 

  Euclidean Algorithm.

  In order to comprehend stepwise execution of the division algorithm below

  (which is exactly as done by hand) we adapt to our representation of (univariate)

  polynomials as lists and take an instructuve example:

    (2*x^4 + 2*x^3 + 2*x^2 + 2*x + 2) : (3*x + 4) = 2/3 * x^2 ...

  is translated to a reversed arrangement of monomials

    (2 + 2*x + 2*x^2 + 2*x^3 + 2*x^4) : (1 + 2*x + 3*x^2) = ...

  and further to lists:

        [  2,   2,   2,   2,   2] : [1, 2, 3] =

      3*[  2,   2,   2,   2,   2]                   *3

      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:

        [  6,   6,   6,   6,   6]                   :

      --[  0,   0,   2,   4,   6] <--........ * [ 2]

      ------------------------                      :

        [  6,   6,   4,   2,   0]             = [ 2]*3

      3*[  6,   6,   4,   2]                           *3

      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:

        [ 18,  18,  12,   6]                           :

      --[  0,   2,   4,   6] <-------........--- * [2]

      ------------------------                         :

        [ 18,  16,   8,   0]                  = [ 6, 2]

      3*[ 18,  16,   8]                                   *3

      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:

        [ 54,  48,  24]                                    :

      --[  8,  16,  24] <------------........-------- * [8]

      ------------------------                             :

        [46, 32, 0]                           = [18, 6,  8]

      =======================================

     27*[  2,   2,   2,   2,   2] : [1, 2, 3] = [18, 6,  8],  remd [46, 32]

  We implement this algorithm below and check it with these and other values.

  First we need some auxiliary functions.

*}

ML {* (* auxiliary functions *)

  fun for_quot_regarding p1 p2 p1' quot remd =

    let 

       val zeros = deg_up p1' - deg_up remd - 1(*length [q]*)

       val max_qot_deg = deg_up p1 - deg_up p2

     in

       if zeros + 1(*length [q]*) + deg_up quot > max_qot_deg (*possible in last step of div*)

       then max_qot_deg - deg_up quot - 1(*length [q]*) else zeros

     end

  (* multiply polynomial p such that it has degree d with d = deg p*)

  fun mult_to_deg d p =

    let val d' = deg_up p

    in if d >= d' then (replicate (d - d') 0) @ p : unipoly else p end;

  (* find a factor n and a quotient q such that n > 0 \<and> n*c1 = q*c2 *)

  fun fact_quot c1 c2 =

    let

      val d = abs (Integer.gcd c1 c2)

      val n = c2 div2 d

      val q = c1 div2 d

    in if n > 0 then (n, q) else (~1 * n, ~1 * q) end;

  infix %+%

  fun p1 %+% p2 =

    if length p1 > length p2

    then map2 (curry op +) p1 (p2 @ (replicate (List.length p1 - List.length p2) 0))

    else map2 (curry op +) (p1 @ (replicate (List.length p2 - List.length p1) 0)) p2;

  infix %*%

  fun p1 %*% p2 =multi_to_uni (uni_to_multi p1 %%*%% uni_to_multi p2);

  ([1,2,1] %*% [1,3,3,1]) = [1, 5, 10, 10, 5, 1];

  replicate 3 0 = [0, 0, 0];

  replicate 0 0 = [];

  val trace_div = Unsynchronized.ref true;

  val trace_div_invariant = Unsynchronized.ref false; (*.. true expects !trace_div = false *)

  fun div_invariant2 p1 p2 n ns zeros q quot remd = 

    (PolyML.makestring p1 ^ " * "  ^ PolyML.makestring (n * ns) ^ " = " ^ 

     PolyML.makestring (mult_to_deg (deg_up p1 - deg_up p2) 

      ((zeros @ [q] @ (quot %* n)))) ^ 

     " ** " ^ PolyML.makestring p2 ^ " ++ " ^ 

     PolyML.makestring ((rev o (mult_to_deg (deg_up p1)) o rev) remd))

  fun writeln_trace p1 p1' p2 quot n q (zeros:int) remd (ns: int) = 

    (if (! trace_div) then

      ( writeln ("  div   " ^ PolyML.makestring p1' ^ " * " ^ PolyML.makestring n);

        writeln  "  div   ~~~~~~~~~~~~~~~~~~~~~~~~~~~";

        writeln ("  div   " ^ PolyML.makestring (p1' %* n));

        writeln ("  div-- " ^ PolyML.makestring (mult_to_deg (deg_up p1') p2 %* q) ^ 

          " <--- " ^ PolyML.makestring p2 ^ " * " ^ PolyML.makestring q);

        writeln  "  div   ---------------------------";

        writeln ("  div   " ^ PolyML.makestring ((p1' %* n)%-%(mult_to_deg (deg_up p1') p2 %* q)) ^ 

          "                            quot = " ^ 

          PolyML.makestring (replicate zeros 0) ^ " @ [" ^ 

          PolyML.makestring q ^ "] @ (" ^ PolyML.makestring n ^ 

          " * " ^ PolyML.makestring quot ^ ")" );

        if (p1' %* n) = ((mult_to_deg (deg_up p1') (p2 %* q)) %+% remd) then () 

        else error ("invariant 1 does not hold: " ^ PolyML.makestring (p1' %* n) ^ " = " ^ 

        PolyML.makestring (mult_to_deg (deg_up p1') p2) ^ " * " ^ PolyML.makestring q ^ " ++ " ^ 

        PolyML.makestring (remd @ []))

      ) else ();

      (*invariant 2: initial p1 *\<^isub>s ns = quot *\<^isub>p p2   +\<^isub>p  remd *)

      if (p1 %* (n * ns)) = 

        ((mult_to_deg (deg_up p1 - deg_up p2) 

          ((replicate zeros 0) @ [q] @ (quot %* n)) %*% p2 %+% remd))

      then

        if (! trace_div_invariant) 

        then writeln ("  div   " ^ div_invariant2 p1 p2 n ns (replicate zeros 0) q quot remd)

        else ()

      else error ("invariant 2 does not hold: " ^ div_invariant2 p1 p2 n ns (replicate zeros 0) q quot remd)

    )

  (* generalized division: divide in Z[x] with multiplication of the dividend 

     such that the quotient is again in Z[x], i.e. has integer coefficients.

     fact_div :: unipoly \<Rightarrow> unipoly \<Rightarrow> (unipoly, unipoly, nat)

     fact_div p1 p2 = (n, q, r)

       assumes p2 \<noteq> [0]  \<and>  deg p1 \<ge> deg p2

       yields p1 *\<^isub>s n = q *\<^isub>p p2   +\<^isub>p   r 

     note 1: p1' changes with iterations, p1 remains equal -- only required for invariant 2 

     note 2: test -- fun %*/% Subscript raised (line 509 of lib-- shows a bug. *)

  infix %*/%

  fun div_int_up p1 p1' p2 quot ns =